diff --git a/Exec/DevTests/MovingTerrain/prob.H b/Exec/DevTests/MovingTerrain/prob.H index c56b3bb7b..6bf576751 100644 --- a/Exec/DevTests/MovingTerrain/prob.H +++ b/Exec/DevTests/MovingTerrain/prob.H @@ -44,11 +44,7 @@ public: const amrex::Real& time) override; void erf_init_rayleigh ( - amrex::Vector& tau, - amrex::Vector& ubar, - amrex::Vector& vbar, - amrex::Vector& wbar, - amrex::Vector& thetabar, + amrex::Vector >& rayleigh_ptrs, amrex::Geometry const& geom, std::unique_ptr& z_phys_cc) override; diff --git a/Exec/DevTests/MovingTerrain/prob.cpp b/Exec/DevTests/MovingTerrain/prob.cpp index dc1859077..bae66ec74 100644 --- a/Exec/DevTests/MovingTerrain/prob.cpp +++ b/Exec/DevTests/MovingTerrain/prob.cpp @@ -132,34 +132,29 @@ Problem::init_custom_pert( void Problem::erf_init_rayleigh( - amrex::Vector& tau, - amrex::Vector& ubar, - amrex::Vector& vbar, - amrex::Vector& wbar, - amrex::Vector& thetabar, - amrex::Geometry const& geom, + amrex::Vector >& rayleigh_ptrs, + amrex::Geometry const& geom, std::unique_ptr& /*z_phys_cc*/) { - // amrex::Error("Should never get here for MovingTerrain problem"); - const int khi = geom.Domain().bigEnd()[2]; - - for (int k = 0; k <= khi; k++) - { - tau[k] = 0.0; - ubar[k] = 0.0; - vbar[k] = 0.0; - wbar[k] = 0.0; - thetabar[k] = 300.0; - } + const int khi = geom.Domain().bigEnd()[2]; - // Damping above k = 60 - for (int k = 60; k <= khi; k++) - { - tau[k] = 10.; // Remember that this gets multiplied by dt - ubar[k] = 0.0; - vbar[k] = 0.0; - wbar[k] = 0.0; - thetabar[k] = 300.0; + for (int k = 0; k <= khi; k++) + { + rayleigh_ptrs[Rayleigh::tau][k] = 0.0; + rayleigh_ptrs[Rayleigh::ubar][k] = 0.0; + rayleigh_ptrs[Rayleigh::vbar][k] = 0.0; + rayleigh_ptrs[Rayleigh::wbar][k] = 0.0; + rayleigh_ptrs[Rayleigh::thetabar][k] = 300.0; + } + + // Damping above k = 60 + for (int k = 60; k <= khi; k++) + { + rayleigh_ptrs[Rayleigh::tau][k] = 10.0; // Remember that this gets multiplied by dt + rayleigh_ptrs[Rayleigh::ubar][k] = 2.0; + rayleigh_ptrs[Rayleigh::vbar][k] = 1.0; + rayleigh_ptrs[Rayleigh::wbar][k] = 0.0; + rayleigh_ptrs[Rayleigh::thetabar][k] = 300.0; } } diff --git a/Exec/Radiation/prob.H b/Exec/Radiation/prob.H index 5c20a28c2..b7eaa308d 100644 --- a/Exec/Radiation/prob.H +++ b/Exec/Radiation/prob.H @@ -47,11 +47,7 @@ public: const SolverChoice& sc) override; void erf_init_rayleigh ( - amrex::Vector& tau, - amrex::Vector& ubar, - amrex::Vector& vbar, - amrex::Vector& wbar, - amrex::Vector& thetabar, + amrex::Vector >& rayleigh_ptrs, amrex::Geometry const& geom, std::unique_ptr& z_phys_cc) override; diff --git a/Exec/Radiation/prob.cpp b/Exec/Radiation/prob.cpp index 912d82b54..ac3c6f662 100644 --- a/Exec/Radiation/prob.cpp +++ b/Exec/Radiation/prob.cpp @@ -214,23 +214,19 @@ Problem::init_custom_pert( void Problem::erf_init_rayleigh( - amrex::Vector& tau, - amrex::Vector& ubar, - amrex::Vector& vbar, - amrex::Vector& wbar, - amrex::Vector& thetabar, + amrex::Vector >& rayleigh_ptrs, amrex::Geometry const& geom, - std::unique_ptr& /*z_phys_cc*/) + std::unique_ptr& /*z_phys_cc*/) { const int khi = geom.Domain().bigEnd()[2]; // We just use these values to test the Rayleigh damping for (int k = 0; k <= khi; k++) { - tau[k] = 1.0; - ubar[k] = 2.0; - vbar[k] = 1.0; - wbar[k] = 0.0; - thetabar[k] = parms.T_0; + rayleigh_ptrs[Rayleigh::tau][k] = 1.0; + rayleigh_ptrs[Rayleigh::ubar][k] = 2.0; + rayleigh_ptrs[Rayleigh::vbar][k] = 1.0; + rayleigh_ptrs[Rayleigh::wbar][k] = 0.0; + rayleigh_ptrs[Rayleigh::thetabar][k] = parms.T_0; } } diff --git a/Exec/RegTests/ScalarAdvDiff/prob.H b/Exec/RegTests/ScalarAdvDiff/prob.H index 4208fc8d5..c724fc2b3 100644 --- a/Exec/RegTests/ScalarAdvDiff/prob.H +++ b/Exec/RegTests/ScalarAdvDiff/prob.H @@ -53,11 +53,7 @@ public: const SolverChoice& sc) override; void erf_init_rayleigh ( - amrex::Vector& tau, - amrex::Vector& ubar, - amrex::Vector& vbar, - amrex::Vector& wbar, - amrex::Vector& thetabar, + amrex::Vector >& rayleigh_ptrs, amrex::Geometry const& geom, std::unique_ptr& z_phys_cc) override; diff --git a/Exec/RegTests/ScalarAdvDiff/prob.cpp b/Exec/RegTests/ScalarAdvDiff/prob.cpp index 221c3b2ce..ba9cbaf86 100644 --- a/Exec/RegTests/ScalarAdvDiff/prob.cpp +++ b/Exec/RegTests/ScalarAdvDiff/prob.cpp @@ -36,24 +36,20 @@ Problem::Problem() void Problem::erf_init_rayleigh( - amrex::Vector& tau, - amrex::Vector& ubar, - amrex::Vector& vbar, - amrex::Vector& wbar, - amrex::Vector& thetabar, - amrex::Geometry const& geom, - std::unique_ptr& z_phys_cc) + amrex::Vector >& rayleigh_ptrs, + amrex::Geometry const& geom, + std::unique_ptr& /*z_phys_cc*/) { const int khi = geom.Domain().bigEnd()[2]; // We just use these values to test the Rayleigh damping for (int k = 0; k <= khi; k++) { - tau[k] = 1.0; - ubar[k] = 2.0; - vbar[k] = 1.0; - wbar[k] = 0.0; - thetabar[k] = parms.T_0; + rayleigh_ptrs[Rayleigh::tau][k] = 1.0; + rayleigh_ptrs[Rayleigh::ubar][k] = 2.0; + rayleigh_ptrs[Rayleigh::vbar][k] = 1.0; + rayleigh_ptrs[Rayleigh::wbar][k] = 0.0; + rayleigh_ptrs[Rayleigh::thetabar][k] = parms.T_0; } } diff --git a/Exec/SquallLine_2D/prob.H b/Exec/SquallLine_2D/prob.H index 9f2fd7122..9de09946a 100644 --- a/Exec/SquallLine_2D/prob.H +++ b/Exec/SquallLine_2D/prob.H @@ -62,11 +62,7 @@ public: amrex::Geometry const& geom) override; void erf_init_rayleigh ( - amrex::Vector& tau, - amrex::Vector& ubar, - amrex::Vector& vbar, - amrex::Vector& wbar, - amrex::Vector& thetabar, + amrex::Vector >& rayleigh_ptrs, amrex::Geometry const& geom, std::unique_ptr& z_phys_cc) override; diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 8c4ccdf81..7cd456357 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -402,25 +402,20 @@ Problem::init_custom_pert ( } void -Problem::erf_init_rayleigh (amrex::Vector& tau, - amrex::Vector& ubar, - amrex::Vector& vbar, - amrex::Vector& wbar, - amrex::Vector& thetabar, - amrex::Geometry const& geom, - std::unique_ptr& /*z_phys_cc*/) +Problem::erf_init_rayleigh( + amrex::Vector >& rayleigh_ptrs, + amrex::Geometry const& geom, + std::unique_ptr& /*z_phys_cc*/) { const int khi = geom.Domain().bigEnd()[2]; // We just use these values to test the Rayleigh damping for (int k = 0; k <= khi; k++) { - tau[k] = 1.0; - ubar[k] = 2.0; - vbar[k] = 1.0; - wbar[k] = 0.0; - thetabar[k] = parms.Theta_0; + rayleigh_ptrs[Rayleigh::tau][k] = 1.0; + rayleigh_ptrs[Rayleigh::ubar][k] = 2.0; + rayleigh_ptrs[Rayleigh::vbar][k] = 1.0; + rayleigh_ptrs[Rayleigh::wbar][k] = 0.0; + rayleigh_ptrs[Rayleigh::thetabar][k] = parms.T_0; } } - - diff --git a/Exec/SuperCell/prob.H b/Exec/SuperCell/prob.H index 5c20a28c2..b7eaa308d 100644 --- a/Exec/SuperCell/prob.H +++ b/Exec/SuperCell/prob.H @@ -47,11 +47,7 @@ public: const SolverChoice& sc) override; void erf_init_rayleigh ( - amrex::Vector& tau, - amrex::Vector& ubar, - amrex::Vector& vbar, - amrex::Vector& wbar, - amrex::Vector& thetabar, + amrex::Vector >& rayleigh_ptrs, amrex::Geometry const& geom, std::unique_ptr& z_phys_cc) override; diff --git a/Exec/SuperCell/prob.cpp b/Exec/SuperCell/prob.cpp index 128fa7471..67a0668a5 100644 --- a/Exec/SuperCell/prob.cpp +++ b/Exec/SuperCell/prob.cpp @@ -214,11 +214,7 @@ Problem::init_custom_pert( void Problem::erf_init_rayleigh( - amrex::Vector& tau, - amrex::Vector& ubar, - amrex::Vector& vbar, - amrex::Vector& wbar, - amrex::Vector& thetabar, + amrex::Vector >& rayleigh_ptrs, amrex::Geometry const& geom, std::unique_ptr& /*z_phys_cc*/) { @@ -227,10 +223,10 @@ Problem::erf_init_rayleigh( // We just use these values to test the Rayleigh damping for (int k = 0; k <= khi; k++) { - tau[k] = 1.0; - ubar[k] = 2.0; - vbar[k] = 1.0; - wbar[k] = 0.0; - thetabar[k] = parms.T_0; + rayleigh_ptrs[Rayleigh::tau][k] = 1.0; + rayleigh_ptrs[Rayleigh::ubar][k] = 2.0; + rayleigh_ptrs[Rayleigh::vbar][k] = 1.0; + rayleigh_ptrs[Rayleigh::wbar][k] = 0.0; + rayleigh_ptrs[Rayleigh::thetabar][k] = parms.T_0; } } diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 22a8dcfc1..abf5e6b6a 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -36,6 +36,10 @@ enum struct Coord { x, y, z }; +enum Rayleigh { + tau, ubar, vbar, wbar, thetabar, nvars +}; + /** * Container holding many of the algorithmic options and parameters */ diff --git a/Source/Diffusion/Diffusion.H b/Source/Diffusion/Diffusion.H index 93aedd486..042750bc3 100644 --- a/Source/Diffusion/Diffusion.H +++ b/Source/Diffusion/Diffusion.H @@ -18,8 +18,6 @@ void DiffusionSrcForMom_N (const amrex::Box& bxx, const amrex::Box& bxy, const a const amrex::Array4& tau12 , const amrex::Array4& tau13 , const amrex::Array4& tau23 , - const amrex::Array4& cons, - const DiffChoice& diffChoice, const amrex::GpuArray& dxInv, const amrex::Array4& mf_m , const amrex::Array4& mf_u , @@ -34,8 +32,7 @@ void DiffusionSrcForMom_T (const amrex::Box& bxx, const amrex::Box& bxy, const a const amrex::Array4& tau12 , const amrex::Array4& tau13, const amrex::Array4& tau21 , const amrex::Array4& tau23, const amrex::Array4& tau31 , const amrex::Array4& tau32, - const amrex::Array4& cons, const amrex::Array4& detJ, - const DiffChoice& diffChoice, + const amrex::Array4& detJ, const amrex::GpuArray& dxInv, const amrex::Array4& mf_m , const amrex::Array4& mf_u , diff --git a/Source/Diffusion/DiffusionSrcForMom_N.cpp b/Source/Diffusion/DiffusionSrcForMom_N.cpp index 88934468c..a5ece6dfa 100644 --- a/Source/Diffusion/DiffusionSrcForMom_N.cpp +++ b/Source/Diffusion/DiffusionSrcForMom_N.cpp @@ -19,8 +19,6 @@ using namespace amrex; * @param[in] tau12 12 stress * @param[in] tau13 13 stress * @param[in] tau23 23 stress - * @param[in] cons conserved cell center quantities - * @param[in] diffChoice container with diffusion parameters * @param[in] dxInv inverse cell size array * @param[in] mf_m map factor at cell center */ @@ -32,7 +30,6 @@ DiffusionSrcForMom_N (const Box& bxx, const Box& bxy , const Box& bxz, const Array4& tau11, const Array4& tau22, const Array4& tau33, const Array4& tau12, const Array4& tau13, const Array4& tau23, - const Array4& cons , const DiffChoice& diffChoice, const GpuArray& dxInv, const Array4& mf_m, const Array4& /*mf_u*/, diff --git a/Source/Diffusion/DiffusionSrcForMom_T.cpp b/Source/Diffusion/DiffusionSrcForMom_T.cpp index eefc92a85..463b74c72 100644 --- a/Source/Diffusion/DiffusionSrcForMom_T.cpp +++ b/Source/Diffusion/DiffusionSrcForMom_T.cpp @@ -22,7 +22,6 @@ using namespace amrex; * @param[in] tau23 23 stress * @param[in] tau31 31 stress * @param[in] tau32 32 stress - * @param[in] cons conserved cell center quantities * @param[in] detJ Jacobian determinant * @param[in] differChoice container with diffusion parameters * @param[in] dxInv inverse cell size array @@ -37,8 +36,7 @@ DiffusionSrcForMom_T (const Box& bxx, const Box& bxy , const Box& bxz, const Array4& tau12, const Array4& tau13, const Array4& tau21, const Array4& tau23, const Array4& tau31, const Array4& tau32, - const Array4& cons , const Array4& detJ , - const DiffChoice& diffChoice, + const Array4& detJ , const GpuArray& dxInv, const Array4& mf_m, const Array4& /*mf_u*/, diff --git a/Source/ERF.H b/Source/ERF.H index ecb67f764..5152d2f75 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -747,20 +747,14 @@ private: static int ng_pres_hse; // Custom source terms - amrex::Vector > h_rhotheta_src; + amrex::Vector< amrex::Vector > h_rhotheta_src; amrex::Vector > d_rhotheta_src; - // Mean quantities and tau for Rayleigh damping - amrex::Vector > h_rayleigh_tau; - amrex::Vector > h_rayleigh_ubar; - amrex::Vector > h_rayleigh_vbar; - amrex::Vector > h_rayleigh_wbar; - amrex::Vector > h_rayleigh_thetabar; - amrex::Vector > d_rayleigh_tau; - amrex::Vector > d_rayleigh_ubar; - amrex::Vector > d_rayleigh_vbar; - amrex::Vector > d_rayleigh_wbar; - amrex::Vector > d_rayleigh_thetabar; + // This is a vector over levels of vectors across quantities of Vectors + amrex::Vector > > h_rayleigh_ptrs; + + // This is a vector over levels of vectors across quantities of DeviceVectors + amrex::Vector > > d_rayleigh_ptrs; amrex::Vector h_havg_density; amrex::Vector h_havg_temperature; diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index b2543e7b7..11527b531 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -18,46 +18,31 @@ ERF::initRayleigh () { AMREX_ALWAYS_ASSERT(solverChoice.use_rayleigh_damping); - h_rayleigh_tau.resize(max_level+1, amrex::Vector(0)); - h_rayleigh_ubar.resize(max_level+1, amrex::Vector(0)); - h_rayleigh_vbar.resize(max_level+1, amrex::Vector(0)); - h_rayleigh_wbar.resize(max_level+1, amrex::Vector(0)); - h_rayleigh_thetabar.resize(max_level+1, amrex::Vector(0)); - d_rayleigh_tau.resize(max_level+1, amrex::Gpu::DeviceVector(0)); - d_rayleigh_ubar.resize(max_level+1, amrex::Gpu::DeviceVector(0)); - d_rayleigh_vbar.resize(max_level+1, amrex::Gpu::DeviceVector(0)); - d_rayleigh_wbar.resize(max_level+1, amrex::Gpu::DeviceVector(0)); - d_rayleigh_thetabar.resize(max_level+1, amrex::Gpu::DeviceVector(0)); + h_rayleigh_ptrs.resize(max_level+1); + d_rayleigh_ptrs.resize(max_level+1); for (int lev = 0; lev <= finest_level; lev++) { + // These have 5 components: tau, ubar, vbar, wbar, thetabar + h_rayleigh_ptrs[lev].resize(Rayleigh::nvars); + d_rayleigh_ptrs[lev].resize(Rayleigh::nvars); + const int zlen_rayleigh = geom[lev].Domain().length(2); - h_rayleigh_tau[lev].resize(zlen_rayleigh, 0.0_rt); - d_rayleigh_tau[lev].resize(zlen_rayleigh, 0.0_rt); - h_rayleigh_ubar[lev].resize(zlen_rayleigh, 0.0_rt); - d_rayleigh_ubar[lev].resize(zlen_rayleigh, 0.0_rt); - h_rayleigh_vbar[lev].resize(zlen_rayleigh, 0.0_rt); - d_rayleigh_vbar[lev].resize(zlen_rayleigh, 0.0_rt); - h_rayleigh_wbar[lev].resize(zlen_rayleigh, 0.0_rt); - d_rayleigh_wbar[lev].resize(zlen_rayleigh, 0.0_rt); - h_rayleigh_thetabar[lev].resize(zlen_rayleigh, 0.0_rt); - d_rayleigh_thetabar[lev].resize(zlen_rayleigh, 0.0_rt); - - prob->erf_init_rayleigh(h_rayleigh_tau[lev], h_rayleigh_ubar[lev], h_rayleigh_vbar[lev], - h_rayleigh_wbar[lev], h_rayleigh_thetabar[lev], geom[lev], - z_phys_cc[lev]); - // Copy from host version to device version - amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_tau[lev].begin(), h_rayleigh_tau[lev].end(), - d_rayleigh_tau[lev].begin()); - amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_ubar[lev].begin(), h_rayleigh_ubar[lev].end(), - d_rayleigh_ubar[lev].begin()); - amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_vbar[lev].begin(), h_rayleigh_vbar[lev].end(), - d_rayleigh_vbar[lev].begin()); - amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_wbar[lev].begin(), h_rayleigh_wbar[lev].end(), - d_rayleigh_wbar[lev].begin()); - amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_thetabar[lev].begin(), h_rayleigh_thetabar[lev].end(), - d_rayleigh_thetabar[lev].begin()); + // Allocate space for these 1D vectors + for (int n = 0; n < Rayleigh::nvars; n++) { + h_rayleigh_ptrs[lev][n].resize(zlen_rayleigh, 0.0_rt); + d_rayleigh_ptrs[lev][n].resize(zlen_rayleigh, 0.0_rt); + } + + // Init the host vectors + prob->erf_init_rayleigh(h_rayleigh_ptrs[lev], geom[lev], z_phys_cc[lev]); + + // Copy from host vectors to device vectors + for (int n = 0; n < Rayleigh::nvars; n++) { + amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_ptrs[lev][n].begin(), h_rayleigh_ptrs[lev][n].end(), + d_rayleigh_ptrs[lev][n].begin()); + } } } @@ -93,29 +78,29 @@ ERF::setRayleighRefFromSounding (bool restarting) for (int k = 0; k <= khi; k++) { const Real z = prob_lo[2] + (k + 0.5) * dx[2]; - h_rayleigh_ubar[lev][k] = interpolate_1d(z_inp_sound, U_inp_sound, z, inp_sound_size); - h_rayleigh_vbar[lev][k] = interpolate_1d(z_inp_sound, V_inp_sound, z, inp_sound_size); - h_rayleigh_wbar[lev][k] = 0.0; - h_rayleigh_thetabar[lev][k] = interpolate_1d(z_inp_sound, theta_inp_sound, z, inp_sound_size); - if (h_rayleigh_tau[lev][k] > 0) { - amrex::Print() << z << ":" << " tau=" << h_rayleigh_tau[lev][k]; - if (solverChoice.rayleigh_damp_U) amrex::Print() << " ubar=" << h_rayleigh_ubar[lev][k]; - if (solverChoice.rayleigh_damp_V) amrex::Print() << " vbar=" << h_rayleigh_vbar[lev][k]; - if (solverChoice.rayleigh_damp_W) amrex::Print() << " wbar=" << h_rayleigh_wbar[lev][k]; - if (solverChoice.rayleigh_damp_T) amrex::Print() << " thetabar=" << h_rayleigh_thetabar[lev][k]; + h_rayleigh_ptrs[lev][Rayleigh::ubar][k] = interpolate_1d(z_inp_sound, U_inp_sound, z, inp_sound_size); + h_rayleigh_ptrs[lev][Rayleigh::vbar][k] = interpolate_1d(z_inp_sound, V_inp_sound, z, inp_sound_size); + h_rayleigh_ptrs[lev][Rayleigh::wbar][k] = Real(0.0); + h_rayleigh_ptrs[lev][Rayleigh::thetabar][k] = interpolate_1d(z_inp_sound, theta_inp_sound, z, inp_sound_size); + if (h_rayleigh_ptrs[lev][Rayleigh::tau][k] > 0) { + amrex::Print() << z << ":" << " tau=" << h_rayleigh_ptrs[lev][Rayleigh::tau][k]; + if (solverChoice.rayleigh_damp_U) amrex::Print() << " ubar = " << h_rayleigh_ptrs[lev][Rayleigh::ubar][k]; + if (solverChoice.rayleigh_damp_V) amrex::Print() << " vbar = " << h_rayleigh_ptrs[lev][Rayleigh::vbar][k]; + if (solverChoice.rayleigh_damp_W) amrex::Print() << " wbar = " << h_rayleigh_ptrs[lev][Rayleigh::wbar][k]; + if (solverChoice.rayleigh_damp_T) amrex::Print() << " thetabar= " << h_rayleigh_ptrs[lev][Rayleigh::thetabar][k]; amrex::Print() << std::endl; } } // Copy from host version to device version - amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_ubar[lev].begin(), h_rayleigh_ubar[lev].end(), - d_rayleigh_ubar[lev].begin()); - amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_vbar[lev].begin(), h_rayleigh_vbar[lev].end(), - d_rayleigh_vbar[lev].begin()); - amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_wbar[lev].begin(), h_rayleigh_wbar[lev].end(), - d_rayleigh_wbar[lev].begin()); - amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_thetabar[lev].begin(), h_rayleigh_thetabar[lev].end(), - d_rayleigh_thetabar[lev].begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_ptrs[lev][Rayleigh::ubar].begin(), h_rayleigh_ptrs[lev][Rayleigh::ubar].end(), + d_rayleigh_ptrs[lev][Rayleigh::ubar].begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_ptrs[lev][Rayleigh::vbar].begin(), h_rayleigh_ptrs[lev][Rayleigh::vbar].end(), + d_rayleigh_ptrs[lev][Rayleigh::vbar].begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_ptrs[lev][Rayleigh::wbar].begin(), h_rayleigh_ptrs[lev][Rayleigh::wbar].end(), + d_rayleigh_ptrs[lev][Rayleigh::wbar].begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_ptrs[lev][Rayleigh::thetabar].begin(), h_rayleigh_ptrs[lev][Rayleigh::thetabar].end(), + d_rayleigh_ptrs[lev][Rayleigh::thetabar].begin()); } } diff --git a/Source/Prob/init_rayleigh_damping.H b/Source/Prob/init_rayleigh_damping.H index d1fd58afe..110e65deb 100644 --- a/Source/Prob/init_rayleigh_damping.H +++ b/Source/Prob/init_rayleigh_damping.H @@ -3,11 +3,7 @@ * on Durran and Klemp 1983 */ void -erf_init_rayleigh (amrex::Vector& tau, - amrex::Vector& ubar, - amrex::Vector& vbar, - amrex::Vector& wbar, - amrex::Vector& thetabar, +erf_init_rayleigh (amrex::Vector >& rayleigh_ptrs, amrex::Geometry const& geom, std::unique_ptr& z_phys_cc) override { @@ -51,20 +47,20 @@ erf_init_rayleigh (amrex::Vector& tau, if (zfrac >= 0) { const amrex::Real sinefac = std::sin(PIoTwo*zfrac); - tau[k] = parms.dampcoef * sinefac * sinefac; - ubar[k] = parms.U_0; - vbar[k] = parms.V_0; - wbar[k] = parms.W_0; - thetabar[k] = parms.T_0; + rayleigh_ptrs[Rayleigh::tau][k] = parms.dampcoef * sinefac * sinefac; + rayleigh_ptrs[Rayleigh::ubar][k] = parms.U_0; + rayleigh_ptrs[Rayleigh::vbar][k] = parms.V_0; + rayleigh_ptrs[Rayleigh::wbar][k] = parms.W_0; + rayleigh_ptrs[Rayleigh::thetabar][k] = parms.T_0; amrex::Print() << " " << zcc[k] << " " << sinefac*sinefac << std::endl; } else { - tau[k] = 0.0; - ubar[k] = 0.0; - vbar[k] = 0.0; - wbar[k] = 0.0; - thetabar[k] = 0.0; + rayleigh_ptrs[Rayleigh::tau][k] = amrex::Real(0.0); + rayleigh_ptrs[Rayleigh::ubar][k] = amrex::Real(0.0); + rayleigh_ptrs[Rayleigh::vbar][k] = amrex::Real(0.0); + rayleigh_ptrs[Rayleigh::wbar][k] = amrex::Real(0.0); + rayleigh_ptrs[Rayleigh::thetabar][k] = amrex::Real(0.0); } } } diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 06b124ab4..6dd24a47a 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -74,11 +74,13 @@ void ERF::advance_dycore(int level, Real* dptr_rhotheta_src = solverChoice.custom_rhotheta_forcing ? d_rhotheta_src[level].data() : nullptr; - Real* dptr_rayleigh_tau = solverChoice.use_rayleigh_damping ? d_rayleigh_tau[level].data() : nullptr; - Real* dptr_rayleigh_ubar = solverChoice.use_rayleigh_damping ? d_rayleigh_ubar[level].data() : nullptr; - Real* dptr_rayleigh_vbar = solverChoice.use_rayleigh_damping ? d_rayleigh_vbar[level].data() : nullptr; - Real* dptr_rayleigh_wbar = solverChoice.use_rayleigh_damping ? d_rayleigh_wbar[level].data() : nullptr; - Real* dptr_rayleigh_thetabar = solverChoice.use_rayleigh_damping ? d_rayleigh_thetabar[level].data() : nullptr; + Vector d_rayleigh_ptrs_at_lev; + d_rayleigh_ptrs_at_lev.resize(Rayleigh::nvars); + d_rayleigh_ptrs_at_lev[Rayleigh::tau] = solverChoice.use_rayleigh_damping ? d_rayleigh_ptrs[level][Rayleigh::tau ].data() : nullptr; + d_rayleigh_ptrs_at_lev[Rayleigh::ubar] = solverChoice.use_rayleigh_damping ? d_rayleigh_ptrs[level][Rayleigh::ubar].data() : nullptr; + d_rayleigh_ptrs_at_lev[Rayleigh::vbar] = solverChoice.use_rayleigh_damping ? d_rayleigh_ptrs[level][Rayleigh::vbar].data() : nullptr; + d_rayleigh_ptrs_at_lev[Rayleigh::wbar] = solverChoice.use_rayleigh_damping ? d_rayleigh_ptrs[level][Rayleigh::wbar].data() : nullptr; + d_rayleigh_ptrs_at_lev[Rayleigh::thetabar] = solverChoice.use_rayleigh_damping ? d_rayleigh_ptrs[level][Rayleigh::thetabar].data() : nullptr; bool l_use_terrain = solverChoice.use_terrain; bool l_use_diff = ( (dc.molec_diff_type != MolecDiffType::None) || diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index 2135af445..a46d2c36a 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -62,11 +62,7 @@ using namespace amrex; * @param[in] mapfac_u map factor at x-faces * @param[in] mapfac_v map factor at y-faces * @param[in] dptr_rhotheta_src custom temperature source term - * @param[in] dptr_rayleigh_tau strength of Rayleigh damping - * @param[in] dptr_rayleigh_ubar reference value for x-velocity used to define Rayleigh damping - * @param[in] dptr_rayleigh_vbar reference value for y-velocity used to define Rayleigh damping - * @param[in] dptr_rayleigh_wbar reference value for z-velocity used to define Rayleigh damping - * @param[in] dptr_rayleigh_thetabar reference value for potential temperature used to define Rayleigh damping + * @param[in] d_rayleigh_ptrs_at_lev Vector of {strength of Rayleigh damping, reference value for xvel/yvel/zvel/theta} used to define Rayleigh damping */ void erf_slow_rhs_inc (int /*level*/, int nrk, @@ -101,9 +97,7 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v, const amrex::Real* dptr_rhotheta_src, - const amrex::Real* dptr_rayleigh_tau, const amrex::Real* dptr_rayleigh_ubar, - const amrex::Real* dptr_rayleigh_vbar, const amrex::Real* dptr_rayleigh_wbar, - const amrex::Real* dptr_rayleigh_thetabar) + const Vector d_rayleigh_dptrs) { BL_PROFILE_REGION("erf_slow_rhs_pre()"); @@ -141,6 +135,12 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, const GpuArray dxInv = geom.InvCellSizeArray(); + Real* tau = d_rayleigh_ptrs_at_lev[Rayleigh::tau]; + Real* ubar = d_rayleigh_ptrs_at_lev[Rayleigh::ubar]; + Real* vbar = d_rayleigh_ptrs_at_lev[Rayleigh::vbar]; + Real* wbar = d_rayleigh_ptrs_at_lev[Rayleigh::wbar]; + Real* thetabar = d_rayleigh_ptrs_at_lev[Rayleigh::thetabar]; + // ************************************************************************* // Combine external forcing terms // ************************************************************************* @@ -668,7 +668,7 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real theta = cell_prim(i,j,k,np); - cell_rhs(i, j, k, n) -= dptr_rayleigh_tau[k] * (theta - dptr_rayleigh_thetabar[k]) * cell_data(i,j,k,nr); + cell_rhs(i, j, k, n) -= tau[k] * (theta - thetabar[k]) * cell_data(i,j,k,nr); }); } @@ -721,9 +721,11 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // x-momentum equation + Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) ); + // Note we do NOT include a pressure gradient here rho_u_rhs(i, j, k) += - solverChoice.abl_pressure_grad[0] - + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i-1,j,k,Rho_comp)) * solverChoice.abl_geo_forcing[0]; + + rho_on_u_face * solverChoice.abl_geo_forcing[0]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (solverChoice.use_coriolis) @@ -735,10 +737,10 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, } // Add Rayleigh damping - if (solverChoice.use_rayleigh_damping && solverChoice.rayleigh_damp_U) + if (use_rayleigh_damping && solverChoice.rayleigh_damp_U) { - Real uu = rho_u(i,j,k) / cell_data(i,j,k,Rho_comp); - rho_u_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (uu - dptr_rayleigh_ubar[k]) * cell_data(i,j,k,Rho_comp); + Real uu = rho_u(i,j,k) / rho_on_u_face; + rho_u_rhs(i, j, k) -= tau[k] * (uu - ubar[k]) * rho_on_u_face; } if (nrk == 1) { @@ -757,9 +759,11 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // y-momentum equation + Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) ); + // Note we do NOT include a pressure gradient here rho_v_rhs(i, j, k) += - solverChoice.abl_pressure_grad[1] - + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i,j-1,k,Rho_comp)) * solverChoice.abl_geo_forcing[1]; + + rho_v_face * solverChoice.abl_geo_forcing[1]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (solverChoice.use_coriolis) @@ -769,12 +773,13 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, } // Add Rayleigh damping - if (solverChoice.use_rayleigh_damping && solverChoice.rayleigh_damp_V) + if (use_rayleigh_damping && solverChoice.rayleigh_damp_V) { - Real vv = rho_v(i,j,k) / cell_data(i,j,k,Rho_comp); - rho_v_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (vv - dptr_rayleigh_vbar[k]) * cell_data(i,j,k,Rho_comp); + Real vv = rho_v(i,j,k) / rho_on_v_face; + rho_v_rhs(i, j, k) -= tau[k] * (vv - vbar[k]) * rho_on_v_face; } + if (nrk == 1) { rho_v_rhs(i,j,k) *= 0.5; rho_v_rhs(i,j,k) += 0.5 / dt * (rho_v(i,j,k) - rho_v_old(i,j,k)); @@ -803,9 +808,11 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // z-momentum equation + Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) ); + // Note we do NOT include a pressure gradient or buoyancy here rho_w_rhs(i, j, k) += - solverChoice.abl_pressure_grad[2] - + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i,j,k-1,Rho_comp)) * solverChoice.abl_geo_forcing[2]; + rho_on_w_face * solverChoice.abl_geo_forcing[2]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (solverChoice.use_coriolis) @@ -815,12 +822,13 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, } // Add Rayleigh damping - if (solverChoice.use_rayleigh_damping && solverChoice.rayleigh_damp_W) + if (use_rayleigh_damping && solverChoice.rayleigh_damp_W) { - Real ww = rho_w(i,j,k) / cell_data(i,j,k,Rho_comp); - rho_w_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (ww - dptr_rayleigh_wbar[k]) * cell_data(i,j,k,Rho_comp); + Real ww = rho_w(i,j,k) / rho_on_w_face; + rho_w_rhs(i, j, k) -= tau[k] * (ww - wbar[k]) * rho_on_w_face; } + if (nrk == 1) { rho_w_rhs(i,j,k) *= 0.5; rho_w_rhs(i,j,k) += 0.5 / dt * (rho_w(i,j,k) - rho_w_old(i,j,k)); diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 0a908f77d..010213f66 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -64,11 +64,7 @@ using namespace amrex; * @param[inout] fr_as_crse YAFluxRegister at level l at level l / l+1 interface * @param[inout] fr_as_fine YAFluxRegister at level l at level l-1 / l interface * @param[in] dptr_rhotheta_src custom temperature source term - * @param[in] dptr_rayleigh_tau strength of Rayleigh damping - * @param[in] dptr_rayleigh_ubar reference value for x-velocity used to define Rayleigh damping - * @param[in] dptr_rayleigh_vbar reference value for y-velocity used to define Rayleigh damping - * @param[in] dptr_rayleigh_wbar reference value for z-velocity used to define Rayleigh damping - * @param[in] dptr_rayleigh_thetabar reference value for potential temperature used to define Rayleigh damping + * @param[in] d_rayleigh_ptrs_at_lev Vector of {strength of Rayleigh damping, reference value for xvel/yvel/zvel/theta} used to define Rayleigh damping */ void erf_slow_rhs_pre (int level, int finest_level, @@ -107,9 +103,7 @@ void erf_slow_rhs_pre (int level, int finest_level, YAFluxRegister* fr_as_crse, YAFluxRegister* fr_as_fine, const amrex::Real* dptr_rhotheta_src, - const amrex::Real* dptr_rayleigh_tau, const amrex::Real* dptr_rayleigh_ubar, - const amrex::Real* dptr_rayleigh_vbar, const amrex::Real* dptr_rayleigh_wbar, - const amrex::Real* dptr_rayleigh_thetabar) + const Vector d_rayleigh_ptrs_at_lev) { BL_PROFILE_REGION("erf_slow_rhs_pre()"); @@ -154,6 +148,12 @@ void erf_slow_rhs_pre (int level, int finest_level, const GpuArray dxInv = geom.InvCellSizeArray(); const Real* dx = geom.CellSize(); + Real* tau = d_rayleigh_ptrs_at_lev[Rayleigh::tau]; + Real* ubar = d_rayleigh_ptrs_at_lev[Rayleigh::ubar]; + Real* vbar = d_rayleigh_ptrs_at_lev[Rayleigh::vbar]; + Real* wbar = d_rayleigh_ptrs_at_lev[Rayleigh::wbar]; + Real* thetabar = d_rayleigh_ptrs_at_lev[Rayleigh::thetabar]; + // ************************************************************************* // Combine external forcing terms // ************************************************************************* @@ -748,7 +748,7 @@ void erf_slow_rhs_pre (int level, int finest_level, ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real theta = cell_prim(i,j,k,np); - cell_rhs(i, j, k, n) -= dptr_rayleigh_tau[k] * (theta - dptr_rayleigh_thetabar[k]) * cell_data(i,j,k,nr); + cell_rhs(i, j, k, n) -= tau[k] * (theta - thetabar[k]) * cell_data(i,j,k,nr); }); } @@ -783,14 +783,14 @@ void erf_slow_rhs_pre (int level, int finest_level, tau12, tau13, tau21, tau23, tau31, tau32, - cell_data, detJ_arr, dc, dxInv, + detJ_arr, dxInv, mf_m, mf_u, mf_v); } else { DiffusionSrcForMom_N(tbx, tby, tbz, rho_u_rhs, rho_v_rhs, rho_w_rhs, tau11, tau22, tau33, tau12, tau13, tau23, - cell_data, dc, dxInv, + dxInv, mf_m, mf_u, mf_v); } } @@ -823,7 +823,7 @@ void erf_slow_rhs_pre (int level, int finest_level, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // x-momentum equation - Real rho_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) ); + Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) ); Real met_h_xi = Compute_h_xi_AtIface (i, j, k, dxInv, z_nd); Real met_h_zeta = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd); @@ -853,7 +853,7 @@ void erf_slow_rhs_pre (int level, int finest_level, } rho_u_rhs(i, j, k) += (-gpx - abl_pressure_grad[0]) / (1.0 + q) - + rho_u_face * abl_geo_forcing[0]; + + rho_on_u_face * abl_geo_forcing[0]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) @@ -866,8 +866,8 @@ void erf_slow_rhs_pre (int level, int finest_level, // Add Rayleigh damping if (use_rayleigh_damping && rayleigh_damp_U) { - Real uu = rho_u(i,j,k) / rho_u_face; - rho_u_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (uu - dptr_rayleigh_ubar[k]) * rho_u_face; + Real uu = rho_u(i,j,k) / rho_on_u_face; + rho_u_rhs(i, j, k) -= tau[k] * (uu - ubar[k]) * rho_on_u_face; } if (l_moving_terrain) { @@ -884,7 +884,7 @@ void erf_slow_rhs_pre (int level, int finest_level, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // x-momentum equation - Real rho_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) ); + Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) ); Real gpx = dxInv[0] * (pp_arr(i,j,k) - pp_arr(i-1,j,k)); gpx *= mf_u(i,j,0); @@ -895,7 +895,7 @@ void erf_slow_rhs_pre (int level, int finest_level, } rho_u_rhs(i, j, k) += (-gpx - abl_pressure_grad[0]) / (1.0 + q) - + rho_u_face * abl_geo_forcing[0]; + + rho_on_u_face * abl_geo_forcing[0]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) @@ -909,7 +909,7 @@ void erf_slow_rhs_pre (int level, int finest_level, if (use_rayleigh_damping && rayleigh_damp_U) { Real uu = rho_u(i,j,k) / cell_data(i,j,k,Rho_comp); - rho_u_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (uu - dptr_rayleigh_ubar[k]) * cell_data(i,j,k,Rho_comp); + rho_u_rhs(i, j, k) -= tau[k] * (uu - ubar[k]) * rho_on_u_face; } }); } // no terrain @@ -926,7 +926,7 @@ void erf_slow_rhs_pre (int level, int finest_level, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // y-momentum equation - Real rho_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) ); + Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) ); Real met_h_eta = Compute_h_eta_AtJface (i, j, k, dxInv, z_nd); Real met_h_zeta = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd); @@ -956,7 +956,7 @@ void erf_slow_rhs_pre (int level, int finest_level, +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j-1,k,PrimQ2_comp) ); } rho_v_rhs(i, j, k) += (-gpy - abl_pressure_grad[1]) / (1.0_rt + q) - + rho_v_face * abl_geo_forcing[1]; + + rho_on_v_face * abl_geo_forcing[1]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) { @@ -967,8 +967,8 @@ void erf_slow_rhs_pre (int level, int finest_level, // Add Rayleigh damping if (use_rayleigh_damping && rayleigh_damp_V) { - Real vv = rho_v(i,j,k) / rho_v_face; - rho_v_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (vv - dptr_rayleigh_vbar[k]) * rho_v_face; + Real vv = rho_v(i,j,k) / rho_on_v_face; + rho_v_rhs(i, j, k) -= tau[k] * (vv - vbar[k]) * rho_on_v_face; } if (l_moving_terrain) { @@ -985,7 +985,7 @@ void erf_slow_rhs_pre (int level, int finest_level, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // y-momentum equation - Real rho_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) ); + Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) ); Real gpy = dxInv[1] * (pp_arr(i,j,k) - pp_arr(i,j-1,k)); gpy *= mf_v(i,j,0); @@ -996,7 +996,7 @@ void erf_slow_rhs_pre (int level, int finest_level, } rho_v_rhs(i, j, k) += (-gpy - abl_pressure_grad[1]) / (1.0_rt + q) - + rho_v_face * abl_geo_forcing[1]; + + rho_on_v_face * abl_geo_forcing[1]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) @@ -1008,8 +1008,8 @@ void erf_slow_rhs_pre (int level, int finest_level, // Add Rayleigh damping if (use_rayleigh_damping && rayleigh_damp_V) { - Real vv = rho_v(i,j,k) / rho_v_face; - rho_v_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (vv - dptr_rayleigh_vbar[k]) * rho_v_face; + Real vv = rho_v(i,j,k) / rho_on_v_face; + rho_v_rhs(i, j, k) -= tau[k] * (vv - vbar[k]) * rho_on_v_face; } }); } // no terrain @@ -1037,7 +1037,7 @@ void erf_slow_rhs_pre (int level, int finest_level, if (l_use_terrain) { ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // z-momentum equation - Real rho_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) ); + Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) ); Real met_h_zeta = Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd); Real gpz = dxInv[2] * ( pp_arr(i,j,k)-pp_arr(i,j,k-1) ) / met_h_zeta; @@ -1047,7 +1047,7 @@ void erf_slow_rhs_pre (int level, int finest_level, +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j,k-1,PrimQ2_comp) ); } rho_w_rhs(i, j, k) += (buoyancy_fab(i,j,k) - gpz - abl_pressure_grad[2]) / (1.0_rt + q) - + rho_w_face * abl_geo_forcing[2]; + + rho_on_w_face * abl_geo_forcing[2]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) @@ -1059,8 +1059,8 @@ void erf_slow_rhs_pre (int level, int finest_level, // Add Rayleigh damping if (use_rayleigh_damping && rayleigh_damp_W) { - Real ww = rho_w(i,j,k) / rho_w_face; - rho_w_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (ww - dptr_rayleigh_wbar[k]) * rho_w_face; + Real ww = rho_w(i,j,k) / rho_on_w_face; + rho_w_rhs(i, j, k) -= tau[k] * (ww - wbar[k]) * rho_on_w_face; } if (l_use_terrain && l_moving_terrain) { @@ -1075,7 +1075,7 @@ void erf_slow_rhs_pre (int level, int finest_level, ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // z-momentum equation - Real rho_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) ); + Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) ); Real gpz = dxInv[2] * ( pp_arr(i,j,k)-pp_arr(i,j,k-1) ); Real q = 0.0; @@ -1084,7 +1084,7 @@ void erf_slow_rhs_pre (int level, int finest_level, +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j,k-1,PrimQ2_comp) ); } rho_w_rhs(i, j, k) += (buoyancy_fab(i,j,k) - gpz - abl_pressure_grad[2]) / (1.0_rt + q) - + rho_w_face * abl_geo_forcing[2]; + + rho_on_w_face * abl_geo_forcing[2]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) @@ -1096,8 +1096,8 @@ void erf_slow_rhs_pre (int level, int finest_level, // Add Rayleigh damping if (use_rayleigh_damping && rayleigh_damp_W) { - Real ww = rho_w(i,j,k) / rho_w_face; - rho_w_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (ww - dptr_rayleigh_wbar[k]) * rho_w_face; + Real ww = rho_w(i,j,k) / rho_on_w_face; + rho_w_rhs(i, j, k) -= tau[k] * (ww - wbar[k]) * rho_on_w_face; } }); } // no terrain diff --git a/Source/TimeIntegration/TI_headers.H b/Source/TimeIntegration/TI_headers.H index dd67dc728..1cd677480 100644 --- a/Source/TimeIntegration/TI_headers.H +++ b/Source/TimeIntegration/TI_headers.H @@ -56,11 +56,7 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk, amrex::YAFluxRegister* fr_as_crse, amrex::YAFluxRegister* fr_as_fine, const amrex::Real* dptr_rhotheta_src, - const amrex::Real* dptr_rayleigh_tau, - const amrex::Real* dptr_rayleigh_ubar, - const amrex::Real* dptr_rayleigh_vbar, - const amrex::Real* dptr_rayleigh_wbar, - const amrex::Real* dptr_rayleigh_thetabar); + const amrex::Vector d_rayleigh_ptrs_at_lev); /** * Function for computing the slow RHS for the evolution equations for the scalars other than density or potential temperature @@ -271,11 +267,7 @@ void erf_slow_rhs_inc (int level, int nrk, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v, const amrex::Real* dptr_rhotheta_src, - const amrex::Real* dptr_rayleigh_tau, - const amrex::Real* dptr_rayleigh_ubar, - const amrex::Real* dptr_rayleigh_vbar, - const amrex::Real* dptr_rayleigh_wbar, - const amrex::Real* dptr_rayleigh_thetabar); + const Vector d_rayleigh_ptrs_at_lev); #endif void ApplySpongeZoneBCs (const SpongeChoice& spongeChoice, diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 78d17c53e..5edcc03c4 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -106,10 +106,7 @@ z_phys_nd_src[level], detJ_cc_src[level], p0_new, mapfac_m[level], mapfac_u[level], mapfac_v[level], fr_as_crse, fr_as_fine, - dptr_rhotheta_src, - dptr_rayleigh_tau, dptr_rayleigh_ubar, - dptr_rayleigh_vbar, dptr_rayleigh_wbar, - dptr_rayleigh_thetabar); + dptr_rhotheta_src, d_rayleigh_ptrs_at_lev); // We define and evolve (rho theta)_0 in order to re-create p_0 in a way that is consistent // with our update of (rho theta) but does NOT maintain dp_0 / dz = -rho_0 g. This is why @@ -202,10 +199,7 @@ z_phys_nd[level], detJ_cc[level], p0, mapfac_m[level], mapfac_u[level], mapfac_v[level], fr_as_crse, fr_as_fine, - dptr_rhotheta_src, - dptr_rayleigh_tau, dptr_rayleigh_ubar, - dptr_rayleigh_vbar, dptr_rayleigh_wbar, - dptr_rayleigh_thetabar); + dptr_rhotheta_src, d_rayleigh_ptrs_at_lev); } #ifdef ERF_USE_NETCDF @@ -386,10 +380,7 @@ fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, z_phys_nd[level], detJ_cc[level], p0, mapfac_m[level], mapfac_u[level], mapfac_v[level], - dptr_rhotheta_src, - dptr_rayleigh_tau, dptr_rayleigh_ubar, - dptr_rayleigh_vbar, dptr_rayleigh_wbar, - dptr_rayleigh_thetabar); + dptr_rhotheta_src, d_rayleigh_ptrs_at_lev); #ifdef ERF_USE_NETCDF diff --git a/Source/prob_common.H b/Source/prob_common.H index 4c9db44cc..cc6342385 100644 --- a/Source/prob_common.H +++ b/Source/prob_common.H @@ -127,11 +127,11 @@ public: amrex::Error("Temperature forcing not defined for "+name()+" problem with terrain"); } else { const int khi = geom.Domain().bigEnd()[2]; - const amrex::Real* prob_lo = geom.ProbLo(); - const auto dx = geom.CellSize(); + // const amrex::Real* prob_lo = geom.ProbLo(); + // const auto dx = geom.CellSize(); for (int k = 0; k <= khi; k++) { - const amrex::Real z_cc = prob_lo[2] + (k+0.5)* dx[2]; + // const amrex::Real z_cc = prob_lo[2] + (k+0.5)* dx[2]; // set RHS term of RhoTheta equation based on time, z_cc here src[k] = 0.0; } @@ -187,20 +187,12 @@ public: /** * Function to define the quantities needed to impose Rayleigh damping * - * @param[out] tau strength of Rayleigh damping - * @param[out] ubar reference value for x-velocity used to define Rayleigh damping - * @param[out] vbar reference value for y-velocity used to define Rayleigh damping - * @param[out] wbar reference value for z-velocity used to define Rayleigh damping - * @param[out] thetabar reference value for potential temperature used to define Rayleigh damping + * @param[out] rayleigh_ptrs = {strength of Rayleigh damping, reference values for xvel/yvel/zvel/theta used to define Rayleigh damping} * @param[in] geom container for geometric information * @param[in] z_phys_cc height coordinate at cell centers */ virtual void - erf_init_rayleigh (amrex::Vector& /*tau*/, - amrex::Vector& /*ubar*/, - amrex::Vector& /*vbar*/, - amrex::Vector& /*wbar*/, - amrex::Vector& /*thetabar*/, + erf_init_rayleigh (amrex::Vector >& /*rayleigh_ptrs*/, amrex::Geometry const& /*geom*/, std::unique_ptr& /*z_phys_cc*/) {