diff --git a/Exec/DevTests/Bomex/input_SAM b/Exec/DevTests/Bomex/input_SAM index f2f1e306d..bf904fce2 100644 --- a/Exec/DevTests/Bomex/input_SAM +++ b/Exec/DevTests/Bomex/input_SAM @@ -56,6 +56,7 @@ erf.init_sounding_ideal = true erf.add_custom_rhotheta_forcing = true erf.add_custom_moisture_forcing = true +erf.add_custom_geostrophic_profile = true erf.add_custom_w_subsidence = true erf.custom_forcing_uses_primitive_vars = false diff --git a/Exec/DevTests/Bomex/prob.H b/Exec/DevTests/Bomex/prob.H index fb5d3bc76..e80a72a5f 100644 --- a/Exec/DevTests/Bomex/prob.H +++ b/Exec/DevTests/Bomex/prob.H @@ -107,6 +107,15 @@ public: const amrex::Geometry& geom, std::unique_ptr& z_phys_cc) override; + void update_geostrophic_profile ( + const amrex::Real& /*time*/, + amrex::Vector& u_geos, + amrex::Gpu::DeviceVector& d_u_geos, + amrex::Vector& v_geos, + amrex::Gpu::DeviceVector& d_v_geos, + const amrex::Geometry& geom, + std::unique_ptr& z_phys_cc) override; + protected: std::string name () override { return "BOMEX"; } diff --git a/Exec/DevTests/Bomex/prob.cpp b/Exec/DevTests/Bomex/prob.cpp index a0e9c8085..d4a556bcb 100644 --- a/Exec/DevTests/Bomex/prob.cpp +++ b/Exec/DevTests/Bomex/prob.cpp @@ -10,7 +10,7 @@ amrex_probinit (const amrex_real* problo, const amrex_real* probhi) return std::make_unique(problo, probhi); } -Problem::Problem (const amrex::Real* problo, const amrex::Real* probhi) +Problem::Problem (const Real* problo, const Real* probhi) { // Parse params ParmParse pp("prob"); @@ -60,28 +60,28 @@ Problem::Problem (const amrex::Real* problo, const amrex::Real* probhi) void Problem::init_custom_pert ( - const amrex::Box& bx, - const amrex::Box& xbx, - const amrex::Box& ybx, - const amrex::Box& zbx, - amrex::Array4 const& /*state*/, - amrex::Array4 const& state_pert, - amrex::Array4 const& x_vel_pert, - amrex::Array4 const& y_vel_pert, - amrex::Array4 const& z_vel_pert, - amrex::Array4 const& /*r_hse*/, - amrex::Array4 const& /*p_hse*/, - amrex::Array4 const& /*z_nd*/, - amrex::Array4 const& /*z_cc*/, - amrex::GeometryData const& geomdata, - amrex::Array4 const& /*mf_m*/, - amrex::Array4 const& /*mf_u*/, - amrex::Array4 const& /*mf_v*/, + const Box& bx, + const Box& xbx, + const Box& ybx, + const Box& zbx, + Array4 const& /*state*/, + Array4 const& state_pert, + Array4 const& x_vel_pert, + Array4 const& y_vel_pert, + Array4 const& z_vel_pert, + Array4 const& /*r_hse*/, + Array4 const& /*p_hse*/, + Array4 const& /*z_nd*/, + Array4 const& /*z_cc*/, + GeometryData const& geomdata, + Array4 const& /*mf_m*/, + Array4 const& /*mf_u*/, + Array4 const& /*mf_v*/, const SolverChoice& sc) { const bool use_moisture = (sc.moisture_type != MoistureType::None); - ParallelForRNG(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept + ParallelForRNG(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const RandomEngine& engine) noexcept { // Geometry const Real* prob_lo = geomdata.ProbLo(); @@ -117,7 +117,7 @@ Problem::init_custom_pert ( }); // Set the x-velocity - ParallelForRNG(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept + ParallelForRNG(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const RandomEngine& engine) noexcept { const Real* prob_lo = geomdata.ProbLo(); const Real* dx = geomdata.CellSize(); @@ -142,7 +142,7 @@ Problem::init_custom_pert ( }); // Set the y-velocity - ParallelForRNG(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept + ParallelForRNG(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const RandomEngine& engine) noexcept { const Real* prob_lo = geomdata.ProbLo(); const Real* dx = geomdata.CellSize(); @@ -167,7 +167,7 @@ Problem::init_custom_pert ( }); // Set the z-velocity - ParallelForRNG(zbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept + ParallelForRNG(zbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const RandomEngine& engine) noexcept { const int dom_lo_z = geomdata.Domain().smallEnd()[2]; const int dom_hi_z = geomdata.Domain().bigEnd()[2]; @@ -190,24 +190,24 @@ Problem::init_custom_pert ( // USER-DEFINED FUNCTION //============================================================================= void -Problem::update_rhotheta_sources (const amrex::Real& /*time*/, - amrex::Vector& src, - amrex::Gpu::DeviceVector& d_src, - const amrex::Geometry& geom, - std::unique_ptr& z_phys_cc) +Problem::update_rhotheta_sources (const Real& /*time*/, + Vector& src, + Gpu::DeviceVector& d_src, + const Geometry& geom, + std::unique_ptr& z_phys_cc) { if (src.empty()) return; - const int khi = geom.Domain().bigEnd()[2]; - const amrex::Real* prob_lo = geom.ProbLo(); - const auto dx = geom.CellSize(); + const int khi = geom.Domain().bigEnd()[2]; + const Real* prob_lo = geom.ProbLo(); + const auto dx = geom.CellSize(); // Note: If z_phys_cc, then use_terrain=1 was set. If the z coordinate // varies in time and or space, then the the height needs to be // calculated at each time step. Here, we assume that only grid // stretching exists. if (z_phys_cc && zlevels.empty()) { - amrex::Print() << "Initializing z levels on stretched grid" << std::endl; + Print() << "Initializing z levels on stretched grid" << std::endl; zlevels.resize(khi+1); reduce_to_max_per_level(zlevels, z_phys_cc); } @@ -232,24 +232,24 @@ Problem::update_rhotheta_sources (const amrex::Real& /*time*/, // USER-DEFINED FUNCTION //============================================================================= void -Problem::update_rhoqt_sources (const amrex::Real& /*time*/, - amrex::Vector& qsrc, - amrex::Gpu::DeviceVector& d_qsrc, - const amrex::Geometry& geom, - std::unique_ptr& z_phys_cc) +Problem::update_rhoqt_sources (const Real& /*time*/, + Vector& qsrc, + Gpu::DeviceVector& d_qsrc, + const Geometry& geom, + std::unique_ptr& z_phys_cc) { if (qsrc.empty()) return; - const int khi = geom.Domain().bigEnd()[2]; - const amrex::Real* prob_lo = geom.ProbLo(); - const auto dx = geom.CellSize(); + const int khi = geom.Domain().bigEnd()[2]; + const Real* prob_lo = geom.ProbLo(); + const auto dx = geom.CellSize(); // Note: If z_phys_cc, then use_terrain=1 was set. If the z coordinate // varies in time and or space, then the the height needs to be // calculated at each time step. Here, we assume that only grid // stretching exists. if (z_phys_cc && zlevels.empty()) { - amrex::Print() << "Initializing z levels on stretched grid" << std::endl; + Print() << "Initializing z levels on stretched grid" << std::endl; zlevels.resize(khi+1); reduce_to_max_per_level(zlevels, z_phys_cc); } @@ -274,24 +274,24 @@ Problem::update_rhoqt_sources (const amrex::Real& /*time*/, // USER-DEFINED FUNCTION //============================================================================= void -Problem::update_w_subsidence (const amrex::Real& /*time*/, - amrex::Vector& wbar, - amrex::Gpu::DeviceVector& d_wbar, - const amrex::Geometry& geom, - std::unique_ptr& z_phys_cc) +Problem::update_w_subsidence (const Real& /*time*/, + Vector& wbar, + Gpu::DeviceVector& d_wbar, + const Geometry& geom, + std::unique_ptr& z_phys_cc) { if (wbar.empty()) return; - const int khi = geom.Domain().bigEnd()[2]; - const amrex::Real* prob_lo = geom.ProbLo(); - const auto dx = geom.CellSize(); + const int khi = geom.Domain().bigEnd()[2]; + const Real* prob_lo = geom.ProbLo(); + const auto dx = geom.CellSize(); // Note: If z_phys_cc, then use_terrain=1 was set. If the z coordinate // varies in time and or space, then the the height needs to be // calculated at each time step. Here, we assume that only grid // stretching exists. if (z_phys_cc && zlevels.empty()) { - amrex::Print() << "Initializing z levels on stretched grid" << std::endl; + Print() << "Initializing z levels on stretched grid" << std::endl; zlevels.resize(khi+1); reduce_to_max_per_level(zlevels, z_phys_cc); } @@ -315,3 +315,46 @@ Problem::update_w_subsidence (const amrex::Real& /*time*/, // Copy from host version to device version amrex::Gpu::copy(amrex::Gpu::hostToDevice, wbar.begin(), wbar.end(), d_wbar.begin()); } + +//============================================================================= +// USER-DEFINED FUNCTION +//============================================================================= +void +Problem::update_geostrophic_profile (const Real& /*time*/, + Vector& u_geos, + Gpu::DeviceVector& d_u_geos, + Vector& v_geos, + Gpu::DeviceVector& d_v_geos, + const Geometry& geom, + std::unique_ptr& z_phys_cc) +{ + if (u_geos.empty()) return; + + const int khi = geom.Domain().bigEnd()[2]; + const Real* prob_lo = geom.ProbLo(); + const auto dx = geom.CellSize(); + + // Note: If z_phys_cc, then use_terrain=1 was set. If the z coordinate + // varies in time and or space, then the the height needs to be + // calculated at each time step. Here, we assume that only grid + // stretching exists. + if (z_phys_cc && zlevels.empty()) { + Print() << "Initializing z levels on stretched grid" << std::endl; + zlevels.resize(khi+1); + reduce_to_max_per_level(zlevels, z_phys_cc); + } + + const Real coriolis = 0.376E-4; + + // Only apply temperature source below nominal inversion height + for (int k = 0; k <= khi; k++) { + const Real z_cc = (z_phys_cc) ? zlevels[k] : prob_lo[2] + (k+0.5)* dx[2]; + const Real u_geo_wind = -10.0 + z_cc * 0.0018; + u_geos[k] = 0; // -coriolis_factor * v_geo_wind + v_geos[k] = coriolis * u_geo_wind; + } + + // Copy from host version to device version + amrex::Gpu::copy(amrex::Gpu::hostToDevice, u_geos.begin(), u_geos.end(), d_u_geos.begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, v_geos.begin(), v_geos.end(), d_v_geos.begin()); +} diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 96416ca99..8bdfd8f11 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -216,6 +216,7 @@ struct SolverChoice { pp.query("add_custom_rhotheta_forcing", custom_rhotheta_forcing); pp.query("add_custom_moisture_forcing", custom_moisture_forcing); pp.query("add_custom_w_subsidence", custom_w_subsidence); + pp.query("add_custom_geostrophic_profile", custom_geostrophic_profile); pp.query("custom_forcing_uses_primitive_vars", custom_forcing_prim_vars); pp.query("Ave_Plane", ave_plane); @@ -403,6 +404,7 @@ struct SolverChoice { bool custom_rhotheta_forcing = false; bool custom_moisture_forcing = false; bool custom_w_subsidence = false; + bool custom_geostrophic_profile = false; bool custom_forcing_prim_vars = false; // Numerical diffusion diff --git a/Source/ERF.H b/Source/ERF.H index 81044419e..b35c694fd 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -818,6 +818,12 @@ private: amrex::Vector< amrex::Vector > h_w_subsid; amrex::Vector > d_w_subsid; + amrex::Vector< amrex::Vector > h_u_geos; + amrex::Vector > d_u_geos; + + amrex::Vector< amrex::Vector > h_v_geos; + amrex::Vector > d_v_geos; + // This is a vector over levels of vectors across quantities of Vectors amrex::Vector > > h_rayleigh_ptrs; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index fac3b0213..64c4f0184 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -677,6 +677,25 @@ ERF::InitData () } } + if (solverChoice.custom_geostrophic_profile) + { + h_u_geos.resize(max_level+1, Vector(0)); + d_u_geos.resize(max_level+1, Gpu::DeviceVector(0)); + h_v_geos.resize(max_level+1, Vector(0)); + d_v_geos.resize(max_level+1, Gpu::DeviceVector(0)); + for (int lev = 0; lev <= finest_level; lev++) { + const int domlen = geom[lev].Domain().length(2); + h_u_geos[lev].resize(domlen, 0.0_rt); + d_u_geos[lev].resize(domlen, 0.0_rt); + h_v_geos[lev].resize(domlen, 0.0_rt); + d_v_geos[lev].resize(domlen, 0.0_rt); + prob->update_geostrophic_profile(t_new[0], + h_u_geos[lev], d_u_geos[lev], + h_v_geos[lev], d_v_geos[lev], + geom[lev], z_phys_cc[lev]); + } + } + if (solverChoice.custom_moisture_forcing) { h_rhoqt_src.resize(max_level+1, Vector(0)); diff --git a/Source/Microphysics/SAM/IceFall.cpp b/Source/Microphysics/SAM/IceFall.cpp index cbf723246..120b7c725 100644 --- a/Source/Microphysics/SAM/IceFall.cpp +++ b/Source/Microphysics/SAM/IceFall.cpp @@ -36,15 +36,10 @@ void SAM::IceFall () { auto rho_array = rho->array(mfi); auto fz_array = fz.array(mfi); - const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4{}; - const auto& box3d = mfi.tilebox(); ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - // Jacobian determinant - Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0; - Real rho_avg, qci_avg; if (k==k_lo) { rho_avg = rho_array(i,j,k); diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 0a6dc286f..70dc2f802 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -59,6 +59,8 @@ void ERF::advance_dycore(int level, Real* dptr_rhotheta_src = solverChoice.custom_rhotheta_forcing ? d_rhotheta_src[level].data() : nullptr; Real* dptr_rhoqt_src = solverChoice.custom_moisture_forcing ? d_rhoqt_src[level].data() : nullptr; Real* dptr_wbar_sub = solverChoice.custom_w_subsidence ? d_w_subsid[level].data() : nullptr; + Real* dptr_u_geos = solverChoice.custom_geostrophic_profile ? d_u_geos[level].data() : nullptr; + Real* dptr_v_geos = solverChoice.custom_geostrophic_profile ? d_v_geos[level].data() : nullptr; Vector d_rayleigh_ptrs_at_lev; d_rayleigh_ptrs_at_lev.resize(Rayleigh::nvars); @@ -210,6 +212,13 @@ void ERF::advance_dycore(int level, fine_geom, z_phys_cc[level]); } + if (solverChoice.custom_geostrophic_profile) { + prob->update_geostrophic_profile(old_time, + h_u_geos[level], d_u_geos[level], + h_v_geos[level], d_v_geos[level], + fine_geom, z_phys_cc[level]); + } + // *********************************************************************************************** // Convert old velocity available on faces to old momentum on faces to be used in time integration // *********************************************************************************************** diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index 7557f58be..74ffb54f7 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -63,6 +63,8 @@ using namespace amrex; * @param[in] mapfac_v map factor at y-faces * @param[in] dptr_rhotheta_src custom temperature source term * @param[in] dptr_rhoqt_src custom moisture source term + * @param[in] dptr_u_geos custom geostrophic wind profile + * @param[in] dptr_v_geos custom geostrophic wind profile * @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 */ @@ -97,6 +99,8 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v, const Real* dptr_rhotheta_src, + const Real* dptr_u_geos, + const Real* dptr_v_geos, const Real* dptr_wbar_sub, const Vector d_rayleigh_dptrs) { @@ -752,6 +756,11 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, rho_u_rhs(i, j, k) += - solverChoice.abl_pressure_grad[0] + rho_on_u_face * solverChoice.abl_geo_forcing[0]; + if (solverChoice.custom_geostrophic_profile) { + rho_u_rhs(i, j, k) += - solverChoice.abl_pressure_grad[0] + + rho_on_u_face * dptr_u_geos[k]; + } + // Add Coriolis forcing (that assumes east is +x, north is +y) if (solverChoice.use_coriolis) { @@ -790,6 +799,11 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, rho_v_rhs(i, j, k) += - solverChoice.abl_pressure_grad[1] + rho_v_face * solverChoice.abl_geo_forcing[1]; + if (solverChoice.custom_geostrophic_profile) { + rho_v_rhs(i, j, k) += - solverChoice.abl_pressure_grad[1] + + rho_v_face * dptr_v_geos[k]; + } + // Add Coriolis forcing (that assumes east is +x, north is +y) if (solverChoice.use_coriolis) { diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 1df46dfba..72090c3d6 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -64,6 +64,8 @@ using namespace amrex; * @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_rhoqt_src custom moisture source term + * @param[in] dptr_u_geos custom geostrophic wind profile + * @param[in] dptr_v_geos custom geostrophic wind profile * @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 */ @@ -103,6 +105,8 @@ void erf_slow_rhs_pre (int level, int finest_level, YAFluxRegister* fr_as_crse, YAFluxRegister* fr_as_fine, const Real* dptr_rhotheta_src, + const Real* dptr_u_geos, + const Real* dptr_v_geos, const Real* dptr_wbar_sub, const Vector d_rayleigh_ptrs_at_lev) { @@ -929,7 +933,8 @@ void erf_slow_rhs_pre (int level, int finest_level, { BL_PROFILE("slow_rhs_pre_xmom"); - auto rayleigh_damp_U = solverChoice.rayleigh_damp_U; + auto rayleigh_damp_U = solverChoice.rayleigh_damp_U; + auto geostrophic_wind = solverChoice.custom_geostrophic_profile; // ***************************************************************************** // TERRAIN VERSION // ***************************************************************************** @@ -970,17 +975,20 @@ 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_on_u_face * abl_geo_forcing[0]; + // Add geostrophic wind + if (geostrophic_wind) { + rho_u_rhs(i, j, k) += rho_on_u_face * dptr_u_geos[k]; + } + // Add Coriolis forcing (that assumes east is +x, north is +y) - if (use_coriolis) - { + if (use_coriolis) { Real rho_v_loc = 0.25 * (rho_v(i,j+1,k) + rho_v(i,j,k) + rho_v(i-1,j+1,k) + rho_v(i-1,j,k)); Real rho_w_loc = 0.25 * (rho_w(i,j,k+1) + rho_w(i,j,k) + rho_w(i,j-1,k+1) + rho_w(i,j-1,k)); rho_u_rhs(i, j, k) += coriolis_factor * (rho_v_loc * sinphi - rho_w_loc * cosphi); } // Add Rayleigh damping - if (use_rayleigh_damping && rayleigh_damp_U) - { + if (use_rayleigh_damping && rayleigh_damp_U) { 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; } @@ -1012,17 +1020,20 @@ 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_on_u_face * abl_geo_forcing[0]; + // Add geostrophic wind + if (geostrophic_wind) { + rho_u_rhs(i, j, k) += rho_on_u_face * dptr_u_geos[k]; + } + // Add Coriolis forcing (that assumes east is +x, north is +y) - if (use_coriolis) - { + if (use_coriolis) { Real rho_v_loc = 0.25 * (rho_v(i,j+1,k) + rho_v(i,j,k) + rho_v(i-1,j+1,k) + rho_v(i-1,j,k)); Real rho_w_loc = 0.25 * (rho_w(i,j,k+1) + rho_w(i,j,k) + rho_w(i,j-1,k+1) + rho_w(i,j-1,k)); rho_u_rhs(i, j, k) += coriolis_factor * (rho_v_loc * sinphi - rho_w_loc * cosphi); } // Add Rayleigh damping - if (use_rayleigh_damping && rayleigh_damp_U) - { + 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) -= tau[k] * (uu - ubar[k]) * rho_on_u_face; } @@ -1032,7 +1043,8 @@ void erf_slow_rhs_pre (int level, int finest_level, { BL_PROFILE("slow_rhs_pre_ymom"); - auto rayleigh_damp_V = solverChoice.rayleigh_damp_V; + auto rayleigh_damp_V = solverChoice.rayleigh_damp_V; + auto geostrophic_wind = solverChoice.custom_geostrophic_profile; // ***************************************************************************** // TERRAIN VERSION // ***************************************************************************** @@ -1070,9 +1082,15 @@ void erf_slow_rhs_pre (int level, int finest_level, q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i,j-1,k,PrimQ1_comp) +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_on_v_face * abl_geo_forcing[1]; + // Add geostrophic wind + if (geostrophic_wind) { + rho_v_rhs(i, j, k) += rho_on_v_face * dptr_v_geos[k]; + } + // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) { Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j-1,k) + rho_u(i,j-1,k)); @@ -1080,8 +1098,7 @@ void erf_slow_rhs_pre (int level, int finest_level, } // Add Rayleigh damping - if (use_rayleigh_damping && rayleigh_damp_V) - { + if (use_rayleigh_damping && rayleigh_damp_V) { 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; } @@ -1113,16 +1130,19 @@ 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_on_v_face * abl_geo_forcing[1]; + // Add geostrophic wind + if (geostrophic_wind) { + rho_v_rhs(i, j, k) += rho_on_v_face * dptr_v_geos[k]; + } + // Add Coriolis forcing (that assumes east is +x, north is +y) - if (use_coriolis) - { + if (use_coriolis) { Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j-1,k) + rho_u(i,j-1,k)); rho_v_rhs(i, j, k) += -coriolis_factor * rho_u_loc * sinphi; } // Add Rayleigh damping - if (use_rayleigh_damping && rayleigh_damp_V) - { + if (use_rayleigh_damping && rayleigh_damp_V) { 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; } @@ -1220,15 +1240,13 @@ void erf_slow_rhs_pre (int level, int finest_level, + rho_on_w_face * abl_geo_forcing[2]; // Add Coriolis forcing (that assumes east is +x, north is +y) - if (use_coriolis) - { + if (use_coriolis) { Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j,k-1) + rho_u(i,j,k-1)); rho_w_rhs(i, j, k) += coriolis_factor * rho_u_loc * cosphi; } // Add Rayleigh damping - if (use_rayleigh_damping && rayleigh_damp_W) - { + if (use_rayleigh_damping && rayleigh_damp_W) { 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; } @@ -1257,15 +1275,13 @@ void erf_slow_rhs_pre (int level, int finest_level, + rho_on_w_face * abl_geo_forcing[2]; // Add Coriolis forcing (that assumes east is +x, north is +y) - if (use_coriolis) - { + if (use_coriolis) { Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j,k-1) + rho_u(i,j,k-1)); rho_w_rhs(i, j, k) += coriolis_factor * rho_u_loc * cosphi; } // Add Rayleigh damping - if (use_rayleigh_damping && rayleigh_damp_W) - { + if (use_rayleigh_damping && rayleigh_damp_W) { 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; } diff --git a/Source/TimeIntegration/TI_headers.H b/Source/TimeIntegration/TI_headers.H index c75524b9b..a0623205f 100644 --- a/Source/TimeIntegration/TI_headers.H +++ b/Source/TimeIntegration/TI_headers.H @@ -55,6 +55,8 @@ 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_u_geos, + const amrex::Real* dptr_v_geos, const amrex::Real* dptr_wbar_sub, const amrex::Vector d_rayleigh_ptrs_at_lev); @@ -268,6 +270,8 @@ 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_u_geos, + const amrex::Real* dptr_v_geos, const amrex::Real* dptr_wbar_sub, const Vector d_rayleigh_ptrs_at_lev); #endif diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 98ce11834..a65cb1184 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -110,7 +110,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_wbar_sub, d_rayleigh_ptrs_at_lev); + dptr_rhotheta_src, dptr_u_geos, dptr_v_geos, dptr_wbar_sub, 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 @@ -203,7 +203,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_wbar_sub, d_rayleigh_ptrs_at_lev); + dptr_rhotheta_src, dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev); } #ifdef ERF_USE_NETCDF @@ -363,7 +363,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_wbar_sub, d_rayleigh_ptrs_at_lev); + dptr_rhotheta_src, dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev); #ifdef ERF_USE_NETCDF diff --git a/Source/prob_common.H b/Source/prob_common.H index 12ef4ac61..4a4c6584b 100644 --- a/Source/prob_common.H +++ b/Source/prob_common.H @@ -216,6 +216,48 @@ public: } /** + * Function to update user-specified temperature source terms that can + * vary with time and height. + * + * @param[in] time current time + * @param[out] u_geos geostrophic wind profile + * @param[out] v_geos geostrophic wind profile + * @param[in] geom container for geometric information + * @param[in] z_phys_cc height coordinate at cell centers + */ + virtual void + update_geostrophic_profile (const amrex::Real& /*time*/, + amrex::Vector& u_geos, + amrex::Gpu::DeviceVector& d_u_geos, + amrex::Vector& v_geos, + amrex::Gpu::DeviceVector& d_v_geos, + const amrex::Geometry& geom, + std::unique_ptr& z_phys_cc) + { + if (u_geos.empty()) return; + + if (z_phys_cc) { + // use_terrain=1 + amrex::Error("Geostrophic wind profile 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(); + for (int k = 0; k <= khi; k++) + { + // 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 + u_geos[k] = 0.0; + v_geos[k] = 0.0; + } + } + + // Copy from host version to device version + amrex::Gpu::copy(amrex::Gpu::hostToDevice, u_geos.begin(), u_geos.end(), d_u_geos.begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, v_geos.begin(), v_geos.end(), d_v_geos.begin()); + } + + /** * Function to perform custom initialization of terrain * * Note: Terrain functionality can also be used to provide grid stretching.