From 17472ad4b93349ca8798c2472f72565fefdc8d93 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Thu, 9 May 2024 16:35:57 -0700 Subject: [PATCH 01/20] SAM Newton Iters (#1611) * update pressure in Newton iterations. * Correction for when qt= tbgmax) { + if(tabs_array(i,j,k) >= tbgmax) { omn = 1.0; delta_qi = qci_array(i,j,k); qci_array(i,j,k) = 0.0; @@ -71,7 +67,7 @@ void SAM::Cloud () { pres_array(i,j,k) /= 100.0; } // Cloud water not permitted (freeze to form ice) - else if(tabs <= tbgmin) { + else if(tabs_array(i,j,k) <= tbgmin) { omn = 0.0; delta_qc = qcl_array(i,j,k); qcl_array(i,j,k) = 0.0; @@ -84,7 +80,7 @@ void SAM::Cloud () { } // Mixed cloud phase (split according to omn) else { - omn = an*tabs-bn; + omn = an*tabs_array(i,j,k)-bn; delta_qc = qcl_array(i,j,k) - qn_array(i,j,k) * omn; delta_qi = qci_array(i,j,k) - qn_array(i,j,k) * (1.0 - omn); qcl_array(i,j,k) = qn_array(i,j,k) * omn; @@ -96,6 +92,10 @@ void SAM::Cloud () { pres_array(i,j,k) /= 100.0; } + // Initial guess for temperature & pressure + Real tabs = tabs_array(i,j,k); + Real pres = pres_array(i,j,k); + // Saturation moisture fractions erf_qsatw(tabs, pres, qsatw); erf_qsati(tabs, pres, qsati); @@ -152,6 +152,13 @@ void SAM::Cloud () { // Update the temperature dtabs = -fff/dfff; tabs = tabs+dtabs; + + // Update the pressure + pres = rho_array(i,j,k) * R_d * tabs + * (1.0 + R_v/R_d * qsat); + pres /= 100.0; + + // Update iteration niter = niter+1; } while(std::abs(dtabs) > tol && niter < 20); @@ -183,7 +190,34 @@ void SAM::Cloud () { // Pressure unit conversion pres_array(i,j,k) /= 100.0; - } // qt > qsat + } + // We cannot blindly relax to qsat, but we can convert qc/qi -> qv + else { + // Changes in each component + delta_qv = qcl_array(i,j,k) + qci_array(i,j,k); + delta_qc = -qcl_array(i,j,k); + delta_qi = -qci_array(i,j,k); + + // Partition the change in non-precipitating q + qv_array(i,j,k) += delta_qv; + qcl_array(i,j,k) = 0.0; + qci_array(i,j,k) = 0.0; + qn_array(i,j,k) = 0.0; + qt_array(i,j,k) = qv_array(i,j,k); + + // Update temperature (endothermic since we evap/sublime) + tabs_array(i,j,k) -= fac_fus * delta_qc + fac_sub * delta_qi; + + // Update pressure + pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) + * (1.0 + R_v/R_d * qv_array(i,j,k)); + + // Update theta from temperature + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); + + // Pressure unit conversion + pres_array(i,j,k) /= 100.0; + } }); } // mfi } From 084af377510ef77a6598c2a380f09e4ec91d2c7f Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Thu, 9 May 2024 18:01:37 -0700 Subject: [PATCH 02/20] enable incompressible to work with multilevel but only OneWay coupling (#1612) * enable incompressible to work with multilevel but only OneWay coupling * code cleanup in ERF_slow_rhs_pre and ERF_slow_rhs_inc * add file to cmake * more changes to sync inc with pre --------- Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> --- CMake/BuildERFExe.cmake | 1 + Exec/Make.ERF | 1 + Source/ERF.H | 16 +- Source/ERF.cpp | 18 +- Source/ERF_make_new_arrays.cpp | 5 + Source/ERF_make_new_level.cpp | 4 + Source/IO/Plotfile.cpp | 16 +- Source/TimeIntegration/ERF_make_tau_terms.cpp | 419 ++++++++++++++++++ Source/TimeIntegration/ERF_slow_rhs_inc.cpp | 359 ++------------- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 329 +------------- Source/TimeIntegration/Make.package | 1 + Source/TimeIntegration/TI_no_substep_fun.H | 4 +- Source/TimeIntegration/TI_slow_headers.H | 36 +- Source/TimeIntegration/TI_slow_rhs_fun.H | 8 +- Source/Utils/ERF_PoissonSolve.cpp | 169 ++++--- Source/Utils/ERF_PoissonSolve_tb.cpp | 337 ++++++-------- Submodules/AMReX | 2 +- 17 files changed, 777 insertions(+), 948 deletions(-) create mode 100644 Source/TimeIntegration/ERF_make_tau_terms.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 6f4bae92e..412c1d926 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -165,6 +165,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/TimeIntegration/ERF_advance_lsm.cpp ${SRC_DIR}/TimeIntegration/ERF_advance_radiation.cpp ${SRC_DIR}/TimeIntegration/ERF_make_fast_coeffs.cpp + ${SRC_DIR}/TimeIntegration/ERF_make_tau_terms.cpp ${SRC_DIR}/TimeIntegration/ERF_slow_rhs_pre.cpp ${SRC_DIR}/TimeIntegration/ERF_slow_rhs_post.cpp ${SRC_DIR}/TimeIntegration/ERF_fast_rhs_N.cpp diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 82311228b..1b52fbb37 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -177,6 +177,7 @@ endif ifeq ($(USE_POISSON_SOLVE), TRUE) DEFINES += -DERF_USE_POISSON_SOLVE + USERSuffix += .INC endif ifeq ($(USE_PARTICLES), TRUE) diff --git a/Source/ERF.H b/Source/ERF.H index 1d44e73ad..5beb392c3 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -123,17 +123,16 @@ public: #ifdef ERF_USE_POISSON_SOLVE // Project the velocities to be divergence-free - void project_velocities ( amrex::Vector& vars); - void project_velocities (int lev_min, int lev_max, amrex::Vector>& vars); + void project_velocities (int lev, amrex::Real dt, amrex::Vector& vars, amrex::MultiFab& p); + + // Project the velocities to be divergence-free with a thin body + void project_velocities_tb (int lev, amrex::Real dt, amrex::Vector& vars, amrex::MultiFab& p); // Define the projection bc's based on the domain bc types amrex::Array get_projection_bc (amrex::Orientation::Side side) const noexcept; bool projection_has_dirichlet (amrex::Array bcs) const; - // Project the velocities to be divergence-free with a thin body - void project_velocities_tb ( amrex::Vector& vars, const amrex::Real dt); - void project_velocities_tb (amrex::Vector>& vars, const amrex::Real dt); #endif // Init (NOT restart or regrid) @@ -605,6 +604,10 @@ private: #endif amrex::Vector > > > mri_integrator_mem; +#ifdef ERF_USE_POISSON_SOLVE + amrex::Vector pp_inc; +#endif + // Vector over levels of routines to impose physical boundary conditions amrex::Vector> physbcs_cons; amrex::Vector> physbcs_u; @@ -840,6 +843,9 @@ private: #endif static int verbose; +#ifdef ERF_USE_POISSON_SOLVE + static int mg_verbose; +#endif // Diagnostic output interval static int sum_interval; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index e1991f4e9..c6215973f 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -38,6 +38,9 @@ int ERF::fixed_mri_dt_ratio = 0; // Dictate verbosity in screen output int ERF::verbose = 0; +#ifdef ERF_USE_POISSON_SOLVE +int ERF::mg_verbose = 0; +#endif // Frequency of diagnostic output int ERF::sum_interval = -1; @@ -182,6 +185,11 @@ ERF::ERF () vars_new.resize(nlevs_max); vars_old.resize(nlevs_max); +#ifdef ERF_USE_POISSON_SOLVE + pp_inc.resize(nlevs_max); +#endif + + rU_new.resize(nlevs_max); rV_new.resize(nlevs_max); rW_new.resize(nlevs_max); @@ -819,7 +827,12 @@ ERF::InitData () // Note -- this projection is only defined for no terrain if (solverChoice.project_initial_velocity) { AMREX_ALWAYS_ASSERT(solverChoice.use_terrain == 0); - project_velocities(0, finest_level, vars_new); + Real dummy_dt = 1.0; + for (int lev = 0; lev <= finest_level; ++lev) + { + project_velocities(lev, dummy_dt, vars_new[lev], pp_inc[lev]); + pp_inc[lev].setVal(0.); + } } } #endif @@ -1228,6 +1241,9 @@ ERF::ReadParameters () // Verbosity pp.query("v", verbose); +#ifdef ERF_USE_POISSON_SOLVE + pp.query("mg_v", mg_verbose); +#endif // Frequency of diagnostic output pp.query("sum_interval", sum_interval); diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 13102fc22..8dbda2bff 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -141,6 +141,11 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, lev_new[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0)); lev_old[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0)); +#ifdef ERF_USE_POISSON_SOLVE + pp_inc[lev].define(ba, dm, 1, 1); + pp_inc[lev].setVal(0.0); +#endif + // ******************************************************************************************** // These are just used for scratch in the time integrator but we might as well define them here // ******************************************************************************************** diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 43f8c7dff..d06f55752 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -499,6 +499,10 @@ ERF::ClearLevel (int lev) rW_new[lev].clear(); rW_old[lev].clear(); +#ifdef ERF_USE_POISSON_SOLVE + pp_inc[lev].clear(); +#endif + // Clears the integrator memory mri_integrator_mem[lev].reset(); diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index e486d5314..1ca5db770 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -80,16 +80,7 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector } #endif -#ifdef ERF_USE_EB - if (containerHasElement(plot_var_names, "volfrac")) { - tmp_plot_names.push_back("volfrac"); - } -#endif - plot_var_names = tmp_plot_names; - - for (int i=0; i plot_var_names) } #endif +#ifdef ERF_USE_POISSON_SOLVE + if (containerHasElement(plot_var_names, "pp_inc")) { + MultiFab::Copy(mf[lev], pp_inc[lev], 0, mf_comp, 1, 0); + mf_comp += 1; + } +#endif + #ifdef ERF_COMPUTE_ERROR // Next, check for error in velocities and if desired, output them -- note we output none or all, not just some if (containerHasElement(plot_var_names, "xvel_err") || diff --git a/Source/TimeIntegration/ERF_make_tau_terms.cpp b/Source/TimeIntegration/ERF_make_tau_terms.cpp new file mode 100644 index 000000000..456697719 --- /dev/null +++ b/Source/TimeIntegration/ERF_make_tau_terms.cpp @@ -0,0 +1,419 @@ +#include +#include +#include +#include + +#include +#include +#include + +using namespace amrex; + +void erf_make_tau_terms (int level, int nrk, + const Gpu::DeviceVector& domain_bcs_type_d, + const Vector& domain_bcs_type_h, + std::unique_ptr& z_phys_nd, + Vector& S_data, + const MultiFab& xvel, + const MultiFab& yvel, + const MultiFab& zvel, + MultiFab& Omega, + MultiFab* Tau11, MultiFab* Tau22, MultiFab* Tau33, + MultiFab* Tau12, MultiFab* Tau13, MultiFab* Tau21, + MultiFab* Tau23, MultiFab* Tau31, MultiFab* Tau32, + MultiFab* SmnSmn, + MultiFab* eddyDiffs, + MultiFab* Hfx3, MultiFab* Diss, + const Geometry geom, + const SolverChoice& solverChoice, + std::unique_ptr& most, + std::unique_ptr& detJ, + std::unique_ptr& mapfac_m, + std::unique_ptr& mapfac_u, + std::unique_ptr& mapfac_v) +{ + BL_PROFILE_REGION("erf_make_tau_terms()"); + + const BCRec* bc_ptr_d = domain_bcs_type_d.data(); + const BCRec* bc_ptr_h = domain_bcs_type_h.data(); + + DiffChoice dc = solverChoice.diffChoice; + TurbChoice tc = solverChoice.turbChoice[level]; + + int start_comp = 0; + int num_comp = 2; + int end_comp = start_comp + num_comp - 1; + + const AdvType l_horiz_adv_type = solverChoice.advChoice.dycore_horiz_adv_type; + const AdvType l_vert_adv_type = solverChoice.advChoice.dycore_vert_adv_type; + const Real l_horiz_upw_frac = solverChoice.advChoice.dycore_horiz_upw_frac; + const Real l_vert_upw_frac = solverChoice.advChoice.dycore_vert_upw_frac; + const bool l_use_terrain = solverChoice.use_terrain; + const bool l_moving_terrain = (solverChoice.terrain_type == TerrainType::Moving); + if (l_moving_terrain) AMREX_ALWAYS_ASSERT (l_use_terrain); + + const bool l_reflux = (solverChoice.coupling_type == CouplingType::TwoWay); + + const bool l_use_diff = ( (dc.molec_diff_type != MolecDiffType::None) || + (tc.les_type != LESType::None) || + (tc.pbl_type != PBLType::None) ); + const bool l_use_constAlpha = ( dc.molec_diff_type == MolecDiffType::ConstantAlpha ); + const bool l_use_turb = ( tc.les_type == LESType::Smagorinsky || + tc.les_type == LESType::Deardorff || + tc.pbl_type == PBLType::MYNN25 || + tc.pbl_type == PBLType::YSU ); + + const bool use_moisture = (solverChoice.moisture_type != MoistureType::None); + const bool use_most = (most != nullptr); + const bool exp_most = (solverChoice.use_explicit_most); + + const Box& domain = geom.Domain(); + const int domhi_z = domain.bigEnd(2); + const int domlo_z = domain.smallEnd(2); + + const GpuArray dxInv = geom.InvCellSizeArray(); + const Real* dx = geom.CellSize(); + + // ***************************************************************************** + // Combine external forcing terms + // ***************************************************************************** + const Array grav{0.0, 0.0, -solverChoice.gravity}; + const GpuArray grav_gpu{grav[0], grav[1], grav[2]}; + + // ***************************************************************************** + // Pre-computed quantities + // ***************************************************************************** + int nvars = S_data[IntVars::cons].nComp(); + const BoxArray& ba = S_data[IntVars::cons].boxArray(); + const DistributionMapping& dm = S_data[IntVars::cons].DistributionMap(); + + std::unique_ptr expr; + std::unique_ptr dflux_x; + std::unique_ptr dflux_y; + std::unique_ptr dflux_z; + + if (l_use_diff) { + expr = std::make_unique(ba , dm, 1, IntVect(1,1,0)); + dflux_x = std::make_unique(convert(ba,IntVect(1,0,0)), dm, nvars, 0); + dflux_y = std::make_unique(convert(ba,IntVect(0,1,0)), dm, nvars, 0); + dflux_z = std::make_unique(convert(ba,IntVect(0,0,1)), dm, nvars, 0); + + // if using constant alpha (mu = rho * alpha), then first divide by the + // reference density -- mu_eff will be scaled by the instantaneous + // local density later when ComputeStress*Visc_*() is called + Real mu_eff = (l_use_constAlpha) ? 2.0 * dc.dynamicViscosity / dc.rho0_trans + : 2.0 * dc.dynamicViscosity; + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Box& valid_bx = mfi.validbox(); + + // Velocities + const Array4 & u = xvel.array(mfi); + const Array4 & v = yvel.array(mfi); + const Array4 & w = zvel.array(mfi); + + // Contravariant velocity + const Array4& omega_arr = Omega.array(mfi); + + // Map factors + const Array4& mf_m = mapfac_m->const_array(mfi); + const Array4& mf_u = mapfac_u->const_array(mfi); + const Array4& mf_v = mapfac_v->const_array(mfi); + + // Eddy viscosity + const Array4& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : Array4{}; + const Array4& cell_data = l_use_constAlpha ? S_data[IntVars::cons].const_array(mfi) : Array4{}; + + // Terrain metrics + const Array4& z_nd = l_use_terrain ? z_phys_nd->const_array(mfi) : Array4{}; + const Array4& detJ_arr = detJ->const_array(mfi); + + //------------------------------------------------------------------------------- + // NOTE: Tile boxes with terrain are not intuitive. The linear combination of + // stress terms requires care. Create a tile box that intersects the + // valid box, then grow the box in x/y. Compute the strain on the local + // FAB over this grown tile box. Compute the stress over the tile box, + // except tau_ii which still needs the halo cells. Finally, write from + // the local FAB to the Tau MF but only on the tile box. + //------------------------------------------------------------------------------- + + //------------------------------------------------------------------------------- + // TODO: Avoid recomputing strain on the first RK stage. One could populate + // the FABs with tau_ij, compute stress, and then write to tau_ij. The + // problem with this approach is you will over-write the needed halo layer + // needed by subsequent tile boxes (particularly S_ii becomes Tau_ii). + //------------------------------------------------------------------------------- + + // Strain/Stress tile boxes + Box bxcc = mfi.tilebox(); + Box tbxxy = mfi.tilebox(IntVect(1,1,0)); + Box tbxxz = mfi.tilebox(IntVect(1,0,1)); + Box tbxyz = mfi.tilebox(IntVect(0,1,1)); + + // We need a halo cell for terrain + bxcc.grow(IntVect(1,1,0)); + tbxxy.grow(IntVect(1,1,0)); + tbxxz.grow(IntVect(1,1,0)); + tbxyz.grow(IntVect(1,1,0)); + + // Expansion rate + Array4 er_arr = expr->array(mfi); + + // Temporary storage for tiling/OMP + FArrayBox S11,S22,S33; + FArrayBox S12,S13,S23; + S11.resize( bxcc,1,The_Async_Arena()); S22.resize(bxcc,1,The_Async_Arena()); S33.resize(bxcc,1,The_Async_Arena()); + S12.resize(tbxxy,1,The_Async_Arena()); S13.resize(tbxxz,1,The_Async_Arena()); S23.resize(tbxyz,1,The_Async_Arena()); + Array4 s11 = S11.array(); Array4 s22 = S22.array(); Array4 s33 = S33.array(); + Array4 s12 = S12.array(); Array4 s13 = S13.array(); Array4 s23 = S23.array(); + + // Symmetric strain/stresses + Array4 tau11 = Tau11->array(mfi); Array4 tau22 = Tau22->array(mfi); Array4 tau33 = Tau33->array(mfi); + Array4 tau12 = Tau12->array(mfi); Array4 tau13 = Tau13->array(mfi); Array4 tau23 = Tau23->array(mfi); + + // Strain magnitude + Array4 SmnSmn_a; + + if (l_use_terrain) { + // Terrain non-symmetric terms + FArrayBox S21,S31,S32; + S21.resize(tbxxy,1,The_Async_Arena()); S31.resize(tbxxz,1,The_Async_Arena()); S32.resize(tbxyz,1,The_Async_Arena()); + Array4 s21 = S21.array(); Array4 s31 = S31.array(); Array4 s32 = S32.array(); + Array4 tau21 = Tau21->array(mfi); Array4 tau31 = Tau31->array(mfi); Array4 tau32 = Tau32->array(mfi); + + // ***************************************************************************** + // Expansion rate compute terrain + // ***************************************************************************** + { + BL_PROFILE("slow_rhs_making_er_T"); + // First create Omega using velocity (not momentum) + Box gbxo = surroundingNodes(bxcc,2); + ParallelFor(gbxo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + omega_arr(i,j,k) = (k == 0) ? 0. : OmegaFromW(i,j,k,w(i,j,k),u,v,z_nd,dxInv); + }); + + ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + + Real met_u_h_zeta_hi = Compute_h_zeta_AtIface(i+1, j , k, dxInv, z_nd); + Real met_u_h_zeta_lo = Compute_h_zeta_AtIface(i , j , k, dxInv, z_nd); + + Real met_v_h_zeta_hi = Compute_h_zeta_AtJface(i , j+1, k, dxInv, z_nd); + Real met_v_h_zeta_lo = Compute_h_zeta_AtJface(i , j , k, dxInv, z_nd); + + Real Omega_hi = omega_arr(i,j,k+1); + Real Omega_lo = omega_arr(i,j,k ); + + Real mfsq = mf_m(i,j,0)*mf_m(i,j,0); + + Real expansionRate = (u(i+1,j ,k)/mf_u(i+1,j,0)*met_u_h_zeta_hi - u(i,j,k)/mf_u(i,j,0)*met_u_h_zeta_lo)*dxInv[0]*mfsq + + (v(i ,j+1,k)/mf_v(i,j+1,0)*met_v_h_zeta_hi - v(i,j,k)/mf_v(i,j,0)*met_v_h_zeta_lo)*dxInv[1]*mfsq + + (Omega_hi - Omega_lo)*dxInv[2]; + + er_arr(i,j,k) = expansionRate / detJ_arr(i,j,k); + }); + } // end profile + + // ***************************************************************************** + // Strain tensor compute terrain + // ***************************************************************************** + { + BL_PROFILE("slow_rhs_making_strain_T"); + ComputeStrain_T(bxcc, tbxxy, tbxxz, tbxyz, domain, + u, v, w, + s11, s22, s33, + s12, s13, + s21, s23, + s31, s32, + z_nd, detJ_arr, bc_ptr_h, dxInv, + mf_m, mf_u, mf_v); + } // profile + + // Populate SmnSmn if using Deardorff (used as diff src in post) + // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF) + if ((nrk==0) && (tc.les_type == LESType::Deardorff)) { + SmnSmn_a = SmnSmn->array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most,exp_most); + }); + } + + // We've updated the strains at all locations including the + // surface. This is required to get the correct strain-rate + // magnitude. Now, update the stress everywhere but the surface + // to retain the values set by MOST. + if (use_most && exp_most) { + // Don't overwrite modeled total stress value at boundary + tbxxz.setSmall(2,1); + tbxyz.setSmall(2,1); + } + + // ***************************************************************************** + // Stress tensor compute terrain + // ***************************************************************************** + { + BL_PROFILE("slow_rhs_making_stress_T"); + + // Remove Halo cells just for tau_ij comps + tbxxy.grow(IntVect(-1,-1,0)); + tbxxz.grow(IntVect(-1,-1,0)); + tbxyz.grow(IntVect(-1,-1,0)); + + if (!l_use_turb) { + ComputeStressConsVisc_T(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, + cell_data, + s11, s22, s33, + s12, s13, + s21, s23, + s31, s32, + er_arr, z_nd, detJ_arr, dxInv); + } else { + ComputeStressVarVisc_T(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, mu_turb, + cell_data, + s11, s22, s33, + s12, s13, + s21, s23, + s31, s32, + er_arr, z_nd, detJ_arr, dxInv); + } + + // Remove halo cells from tau_ii but extend across valid_box bdry + bxcc.grow(IntVect(-1,-1,0)); + if (bxcc.smallEnd(0) == valid_bx.smallEnd(0)) bxcc.growLo(0, 1); + if (bxcc.bigEnd(0) == valid_bx.bigEnd(0)) bxcc.growHi(0, 1); + if (bxcc.smallEnd(1) == valid_bx.smallEnd(1)) bxcc.growLo(1, 1); + if (bxcc.bigEnd(1) == valid_bx.bigEnd(1)) bxcc.growHi(1, 1); + + // Copy from temp FABs back to tau + ParallelFor(bxcc, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + tau11(i,j,k) = s11(i,j,k); + tau22(i,j,k) = s22(i,j,k); + tau33(i,j,k) = s33(i,j,k); + }); + + ParallelFor(tbxxy, tbxxz, tbxyz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + tau12(i,j,k) = s12(i,j,k); + tau21(i,j,k) = s21(i,j,k); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + tau13(i,j,k) = s13(i,j,k); + tau31(i,j,k) = s31(i,j,k); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + tau23(i,j,k) = s23(i,j,k); + tau32(i,j,k) = s32(i,j,k); + }); + } // end profile + + } else { + + // ***************************************************************************** + // Expansion rate compute no terrain + // ***************************************************************************** + { + BL_PROFILE("slow_rhs_making_er_N"); + ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + Real mfsq = mf_m(i,j,0)*mf_m(i,j,0); + er_arr(i,j,k) = (u(i+1, j , k )/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0]*mfsq + + (v(i , j+1, k )/mf_v(i,j+1,0) - v(i, j, k)/mf_v(i,j,0))*dxInv[1]*mfsq + + (w(i , j , k+1) - w(i, j, k))*dxInv[2]; + }); + } // end profile + + + // ***************************************************************************** + // Strain tensor compute no terrain + // ***************************************************************************** + { + BL_PROFILE("slow_rhs_making_strain_N"); + ComputeStrain_N(bxcc, tbxxy, tbxxz, tbxyz, domain, + u, v, w, + s11, s22, s33, + s12, s13, s23, + bc_ptr_h, dxInv, + mf_m, mf_u, mf_v); + } // end profile + + // Populate SmnSmn if using Deardorff (used as diff src in post) + // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF) + if ((nrk==0) && (tc.les_type == LESType::Deardorff)) { + SmnSmn_a = SmnSmn->array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most,exp_most); + }); + } + + // We've updated the strains at all locations including the + // surface. This is required to get the correct strain-rate + // magnitude. Now, update the stress everywhere but the surface + // to retain the values set by MOST. + if (use_most && exp_most) { + // Don't overwrite modeled total stress value at boundary + tbxxz.setSmall(2,1); + tbxyz.setSmall(2,1); + } + + // ***************************************************************************** + // Stress tensor compute no terrain + // ***************************************************************************** + { + BL_PROFILE("slow_rhs_making_stress_N"); + + // Remove Halo cells just for tau_ij comps + tbxxy.grow(IntVect(-1,-1,0)); + tbxxz.grow(IntVect(-1,-1,0)); + tbxyz.grow(IntVect(-1,-1,0)); + + if (!l_use_turb) { + ComputeStressConsVisc_N(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, + cell_data, + s11, s22, s33, + s12, s13, s23, + er_arr); + } else { + ComputeStressVarVisc_N(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, mu_turb, + cell_data, + s11, s22, s33, + s12, s13, s23, + er_arr); + } + + // Remove halo cells from tau_ii but extend across valid_box bdry + bxcc.grow(IntVect(-1,-1,0)); + if (bxcc.smallEnd(0) == valid_bx.smallEnd(0)) bxcc.growLo(0, 1); + if (bxcc.bigEnd(0) == valid_bx.bigEnd(0)) bxcc.growHi(0, 1); + if (bxcc.smallEnd(1) == valid_bx.smallEnd(1)) bxcc.growLo(1, 1); + if (bxcc.bigEnd(1) == valid_bx.bigEnd(1)) bxcc.growHi(1, 1); + + // Copy from temp FABs back to tau + ParallelFor(bxcc, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + tau11(i,j,k) = s11(i,j,k); + tau22(i,j,k) = s22(i,j,k); + tau33(i,j,k) = s33(i,j,k); + }); + ParallelFor(tbxxy, tbxxz, tbxyz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + tau12(i,j,k) = s12(i,j,k); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + tau13(i,j,k) = s13(i,j,k); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + tau23(i,j,k) = s23(i,j,k); + }); + } // end profile + } // l_use_terrain + } // MFIter + } // l_use_diff +} diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index c6df98725..45a096757 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -55,6 +55,7 @@ using namespace amrex; * @param[in] mapfac_m map factor at cell centers * @param[in] mapfac_u map factor at x-faces * @param[in] mapfac_v map factor at y-faces + * @param[in] pp_inc perturbational pressure for incompressible flows */ void erf_slow_rhs_inc (int level, int nrk, @@ -89,6 +90,7 @@ void erf_slow_rhs_inc (int level, int nrk, std::unique_ptr& az, std::unique_ptr& detJ, const MultiFab* p0, + const MultiFab& pp_inc, std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v) @@ -154,333 +156,17 @@ void erf_slow_rhs_inc (int level, int nrk, std::unique_ptr dflux_z; if (l_use_diff) { - expr = std::make_unique(ba , dm, 1, IntVect(1,1,0)); + erf_make_tau_terms(level,nrk,domain_bcs_type_d,domain_bcs_type_h,z_phys_nd, + S_data,xvel,yvel,zvel,Omega, + Tau11,Tau22,Tau33,Tau12,Tau13,Tau21,Tau23,Tau31,Tau32, + SmnSmn,eddyDiffs,Hfx3,Diss,geom,solverChoice,most, + detJ,mapfac_m,mapfac_u,mapfac_v); + dflux_x = std::make_unique(convert(ba,IntVect(1,0,0)), dm, nvars, 0); dflux_y = std::make_unique(convert(ba,IntVect(0,1,0)), dm, nvars, 0); dflux_z = std::make_unique(convert(ba,IntVect(0,0,1)), dm, nvars, 0); - - // if using constant alpha (mu = rho * alpha), then first divide by the - // reference density -- mu_eff will be scaled by the instantaneous - // local density later when ComputeStress*Visc_*() is called - Real mu_eff = (l_use_constAlpha) ? 2.0 * dc.dynamicViscosity / dc.rho0_trans - : 2.0 * dc.dynamicViscosity; - -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); - const Box& valid_bx = mfi.validbox(); - - // Velocities - const Array4 & u = xvel.array(mfi); - const Array4 & v = yvel.array(mfi); - const Array4 & w = zvel.array(mfi); - - // Contravariant velocity - const Array4& omega_arr = Omega.array(mfi); - - // Map factors - const Array4& mf_m = mapfac_m->const_array(mfi); - const Array4& mf_u = mapfac_u->const_array(mfi); - const Array4& mf_v = mapfac_v->const_array(mfi); - - // Eddy viscosity - const Array4& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : Array4{}; - const Array4& cell_data = l_use_constAlpha ? S_data[IntVars::cons].const_array(mfi) : Array4{}; - - // Terrain metrics - const Array4& z_nd = l_use_terrain ? z_phys_nd->const_array(mfi) : Array4{}; - const Array4& detJ_arr = detJ->const_array(mfi); - - //------------------------------------------------------------------------------- - // NOTE: Tile boxes with terrain are not intuitive. The linear combination of - // stress terms requires care. Create a tile box that intersects the - // valid box, then grow the box in x/y. Compute the strain on the local - // FAB over this grown tile box. Compute the stress over the tile box, - // except tau_ii which still needs the halo cells. Finally, write from - // the local FAB to the Tau MF but only on the tile box. - //------------------------------------------------------------------------------- - - //------------------------------------------------------------------------------- - // TODO: Avoid recomputing strain on the first RK stage. One could populate - // the FABs with tau_ij, compute stress, and then write to tau_ij. The - // problem with this approach is you will over-write the needed halo layer - // needed by subsequent tile boxes (particularly S_ii becomes Tau_ii). - //------------------------------------------------------------------------------- - - // Strain/Stress tile boxes - Box bxcc = mfi.tilebox(); - Box tbxxy = mfi.tilebox(IntVect(1,1,0)); - Box tbxxz = mfi.tilebox(IntVect(1,0,1)); - Box tbxyz = mfi.tilebox(IntVect(0,1,1)); - - // We need a halo cell for terrain - bxcc.grow(IntVect(1,1,0)); - tbxxy.grow(IntVect(1,1,0)); - tbxxz.grow(IntVect(1,1,0)); - tbxyz.grow(IntVect(1,1,0)); - - // Expansion rate - Array4 er_arr = expr->array(mfi); - - // Temporary storage for tiling/OMP - FArrayBox S11,S22,S33; - FArrayBox S12,S13,S23; - S11.resize( bxcc,1,The_Async_Arena()); S22.resize( bxcc,1,The_Async_Arena()); S33.resize( bxcc,1,The_Async_Arena()); - S12.resize(tbxxy,1,The_Async_Arena()); S13.resize(tbxxz,1,The_Async_Arena()); S23.resize(tbxyz,1,The_Async_Arena()); - Array4 s11 = S11.array(); Array4 s22 = S22.array(); Array4 s33 = S33.array(); - Array4 s12 = S12.array(); Array4 s13 = S13.array(); Array4 s23 = S23.array(); - - // Symmetric strain/stresses - Array4 tau11 = Tau11->array(mfi); Array4 tau22 = Tau22->array(mfi); Array4 tau33 = Tau33->array(mfi); - Array4 tau12 = Tau12->array(mfi); Array4 tau13 = Tau13->array(mfi); Array4 tau23 = Tau23->array(mfi); - - // Strain magnitude - Array4 SmnSmn_a; - - if (l_use_terrain) { - // Terrain non-symmetric terms - FArrayBox S21,S31,S32; - S21.resize(tbxxy,1,The_Async_Arena()); S31.resize(tbxxz,1,The_Async_Arena()); S32.resize(tbxyz,1,The_Async_Arena()); - Array4 s21 = S21.array(); Array4 s31 = S31.array(); Array4 s32 = S32.array(); - Array4 tau21 = Tau21->array(mfi); Array4 tau31 = Tau31->array(mfi); Array4 tau32 = Tau32->array(mfi); - - // ***************************************************************************** - // Expansion rate compute terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_er_T"); - // First create Omega using velocity (not momentum) - Box gbxo = surroundingNodes(bxcc,2); - ParallelFor(gbxo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - omega_arr(i,j,k) = (k == 0) ? 0. : OmegaFromW(i,j,k,w(i,j,k),u,v,z_nd,dxInv); - }); - - ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - - Real met_u_h_zeta_hi = Compute_h_zeta_AtIface(i+1, j , k, dxInv, z_nd); - Real met_u_h_zeta_lo = Compute_h_zeta_AtIface(i , j , k, dxInv, z_nd); - - Real met_v_h_zeta_hi = Compute_h_zeta_AtJface(i , j+1, k, dxInv, z_nd); - Real met_v_h_zeta_lo = Compute_h_zeta_AtJface(i , j , k, dxInv, z_nd); - - Real Omega_hi = omega_arr(i,j,k+1); - Real Omega_lo = omega_arr(i,j,k ); - - Real mfsq = mf_m(i,j,0)*mf_m(i,j,0); - - Real expansionRate = (u(i+1,j ,k)/mf_u(i+1,j,0)*met_u_h_zeta_hi - u(i,j,k)/mf_u(i,j,0)*met_u_h_zeta_lo)*dxInv[0]*mfsq + - (v(i ,j+1,k)/mf_v(i,j+1,0)*met_v_h_zeta_hi - v(i,j,k)/mf_v(i,j,0)*met_v_h_zeta_lo)*dxInv[1]*mfsq + - (Omega_hi - Omega_lo)*dxInv[2]; - - er_arr(i,j,k) = expansionRate / detJ_arr(i,j,k); - }); - } // end profile - - // ***************************************************************************** - // Strain tensor compute terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_strain_T"); - ComputeStrain_T(bxcc, tbxxy, tbxxz, tbxyz, domain, - u, v, w, - s11, s22, s33, - s12, s13, - s21, s23, - s31, s32, - z_nd, detJ_arr, bc_ptr_h, dxInv, - mf_m, mf_u, mf_v); - } // profile - - // Populate SmnSmn if using Deardorff (used as diff src in post) - // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF) - if ((nrk==0) && (tc.les_type == LESType::Deardorff)) { - SmnSmn_a = SmnSmn->array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most); - }); - } - - // We've updated the strains at all locations including the - // surface. This is required to get the correct strain-rate - // magnitude. Now, update the stress everywhere but the surface - // to retain the values set by MOST. - if (use_most && exp_most) { - // Don't overwrite modeled total stress value at boundary - tbxxz.setSmall(2,1); - tbxyz.setSmall(2,1); - } - - // ***************************************************************************** - // Stress tensor compute terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_stress_T"); - - // Remove Halo cells just for tau_ij comps - tbxxy.grow(IntVect(-1,-1,0)); - tbxxz.grow(IntVect(-1,-1,0)); - tbxyz.grow(IntVect(-1,-1,0)); - - if (!l_use_turb) { - ComputeStressConsVisc_T(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, - cell_data, - s11, s22, s33, - s12, s13, - s21, s23, - s31, s32, - er_arr, z_nd, detJ_arr, dxInv); - } else { - ComputeStressVarVisc_T(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, mu_turb, - cell_data, - s11, s22, s33, - s12, s13, - s21, s23, - s31, s32, - er_arr, z_nd, detJ_arr, dxInv); - } - - // Remove halo cells from tau_ii but extend across valid_box bdry - bxcc.grow(IntVect(-1,-1,0)); - if (bxcc.smallEnd(0) == valid_bx.smallEnd(0)) bxcc.growLo(0, 1); - if (bxcc.bigEnd(0) == valid_bx.bigEnd(0)) bxcc.growHi(0, 1); - if (bxcc.smallEnd(1) == valid_bx.smallEnd(1)) bxcc.growLo(1, 1); - if (bxcc.bigEnd(1) == valid_bx.bigEnd(1)) bxcc.growHi(1, 1); - - // Copy from temp FABs back to tau - ParallelFor(bxcc, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau11(i,j,k) = s11(i,j,k); - tau22(i,j,k) = s22(i,j,k); - tau33(i,j,k) = s33(i,j,k); - }); - - ParallelFor(tbxxy, tbxxz, tbxyz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau12(i,j,k) = s12(i,j,k); - tau21(i,j,k) = s21(i,j,k); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau13(i,j,k) = s13(i,j,k); - tau31(i,j,k) = s31(i,j,k); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau23(i,j,k) = s23(i,j,k); - tau32(i,j,k) = s32(i,j,k); - }); - } // end profile - - } else { - - // ***************************************************************************** - // Expansion rate compute no terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_er_N"); - ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real mfsq = mf_m(i,j,0)*mf_m(i,j,0); - er_arr(i,j,k) = (u(i+1, j , k )/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0]*mfsq + - (v(i , j+1, k )/mf_v(i,j+1,0) - v(i, j, k)/mf_v(i,j,0))*dxInv[1]*mfsq + - (w(i , j , k+1) - w(i, j, k))*dxInv[2]; - }); - } // end profile - - - // ***************************************************************************** - // Strain tensor compute no terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_strain_N"); - ComputeStrain_N(bxcc, tbxxy, tbxxz, tbxyz, domain, - u, v, w, - s11, s22, s33, - s12, s13, s23, - bc_ptr_h, dxInv, - mf_m, mf_u, mf_v); - } // end profile - - // Populate SmnSmn if using Deardorff (used as diff src in post) - // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF) - if ((nrk==0) && (tc.les_type == LESType::Deardorff)) { - SmnSmn_a = SmnSmn->array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most); - }); - } - - // We've updated the strains at all locations including the - // surface. This is required to get the correct strain-rate - // magnitude. Now, update the stress everywhere but the surface - // to retain the values set by MOST. - if (use_most && exp_most) { - // Don't overwrite modeled total stress value at boundary - tbxxz.setSmall(2,1); - tbxyz.setSmall(2,1); - } - - // ***************************************************************************** - // Stress tensor compute no terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_stress_N"); - - // Remove Halo cells just for tau_ij comps - tbxxy.grow(IntVect(-1,-1,0)); - tbxxz.grow(IntVect(-1,-1,0)); - tbxyz.grow(IntVect(-1,-1,0)); - - if (!l_use_turb) { - ComputeStressConsVisc_N(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, - cell_data, - s11, s22, s33, - s12, s13, s23, - er_arr); - } else { - ComputeStressVarVisc_N(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, mu_turb, - cell_data, - s11, s22, s33, - s12, s13, s23, - er_arr); - } - - // Remove halo cells from tau_ii but extend across valid_box bdry - bxcc.grow(IntVect(-1,-1,0)); - if (bxcc.smallEnd(0) == valid_bx.smallEnd(0)) bxcc.growLo(0, 1); - if (bxcc.bigEnd(0) == valid_bx.bigEnd(0)) bxcc.growHi(0, 1); - if (bxcc.smallEnd(1) == valid_bx.smallEnd(1)) bxcc.growLo(1, 1); - if (bxcc.bigEnd(1) == valid_bx.bigEnd(1)) bxcc.growHi(1, 1); - - // Copy from temp FABs back to tau - ParallelFor(bxcc, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau11(i,j,k) = s11(i,j,k); - tau22(i,j,k) = s22(i,j,k); - tau33(i,j,k) = s33(i,j,k); - }); - ParallelFor(tbxxy, tbxxz, tbxyz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau12(i,j,k) = s12(i,j,k); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau13(i,j,k) = s13(i,j,k); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau23(i,j,k) = s23(i,j,k); - }); - } // end profile - } // l_use_terrain - } // MFIter } // l_use_diff - // HACK FOR PRINTING - // S_rhs[IntVars::cons].setVal(0.); - // ************************************************************************* // Define updates and fluxes in the current RK stage // ************************************************************************* @@ -494,6 +180,7 @@ void erf_slow_rhs_inc (int level, int nrk, #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif + { std::array flux; for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) @@ -544,6 +231,9 @@ void erf_slow_rhs_inc (int level, int nrk, const Array4& ymom_src_arr = ymom_src.const_array(mfi); const Array4& zmom_src_arr = zmom_src.const_array(mfi); + // Perturbational pressure + const Array4& pp_arr = pp_inc.const_array(mfi); + // Map factors const Array4& mf_m = mapfac_m->const_array(mfi); const Array4& mf_u = mapfac_u->const_array(mfi); @@ -562,7 +252,7 @@ void erf_slow_rhs_inc (int level, int nrk, const Array4& detJ_arr = detJ->const_array(mfi); // Base state - const Array4& p0_arr = p0->const_array(mfi); + // const Array4& p0_arr = p0->const_array(mfi); // ***************************************************************************** // ***************************************************************************** @@ -662,7 +352,7 @@ void erf_slow_rhs_inc (int level, int nrk, int n_start = amrex::max(start_comp,RhoTheta_comp); int n_comp = end_comp - n_start + 1; - DiffusionSrcForState_N(bx, domain, n_start, n_comp, u, v, + DiffusionSrcForState_N(bx, domain, n_start, n_comp, exp_most, u, v, cell_data, cell_prim, cell_rhs, diffflux_x, diffflux_y, diffflux_z, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, @@ -728,8 +418,9 @@ void erf_slow_rhs_inc (int level, int nrk, ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // x-momentum equation - // Note we do NOT include a pressure gradient here - rho_u_rhs(i, j, k) += xmom_src_arr(i,j,k) - solverChoice.abl_pressure_grad[0]; + Real gpx = dxInv[0] * (pp_arr(i,j,k) - pp_arr(i-1,j,k)); + + rho_u_rhs(i, j, k) += xmom_src_arr(i,j,k) - gpx - solverChoice.abl_pressure_grad[0]; if (nrk == 1) { rho_u_rhs(i,j,k) *= 0.5; @@ -740,8 +431,9 @@ void erf_slow_rhs_inc (int level, int nrk, ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // y-momentum equation - // Note we do NOT include a pressure gradient here - rho_v_rhs(i, j, k) += ymom_src_arr(i,j,k) - solverChoice.abl_pressure_grad[1]; + Real gpy = dxInv[1] * (pp_arr(i,j,k) - pp_arr(i,j-1,k)); + + rho_v_rhs(i, j, k) += ymom_src_arr(i,j,k) - gpy - solverChoice.abl_pressure_grad[1]; if (nrk == 1) { rho_v_rhs(i,j,k) *= 0.5; @@ -795,16 +487,15 @@ void erf_slow_rhs_inc (int level, int nrk, ParallelFor(tbz, [=] 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) ); + Real gpz = dxInv[2] * (pp_arr(i,j,k) - pp_arr(i,j,k-1)); - // Note we do NOT include a pressure gradient here - // HOWEVER: there may be a buoyancy term inside of zmom_src_arr ... - rho_w_rhs(i, j, k) += zmom_src_arr(i,j,k) - solverChoice.abl_pressure_grad[2]; + rho_w_rhs(i, j, k) += zmom_src_arr(i,j,k) - gpz - solverChoice.abl_pressure_grad[2]; - if (nrk == 1) { + 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)); - } + } }); } // mfi + } // omp } diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index b613b79cb..952d991a7 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -54,6 +54,7 @@ using namespace amrex; * @param[in] az area fractions on z-faces * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false) * @param[in] p0 Reference (hydrostatically stratified) pressure + * @param[in] pp_inc Perturbational pressure only used for incompressible flow * @param[in] mapfac_m map factor at cell centers * @param[in] mapfac_u map factor at x-faces * @param[in] mapfac_v map factor at y-faces @@ -94,6 +95,9 @@ void erf_slow_rhs_pre (int level, int finest_level, std::unique_ptr& az, std::unique_ptr& detJ, const MultiFab* p0, +#ifdef ERF_USE_POISSON_SOLVE + const MultiFab& pp_inc, +#endif std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v, @@ -171,328 +175,15 @@ void erf_slow_rhs_pre (int level, int finest_level, std::unique_ptr dflux_z; if (l_use_diff) { - expr = std::make_unique(ba , dm, 1, IntVect(1,1,0)); + erf_make_tau_terms(level,nrk,domain_bcs_type_d,domain_bcs_type_h,z_phys_nd, + S_data,xvel,yvel,zvel,Omega, + Tau11,Tau22,Tau33,Tau12,Tau13,Tau21,Tau23,Tau31,Tau32, + SmnSmn,eddyDiffs,Hfx3,Diss,geom,solverChoice,most, + detJ,mapfac_m,mapfac_u,mapfac_v); + dflux_x = std::make_unique(convert(ba,IntVect(1,0,0)), dm, nvars, 0); dflux_y = std::make_unique(convert(ba,IntVect(0,1,0)), dm, nvars, 0); dflux_z = std::make_unique(convert(ba,IntVect(0,0,1)), dm, nvars, 0); - - // if using constant alpha (mu = rho * alpha), then first divide by the - // reference density -- mu_eff will be scaled by the instantaneous - // local density later when ComputeStress*Visc_*() is called - Real mu_eff = (l_use_constAlpha) ? 2.0 * dc.dynamicViscosity / dc.rho0_trans - : 2.0 * dc.dynamicViscosity; - -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); - const Box& valid_bx = mfi.validbox(); - - // Velocities - const Array4 & u = xvel.array(mfi); - const Array4 & v = yvel.array(mfi); - const Array4 & w = zvel.array(mfi); - - // Contravariant velocity - const Array4& omega_arr = Omega.array(mfi); - - // Map factors - const Array4& mf_m = mapfac_m->const_array(mfi); - const Array4& mf_u = mapfac_u->const_array(mfi); - const Array4& mf_v = mapfac_v->const_array(mfi); - - // Eddy viscosity - const Array4& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : Array4{}; - const Array4& cell_data = l_use_constAlpha ? S_data[IntVars::cons].const_array(mfi) : Array4{}; - - // Terrain metrics - const Array4& z_nd = l_use_terrain ? z_phys_nd->const_array(mfi) : Array4{}; - const Array4& detJ_arr = detJ->const_array(mfi); - - //------------------------------------------------------------------------------- - // NOTE: Tile boxes with terrain are not intuitive. The linear combination of - // stress terms requires care. Create a tile box that intersects the - // valid box, then grow the box in x/y. Compute the strain on the local - // FAB over this grown tile box. Compute the stress over the tile box, - // except tau_ii which still needs the halo cells. Finally, write from - // the local FAB to the Tau MF but only on the tile box. - //------------------------------------------------------------------------------- - - //------------------------------------------------------------------------------- - // TODO: Avoid recomputing strain on the first RK stage. One could populate - // the FABs with tau_ij, compute stress, and then write to tau_ij. The - // problem with this approach is you will over-write the needed halo layer - // needed by subsequent tile boxes (particularly S_ii becomes Tau_ii). - //------------------------------------------------------------------------------- - - // Strain/Stress tile boxes - Box bxcc = mfi.tilebox(); - Box tbxxy = mfi.tilebox(IntVect(1,1,0)); - Box tbxxz = mfi.tilebox(IntVect(1,0,1)); - Box tbxyz = mfi.tilebox(IntVect(0,1,1)); - - // We need a halo cell for terrain - bxcc.grow(IntVect(1,1,0)); - tbxxy.grow(IntVect(1,1,0)); - tbxxz.grow(IntVect(1,1,0)); - tbxyz.grow(IntVect(1,1,0)); - - // Expansion rate - Array4 er_arr = expr->array(mfi); - - // Temporary storage for tiling/OMP - FArrayBox S11,S22,S33; - FArrayBox S12,S13,S23; - S11.resize( bxcc,1,The_Async_Arena()); S22.resize(bxcc,1,The_Async_Arena()); S33.resize(bxcc,1,The_Async_Arena()); - S12.resize(tbxxy,1,The_Async_Arena()); S13.resize(tbxxz,1,The_Async_Arena()); S23.resize(tbxyz,1,The_Async_Arena()); - Array4 s11 = S11.array(); Array4 s22 = S22.array(); Array4 s33 = S33.array(); - Array4 s12 = S12.array(); Array4 s13 = S13.array(); Array4 s23 = S23.array(); - - // Symmetric strain/stresses - Array4 tau11 = Tau11->array(mfi); Array4 tau22 = Tau22->array(mfi); Array4 tau33 = Tau33->array(mfi); - Array4 tau12 = Tau12->array(mfi); Array4 tau13 = Tau13->array(mfi); Array4 tau23 = Tau23->array(mfi); - - // Strain magnitude - Array4 SmnSmn_a; - - if (l_use_terrain) { - // Terrain non-symmetric terms - FArrayBox S21,S31,S32; - S21.resize(tbxxy,1,The_Async_Arena()); S31.resize(tbxxz,1,The_Async_Arena()); S32.resize(tbxyz,1,The_Async_Arena()); - Array4 s21 = S21.array(); Array4 s31 = S31.array(); Array4 s32 = S32.array(); - Array4 tau21 = Tau21->array(mfi); Array4 tau31 = Tau31->array(mfi); Array4 tau32 = Tau32->array(mfi); - - // ***************************************************************************** - // Expansion rate compute terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_er_T"); - // First create Omega using velocity (not momentum) - Box gbxo = surroundingNodes(bxcc,2); - ParallelFor(gbxo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - omega_arr(i,j,k) = (k == 0) ? 0. : OmegaFromW(i,j,k,w(i,j,k),u,v,z_nd,dxInv); - }); - - ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - - Real met_u_h_zeta_hi = Compute_h_zeta_AtIface(i+1, j , k, dxInv, z_nd); - Real met_u_h_zeta_lo = Compute_h_zeta_AtIface(i , j , k, dxInv, z_nd); - - Real met_v_h_zeta_hi = Compute_h_zeta_AtJface(i , j+1, k, dxInv, z_nd); - Real met_v_h_zeta_lo = Compute_h_zeta_AtJface(i , j , k, dxInv, z_nd); - - Real Omega_hi = omega_arr(i,j,k+1); - Real Omega_lo = omega_arr(i,j,k ); - - Real mfsq = mf_m(i,j,0)*mf_m(i,j,0); - - Real expansionRate = (u(i+1,j ,k)/mf_u(i+1,j,0)*met_u_h_zeta_hi - u(i,j,k)/mf_u(i,j,0)*met_u_h_zeta_lo)*dxInv[0]*mfsq + - (v(i ,j+1,k)/mf_v(i,j+1,0)*met_v_h_zeta_hi - v(i,j,k)/mf_v(i,j,0)*met_v_h_zeta_lo)*dxInv[1]*mfsq + - (Omega_hi - Omega_lo)*dxInv[2]; - - er_arr(i,j,k) = expansionRate / detJ_arr(i,j,k); - }); - } // end profile - - // ***************************************************************************** - // Strain tensor compute terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_strain_T"); - ComputeStrain_T(bxcc, tbxxy, tbxxz, tbxyz, domain, - u, v, w, - s11, s22, s33, - s12, s13, - s21, s23, - s31, s32, - z_nd, detJ_arr, bc_ptr_h, dxInv, - mf_m, mf_u, mf_v); - } // profile - - // Populate SmnSmn if using Deardorff (used as diff src in post) - // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF) - if ((nrk==0) && (tc.les_type == LESType::Deardorff)) { - SmnSmn_a = SmnSmn->array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most,exp_most); - }); - } - - // We've updated the strains at all locations including the - // surface. This is required to get the correct strain-rate - // magnitude. Now, update the stress everywhere but the surface - // to retain the values set by MOST. - if (use_most && exp_most) { - // Don't overwrite modeled total stress value at boundary - tbxxz.setSmall(2,1); - tbxyz.setSmall(2,1); - } - - // ***************************************************************************** - // Stress tensor compute terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_stress_T"); - - // Remove Halo cells just for tau_ij comps - tbxxy.grow(IntVect(-1,-1,0)); - tbxxz.grow(IntVect(-1,-1,0)); - tbxyz.grow(IntVect(-1,-1,0)); - - if (!l_use_turb) { - ComputeStressConsVisc_T(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, - cell_data, - s11, s22, s33, - s12, s13, - s21, s23, - s31, s32, - er_arr, z_nd, detJ_arr, dxInv); - } else { - ComputeStressVarVisc_T(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, mu_turb, - cell_data, - s11, s22, s33, - s12, s13, - s21, s23, - s31, s32, - er_arr, z_nd, detJ_arr, dxInv); - } - - // Remove halo cells from tau_ii but extend across valid_box bdry - bxcc.grow(IntVect(-1,-1,0)); - if (bxcc.smallEnd(0) == valid_bx.smallEnd(0)) bxcc.growLo(0, 1); - if (bxcc.bigEnd(0) == valid_bx.bigEnd(0)) bxcc.growHi(0, 1); - if (bxcc.smallEnd(1) == valid_bx.smallEnd(1)) bxcc.growLo(1, 1); - if (bxcc.bigEnd(1) == valid_bx.bigEnd(1)) bxcc.growHi(1, 1); - - // Copy from temp FABs back to tau - ParallelFor(bxcc, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau11(i,j,k) = s11(i,j,k); - tau22(i,j,k) = s22(i,j,k); - tau33(i,j,k) = s33(i,j,k); - }); - - ParallelFor(tbxxy, tbxxz, tbxyz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau12(i,j,k) = s12(i,j,k); - tau21(i,j,k) = s21(i,j,k); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau13(i,j,k) = s13(i,j,k); - tau31(i,j,k) = s31(i,j,k); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau23(i,j,k) = s23(i,j,k); - tau32(i,j,k) = s32(i,j,k); - }); - } // end profile - - } else { - - // ***************************************************************************** - // Expansion rate compute no terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_er_N"); - ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real mfsq = mf_m(i,j,0)*mf_m(i,j,0); - er_arr(i,j,k) = (u(i+1, j , k )/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0]*mfsq + - (v(i , j+1, k )/mf_v(i,j+1,0) - v(i, j, k)/mf_v(i,j,0))*dxInv[1]*mfsq + - (w(i , j , k+1) - w(i, j, k))*dxInv[2]; - }); - } // end profile - - - // ***************************************************************************** - // Strain tensor compute no terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_strain_N"); - ComputeStrain_N(bxcc, tbxxy, tbxxz, tbxyz, domain, - u, v, w, - s11, s22, s33, - s12, s13, s23, - bc_ptr_h, dxInv, - mf_m, mf_u, mf_v); - } // end profile - - // Populate SmnSmn if using Deardorff (used as diff src in post) - // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF) - if ((nrk==0) && (tc.les_type == LESType::Deardorff)) { - SmnSmn_a = SmnSmn->array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most,exp_most); - }); - } - - // We've updated the strains at all locations including the - // surface. This is required to get the correct strain-rate - // magnitude. Now, update the stress everywhere but the surface - // to retain the values set by MOST. - if (use_most && exp_most) { - // Don't overwrite modeled total stress value at boundary - tbxxz.setSmall(2,1); - tbxyz.setSmall(2,1); - } - - // ***************************************************************************** - // Stress tensor compute no terrain - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_stress_N"); - - // Remove Halo cells just for tau_ij comps - tbxxy.grow(IntVect(-1,-1,0)); - tbxxz.grow(IntVect(-1,-1,0)); - tbxyz.grow(IntVect(-1,-1,0)); - - if (!l_use_turb) { - ComputeStressConsVisc_N(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, - cell_data, - s11, s22, s33, - s12, s13, s23, - er_arr); - } else { - ComputeStressVarVisc_N(bxcc, tbxxy, tbxxz, tbxyz, mu_eff, mu_turb, - cell_data, - s11, s22, s33, - s12, s13, s23, - er_arr); - } - - // Remove halo cells from tau_ii but extend across valid_box bdry - bxcc.grow(IntVect(-1,-1,0)); - if (bxcc.smallEnd(0) == valid_bx.smallEnd(0)) bxcc.growLo(0, 1); - if (bxcc.bigEnd(0) == valid_bx.bigEnd(0)) bxcc.growHi(0, 1); - if (bxcc.smallEnd(1) == valid_bx.smallEnd(1)) bxcc.growLo(1, 1); - if (bxcc.bigEnd(1) == valid_bx.bigEnd(1)) bxcc.growHi(1, 1); - - // Copy from temp FABs back to tau - ParallelFor(bxcc, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau11(i,j,k) = s11(i,j,k); - tau22(i,j,k) = s22(i,j,k); - tau33(i,j,k) = s33(i,j,k); - }); - ParallelFor(tbxxy, tbxxz, tbxyz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau12(i,j,k) = s12(i,j,k); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau13(i,j,k) = s13(i,j,k); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau23(i,j,k) = s23(i,j,k); - }); - } // end profile - } // l_use_terrain - } // MFIter } // l_use_diff // ***************************************************************************** diff --git a/Source/TimeIntegration/Make.package b/Source/TimeIntegration/Make.package index 2641d4633..8a845a8cc 100644 --- a/Source/TimeIntegration/Make.package +++ b/Source/TimeIntegration/Make.package @@ -5,6 +5,7 @@ CEXE_sources += ERF_advance_dycore.cpp CEXE_sources += ERF_advance_microphysics.cpp CEXE_sources += ERF_advance_lsm.cpp CEXE_sources += ERF_make_fast_coeffs.cpp +CEXE_sources += ERF_make_tau_terms.cpp CEXE_sources += ERF_slow_rhs_pre.cpp CEXE_sources += ERF_slow_rhs_post.cpp CEXE_sources += ERF_fast_rhs_N.cpp diff --git a/Source/TimeIntegration/TI_no_substep_fun.H b/Source/TimeIntegration/TI_no_substep_fun.H index 5599db59b..bebba6845 100644 --- a/Source/TimeIntegration/TI_no_substep_fun.H +++ b/Source/TimeIntegration/TI_no_substep_fun.H @@ -139,9 +139,9 @@ bool have_tb = (thin_xforce[0] || thin_yforce[0] || thin_zforce[0]); if (solverChoice.project_every_stage || (nrk==1)) { if (have_tb) { - project_velocities_tb(S_sum, slow_dt); + project_velocities_tb(level, slow_dt, S_sum, pp_inc[level]); } else { - project_velocities(S_sum); + project_velocities(level, slow_dt, S_sum, pp_inc[level]); } } } diff --git a/Source/TimeIntegration/TI_slow_headers.H b/Source/TimeIntegration/TI_slow_headers.H index 7d8621dbf..780d74245 100644 --- a/Source/TimeIntegration/TI_slow_headers.H +++ b/Source/TimeIntegration/TI_slow_headers.H @@ -20,6 +20,36 @@ #include #endif +void erf_make_tau_terms (int level, int nrk, + const amrex::Gpu::DeviceVector& domain_bcs_type_d, + const amrex::Vector& domain_bcs_type, + std::unique_ptr& z_phys_nd, + amrex::Vector& S_data, + const amrex::MultiFab& xvel, + const amrex::MultiFab& yvel, + const amrex::MultiFab& zvel, + amrex::MultiFab& Omega, + amrex::MultiFab* Tau11, + amrex::MultiFab* Tau22, + amrex::MultiFab* Tau33, + amrex::MultiFab* Tau12, + amrex::MultiFab* Tau13, + amrex::MultiFab* Tau21, + amrex::MultiFab* Tau23, + amrex::MultiFab* Tau31, + amrex::MultiFab* Tau32, + amrex::MultiFab* SmnSmn, + amrex::MultiFab* eddyDiffs, + amrex::MultiFab* Hfx3, + amrex::MultiFab* Diss, + const amrex::Geometry geom, + const SolverChoice& solverChoice, + std::unique_ptr& most, + std::unique_ptr& dJ, + std::unique_ptr& mapfac_m, + std::unique_ptr& mapfac_u, + std::unique_ptr& mapfac_v); + /** * Function for computing the slow RHS for the evolution equations for the density, potential temperature and momentum. * @@ -63,6 +93,9 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk, std::unique_ptr& az, std::unique_ptr& dJ, const amrex::MultiFab* p0, +#ifdef ERF_USE_POISSON_SOLVE + const amrex::MultiFab& pp_inc, +#endif std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v, @@ -171,6 +204,7 @@ void erf_slow_rhs_inc (int level, int nrk, const amrex::MultiFab* p0, std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, - std::unique_ptr& mapfac_v); + std::unique_ptr& mapfac_v, + const amrex::MultiFab& pp_inc); #endif #endif diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 6868ff216..41c2a5665 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -127,6 +127,9 @@ Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx3, Diss, fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, z_phys_nd_src[level], ax_src[level], ay_src[level], az_src[level], detJ_cc_src[level], p0_new, +#ifdef ERF_USE_POISSON_SOLVE + pp_inc[level], +#endif mapfac_m[level], mapfac_u[level], mapfac_v[level], #ifdef ERF_USE_EB EBFactory(level), @@ -227,6 +230,9 @@ Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx3, Diss, fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0, +#ifdef ERF_USE_POISSON_SOLVE + pp_inc[level], +#endif mapfac_m[level], mapfac_u[level], mapfac_v[level], #ifdef ERF_USE_EB EBFactory(level), @@ -416,7 +422,7 @@ Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx3, Diss, fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0, - mapfac_m[level], mapfac_u[level], mapfac_v[level]); + mapfac_m[level], mapfac_u[level], mapfac_v[level], pp_inc[level]); add_thin_body_sources(xmom_src, ymom_src, zmom_src, xflux_imask[level], yflux_imask[level], zflux_imask[level], diff --git a/Source/Utils/ERF_PoissonSolve.cpp b/Source/Utils/ERF_PoissonSolve.cpp index dd9314670..80d9d5a4f 100644 --- a/Source/Utils/ERF_PoissonSolve.cpp +++ b/Source/Utils/ERF_PoissonSolve.cpp @@ -1,6 +1,6 @@ #include "ERF.H" #include -#include +#include #include "Utils.H" #ifdef ERF_USE_POISSON_SOLVE @@ -12,6 +12,8 @@ using namespace amrex; * if we want to enforce incompressibility of the initial conditions */ +using BCType = LinOpBCType; + Array ERF::get_projection_bc (Orientation::Side side) const noexcept { @@ -44,92 +46,105 @@ bool ERF::projection_has_dirichlet (Array bcs) const /** * Project the single-level velocity field to enforce incompressibility - * Note this assumes that we are projecting the level 0 field using domain bcs + * Note that the level may or may not be level 0. */ -void ERF::project_velocities (Vector& vmf) -{ - Vector> tmpmf(1); - for (auto& mf : vmf) { - tmpmf[0].emplace_back(mf, amrex::make_alias, 0, mf.nComp()); - } - project_velocities(0,0,tmpmf); -} - -/** - * Project the multi-level velocity field to enforce incompressibility - */ -void -ERF::project_velocities (int lev_min, int lev_max, Vector>& vars) +void ERF::project_velocities (int lev, Real l_dt, Vector& vmf, MultiFab& pmf) { BL_PROFILE("ERF::project_velocities()"); + AMREX_ALWAYS_ASSERT(!solverChoice.use_terrain); - const int nlevs = lev_max - lev_min + 1; - - // Use the default settings + // Make sure the solver only sees the levels over which we are solving LPInfo info; - MLPoisson mlpoisson(geom, grids, dmap, info); + Vector ba_tmp; ba_tmp.push_back(vmf[Vars::cons].boxArray()); + Vector dm_tmp; dm_tmp.push_back(vmf[Vars::cons].DistributionMap()); + Vector geom_tmp; geom_tmp.push_back(geom[lev]); + // amrex::Print() << "AT LEVEL " << lev << " BA FOR POISSON SOLVE " << vmf[Vars::cons].boxArray() << std::endl; + + MLABecLaplacian mlabec(geom_tmp, ba_tmp, dm_tmp, info); + + // + // This will hold (1/rho) on faces + // + Array inv_rho; + + // + // The operator is (alpha A - beta del dot B grad) phi = RHS + // Here we set alpha to 0 and beta to -1 + // Then b is (dt/rho) + // + mlabec.setScalars(0.0, -1.0); + inv_rho[0].define(vmf[Vars::xvel].boxArray(),dm_tmp[0],1,0,MFInfo()); + inv_rho[1].define(vmf[Vars::yvel].boxArray(),dm_tmp[0],1,0,MFInfo()); + inv_rho[2].define(vmf[Vars::zvel].boxArray(),dm_tmp[0],1,0,MFInfo()); + +#if 0 + MultiFab density(vmf[Vars::cons], make_alias, 1, 1); + density.FillBoundary(geom_tmp[0].periodicity()); + amrex::average_cellcenter_to_face(GetArrOfPtrs(inv_rho), density, geom_tmp[0]); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + inv_rho[idim].invert(l_dt, 0); + } +#else + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + inv_rho[idim].setVal(l_dt); + } +#endif + + mlabec.setBCoeffs(0, GetArrOfConstPtrs(inv_rho)); auto bclo = get_projection_bc(Orientation::low); auto bchi = get_projection_bc(Orientation::high); bool need_adjust_rhs = (projection_has_dirichlet(bclo) || projection_has_dirichlet(bchi)) ? false : true; - mlpoisson.setDomainBC(bclo, bchi); + mlabec.setDomainBC(bclo, bchi); - for (int ilev = lev_min; ilev <= lev_max; ++ilev) { - mlpoisson.setLevelBC(0, nullptr); + if (lev > 0) { + mlabec.setCoarseFineBC(nullptr, ref_ratio[lev-1], LinOpBCType::Neumann); } + mlabec.setLevelBC(0, nullptr); Vector rhs; Vector phi; Vector > fluxes; - rhs.resize(nlevs); - phi.resize(nlevs); - fluxes.resize(nlevs); + rhs.resize(1); + phi.resize(1); + fluxes.resize(1); - for (int ilev = lev_min; ilev <= lev_max; ++ilev) { - rhs[ilev-lev_min].define(grids[ilev], dmap[ilev], 1, 0); - phi[ilev-lev_min].define(grids[ilev], dmap[ilev], 1, 1); - rhs[ilev-lev_min].setVal(0.0); - phi[ilev-lev_min].setVal(0.0); + rhs[0].define(ba_tmp[0], dm_tmp[0], 1, 0); + phi[0].define(ba_tmp[0], dm_tmp[0], 1, 0); + rhs[0].setVal(0.0); + phi[0].setVal(0.0); - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - fluxes[ilev-lev_min][idim].define( - convert(grids[ilev], IntVect::TheDimensionVector(idim)), - dmap[ilev], 1, 0); - } + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + fluxes[0][idim].define(convert(ba_tmp[0], IntVect::TheDimensionVector(idim)), dm_tmp[0], 1, 0); } // Define a single Array of MultiFabs to hold the velocity components // Then define the RHS to be the divergence of the current velocities Array u; - for (int ilev = lev_min; ilev <= lev_max; ++ilev) + u[0] = &(vmf[Vars::xvel]); + u[1] = &(vmf[Vars::yvel]); + u[2] = &(vmf[Vars::zvel]); + computeDivergence(rhs[0], u, geom_tmp[0]); + + // If all Neumann BCs, adjust RHS to make sure we can converge + if (need_adjust_rhs) { - u[0] = &(vars[ilev][Vars::xvel]); - u[1] = &(vars[ilev][Vars::yvel]); - u[2] = &(vars[ilev][Vars::zvel]); - computeDivergence(rhs[ilev-lev_min], u, geom[ilev]); - - // If all Neumann BCs, adjust RHS to make sure we can converge - if (need_adjust_rhs) { - Real offset = volWgtSumMF(ilev, rhs[ilev], 0, *mapfac_m[ilev], false, false); - amrex::Print() << "Poisson solvability offset = " << offset << std::endl; - rhs[ilev-lev_min].plus(-offset, 0, 1); - } + Real offset = volWgtSumMF(lev, rhs[0], 0, *mapfac_m[lev], false, false); + amrex::Print() << "Poisson solvability offset = " << offset << std::endl; + rhs[0].plus(-offset, 0, 1); } // Initialize phi to 0 - for (int ilev = lev_min; ilev <= lev_max; ++ilev) - { - phi[ilev-lev_min].setVal(0.0); - } + phi[0].setVal(0.0); - MLMG mlmg(mlpoisson); + MLMG mlmg(mlabec); int max_iter = 100; mlmg.setMaxIter(max_iter); - int verbose = 1; - mlmg.setVerbose(verbose); - //mlmg.setBottomVerbose(verbose); + mlmg.setVerbose(mg_verbose); + mlmg.setBottomVerbose(0); mlmg.solve(GetVecOfPtrs(phi), GetVecOfConstPtrs(rhs), @@ -138,39 +153,23 @@ ERF::project_velocities (int lev_min, int lev_max, Vector>& var mlmg.getFluxes(GetVecOfArrOfPtrs(fluxes)); - // Subtract grad(phi) from the velocity components - Real beta = 1.0; - for (int ilev = lev_min; ilev <= lev_max; ++ilev) { - MultiFab::Saxpy(vars[ilev][Vars::xvel], beta, fluxes[ilev-lev_min][0], 0,0,1,0); - MultiFab::Saxpy(vars[ilev][Vars::yvel], beta, fluxes[ilev-lev_min][1], 0,0,1,0); - MultiFab::Saxpy(vars[ilev][Vars::zvel], beta, fluxes[ilev-lev_min][2], 0,0,1,0); - } + // Update pressure variable with phi -- note that phi is change in pressure, not the full pressure + MultiFab::Saxpy(pmf, 1.0, phi[0],0,0,1,0); + pmf.FillBoundary(geom[lev].periodicity()); - // Average down the velocity from finest to coarsest to ensure consistency across levels - Array u_fine; - Array u_crse; - for (int ilev = lev_min+1; ilev <= lev_max; ++ilev) - { - IntVect rr = geom[ilev].Domain().size() / geom[ilev-1].Domain().size(); - u_fine[0] = &(vars[ilev ][Vars::xvel]); - u_fine[1] = &(vars[ilev ][Vars::yvel]); - u_fine[2] = &(vars[ilev ][Vars::zvel]); - u_crse[0] = &(vars[ilev-1][Vars::xvel]); - u_crse[1] = &(vars[ilev-1][Vars::yvel]); - u_crse[2] = &(vars[ilev-1][Vars::zvel]); - average_down_faces(u_fine, u_crse, rr, geom[ilev-1]); - } + // Subtract (dt/rho) grad(phi) from the velocity components + MultiFab::Add(vmf[Vars::xvel], fluxes[0][0], 0,0,1,0); + MultiFab::Add(vmf[Vars::yvel], fluxes[0][1], 0,0,1,0); + MultiFab::Add(vmf[Vars::zvel], fluxes[0][2], 0,0,1,0); #if 1 // Confirm that the velocity is now divergence free - for (int ilev = lev_min; ilev <= lev_max; ++ilev) - { - u[0] = &(vars[ilev][Vars::xvel]); - u[1] = &(vars[ilev][Vars::yvel]); - u[2] = &(vars[ilev][Vars::zvel]); - computeDivergence(rhs[ilev], u, geom[ilev]); - Print() << "Max norm of divergence after solve at level " << ilev << " : " << rhs[ilev].norm0() << std::endl; - } + u[0] = &(vmf[Vars::xvel]); + u[1] = &(vmf[Vars::yvel]); + u[2] = &(vmf[Vars::zvel]); + rhs[0].setVal(0.0); + computeDivergence(rhs[0], u, geom_tmp[0]); + Print() << "Max norm of divergence after solve at level " << lev << " : " << rhs[0].norm0() << std::endl; #endif } #endif diff --git a/Source/Utils/ERF_PoissonSolve_tb.cpp b/Source/Utils/ERF_PoissonSolve_tb.cpp index 1d364163e..cdeb207fc 100644 --- a/Source/Utils/ERF_PoissonSolve_tb.cpp +++ b/Source/Utils/ERF_PoissonSolve_tb.cpp @@ -11,26 +11,16 @@ using namespace amrex; * Project the single-level velocity field to enforce incompressibility with a * thin body */ -void ERF::project_velocities_tb (Vector& vmf, const Real dt) -{ - Vector> tmpmf(1); - for (auto& mf : vmf) { - tmpmf[0].emplace_back(mf, amrex::make_alias, 0, mf.nComp()); - } - project_velocities_tb(tmpmf, dt); -} - -/** - * Project the multi-level velocity field to enforce incompressibility with a - * thin body; we iterate on the pressure solve to obtain the body force that - * satisfies no penetration - */ -void -ERF::project_velocities_tb (Vector>& vars, const Real dt) +void ERF::project_velocities_tb (int lev, Real l_dt, Vector& vmf, MultiFab& pmf) { BL_PROFILE("ERF::project_velocities()"); + AMREX_ALWAYS_ASSERT(!solverChoice.use_terrain); - const int nlevs = geom.size(); + // Make sure the solver only sees the levels over which we are solving + Vector ba_tmp; ba_tmp.push_back(vmf[Vars::cons].boxArray()); + Vector dm_tmp; dm_tmp.push_back(vmf[Vars::cons].DistributionMap()); + Vector geom_tmp; geom_tmp.push_back(geom[lev]); + amrex::Print() << "AT LEVEL " << lev << " BA FOR POISSON SOLVE " << vmf[Vars::cons].boxArray() << std::endl; // Use the default settings LPInfo info; @@ -44,7 +34,7 @@ ERF::project_velocities_tb (Vector>& vars, const Real dt) #endif { // Use the default settings - p_mlpoisson = std::make_unique(geom, grids, dmap, info); + p_mlpoisson = std::make_unique(geom_tmp, ba_tmp, dm_tmp, info); } auto bclo = get_projection_bc(Orientation::low); @@ -52,10 +42,12 @@ ERF::project_velocities_tb (Vector>& vars, const Real dt) bool need_adjust_rhs = (projection_has_dirichlet(bclo) || projection_has_dirichlet(bchi)) ? false : true; p_mlpoisson->setDomainBC(bclo, bchi); - for (int ilev = 0; ilev < nlevs; ++ilev) { - p_mlpoisson->setLevelBC(0, nullptr); + if (lev > 0) { + p_mlpoisson->setCoarseFineBC(nullptr, ref_ratio[lev-1], LinOpBCType::Neumann); } + p_mlpoisson->setLevelBC(0, nullptr); + Vector rhs; Vector phi; Vector > fluxes; @@ -65,100 +57,86 @@ ERF::project_velocities_tb (Vector>& vars, const Real dt) // Used to pass array of const MFs to ComputeDivergence Array u; - rhs.resize(nlevs); - phi.resize(nlevs); - fluxes.resize(nlevs); - deltaf.resize(nlevs); - u_plus_dtdf.resize(nlevs); + rhs.resize(1); + phi.resize(1); + fluxes.resize(1); + deltaf.resize(1); + u_plus_dtdf.resize(1); - for (int ilev = 0; ilev < nlevs; ++ilev) { - rhs[ilev].define(grids[ilev], dmap[ilev], 1, 0); - phi[ilev].define(grids[ilev], dmap[ilev], 1, 1); - rhs[ilev].setVal(0.0); - phi[ilev].setVal(0.0); + rhs[0].define(ba_tmp[0], dm_tmp[0], 1, 0); + phi[0].define(ba_tmp[0], dm_tmp[0], 1, 0); + rhs[0].setVal(0.0); + phi[0].setVal(0.0); - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - fluxes[ilev][idim].define( - convert(grids[ilev], IntVect::TheDimensionVector(idim)), - dmap[ilev], 1, 0); - u_plus_dtdf[ilev][idim].define( - convert(grids[ilev], IntVect::TheDimensionVector(idim)), - dmap[ilev], 1, 0); - deltaf[ilev][idim].define( - convert(grids[ilev], IntVect::TheDimensionVector(idim)), - dmap[ilev], 1, 0); - deltaf[ilev][idim].setVal(0.0); // start with f^* == f^{n-1} - } + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + fluxes[0][idim].define(convert(ba_tmp[0], IntVect::TheDimensionVector(idim)), dm_tmp[0], 1, 0); + u_plus_dtdf[0][idim].define(convert(ba_tmp[0], IntVect::TheDimensionVector(idim)), dm_tmp[0], 1, 0); + + deltaf[0][idim].define(convert(ba_tmp[0], IntVect::TheDimensionVector(idim)), dm_tmp[0], 1, 0); + deltaf[0][idim].setVal(0.0); // start with f^* == f^{n-1} } // DEBUG - u[0] = &(vars[0][Vars::xvel]); - u[1] = &(vars[0][Vars::yvel]); - u[2] = &(vars[0][Vars::zvel]); + u[0] = &(vmf[Vars::xvel]); + u[1] = &(vmf[Vars::yvel]); + u[2] = &(vmf[Vars::zvel]); computeDivergence(rhs[0], u, geom[0]); Print() << "Max norm of divergence before solve at level 0 : " << rhs[0].norm0() << std::endl; for (int itp = 0; itp < solverChoice.ncorr; ++itp) { - for (int ilev = 0; ilev < nlevs; ++ilev) - { - // Calculate u + dt*deltaf - for (int idim = 0; idim < 3; ++idim) { - MultiFab::Copy(u_plus_dtdf[ilev][idim], deltaf[ilev][idim], 0, 0, 1, 0); - u_plus_dtdf[ilev][0].mult(-dt,0,1,0); - } - MultiFab::Add(u_plus_dtdf[ilev][0], vars[ilev][Vars::xvel], 0, 0, 1, 0); - MultiFab::Add(u_plus_dtdf[ilev][1], vars[ilev][Vars::yvel], 0, 0, 1, 0); - MultiFab::Add(u_plus_dtdf[ilev][2], vars[ilev][Vars::zvel], 0, 0, 1, 0); - - u[0] = &(u_plus_dtdf[ilev][0]); - u[1] = &(u_plus_dtdf[ilev][1]); - u[2] = &(u_plus_dtdf[ilev][2]); - computeDivergence(rhs[ilev], u, geom[ilev]); - - // DEBUG - if (itp==0) { - for (MFIter mfi(rhs[ilev], TilingIfNotGPU()); mfi.isValid(); ++mfi) + // Calculate u + dt*deltaf + for (int idim = 0; idim < 3; ++idim) { + MultiFab::Copy(u_plus_dtdf[0][idim], deltaf[0][idim], 0, 0, 1, 0); + u_plus_dtdf[0][0].mult(-l_dt,0,1,0); + } + MultiFab::Add(u_plus_dtdf[0][0], vmf[Vars::xvel], 0, 0, 1, 0); + MultiFab::Add(u_plus_dtdf[0][1], vmf[Vars::yvel], 0, 0, 1, 0); + MultiFab::Add(u_plus_dtdf[0][2], vmf[Vars::zvel], 0, 0, 1, 0); + + u[0] = &(u_plus_dtdf[0][0]); + u[1] = &(u_plus_dtdf[0][1]); + u[2] = &(u_plus_dtdf[0][2]); + computeDivergence(rhs[0], u, geom_tmp[0]); + + // DEBUG + if (itp==0) { + for (MFIter mfi(rhs[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& divU = rhs[0].const_array(mfi); + const Array4& uarr = vmf[Vars::xvel].const_array(mfi); + const Array4& varr = vmf[Vars::yvel].const_array(mfi); + const Array4& warr = vmf[Vars::zvel].const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const Box& bx = mfi.tilebox(); - const Array4& divU = rhs[ilev].const_array(mfi); - const Array4& uarr = vars[ilev][Vars::xvel].const_array(mfi); - const Array4& varr = vars[ilev][Vars::yvel].const_array(mfi); - const Array4& warr = vars[ilev][Vars::zvel].const_array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - if ((i>=120) && (i<=139) && (j==0) && ((k>=127)&&(k<=128))) { - amrex::AllPrint() << "before project div"<=120) && (i<=139) && (j==0) && ((k>=127)&&(k<=128))) { + amrex::AllPrint() << "before project div"<& divU = rhs[0].const_array(mfi); + const Array4& uarr = u[0]->const_array(mfi); + const Array4& varr = u[1]->const_array(mfi); + const Array4& warr = u[2]->const_array(mfi); + const Array4& fzarr = thin_zforce[0]->const_array(mfi); + const Array4& dfzarr = deltaf[0][2].const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const Box& bx = mfi.tilebox(); - const Array4& divU = rhs[ilev].const_array(mfi); - const Array4& uarr = u[0]->const_array(mfi); - const Array4& varr = u[1]->const_array(mfi); - const Array4& warr = u[2]->const_array(mfi); - const Array4& fzarr = thin_zforce[ilev]->const_array(mfi); - const Array4& dfzarr = deltaf[ilev][2].const_array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - if ((i>=120) && (i<=139) && (j==0) && ((k>=127)&&(k<=128))) { - amrex::AllPrint() << "after project div"<=120) && (i<=139) && (j==0) && ((k>=127)&&(k<=128))) { + amrex::AllPrint() << "after project div"< Date: Thu, 9 May 2024 18:53:14 -0700 Subject: [PATCH 03/20] fix compilation with POISSON_SOLVE (#1613) --- Source/TimeIntegration/ERF_slow_rhs_inc.cpp | 1 - Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 1 - Source/TimeIntegration/TI_slow_headers.H | 4 ++-- Source/TimeIntegration/TI_slow_rhs_fun.H | 4 ++-- 4 files changed, 4 insertions(+), 6 deletions(-) diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index 45a096757..076eef8d9 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -122,7 +122,6 @@ void erf_slow_rhs_inc (int level, int nrk, const bool l_use_diff = ( (dc.molec_diff_type != MolecDiffType::None) || (tc.les_type != LESType::None) || (tc.pbl_type != PBLType::None) ); - const bool l_use_constAlpha = ( dc.molec_diff_type == MolecDiffType::ConstantAlpha ); const bool l_use_turb = ( tc.les_type == LESType::Smagorinsky || tc.les_type == LESType::Deardorff || tc.pbl_type == PBLType::MYNN25 || diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 952d991a7..9fb3a964c 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -139,7 +139,6 @@ void erf_slow_rhs_pre (int level, int finest_level, const bool l_use_diff = ( (dc.molec_diff_type != MolecDiffType::None) || (tc.les_type != LESType::None) || (tc.pbl_type != PBLType::None) ); - const bool l_use_constAlpha = ( dc.molec_diff_type == MolecDiffType::ConstantAlpha ); const bool l_use_turb = ( tc.les_type == LESType::Smagorinsky || tc.les_type == LESType::Deardorff || tc.pbl_type == PBLType::MYNN25 || diff --git a/Source/TimeIntegration/TI_slow_headers.H b/Source/TimeIntegration/TI_slow_headers.H index 780d74245..c3d350429 100644 --- a/Source/TimeIntegration/TI_slow_headers.H +++ b/Source/TimeIntegration/TI_slow_headers.H @@ -202,9 +202,9 @@ void erf_slow_rhs_inc (int level, int nrk, std::unique_ptr& az, std::unique_ptr& dJ, const amrex::MultiFab* p0, + const amrex::MultiFab& pp_inc, std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, - std::unique_ptr& mapfac_v, - const amrex::MultiFab& pp_inc); + std::unique_ptr& mapfac_v); #endif #endif diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 41c2a5665..23694da64 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -421,8 +421,8 @@ Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(), Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx3, Diss, fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, - z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0, - mapfac_m[level], mapfac_u[level], mapfac_v[level], pp_inc[level]); + z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0, pp_inc[level], + mapfac_m[level], mapfac_u[level], mapfac_v[level]); add_thin_body_sources(xmom_src, ymom_src, zmom_src, xflux_imask[level], yflux_imask[level], zflux_imask[level], From 871395c48c7c5291603cad2358b0972626c2ed20 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 10 May 2024 08:51:08 -0700 Subject: [PATCH 04/20] Remove prints (#1614) * remove print statements * remove print statements --- Source/Utils/ERF_PoissonSolve.cpp | 4 ++-- Source/Utils/ERF_PoissonSolve_tb.cpp | 12 +++++++++--- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/Source/Utils/ERF_PoissonSolve.cpp b/Source/Utils/ERF_PoissonSolve.cpp index 80d9d5a4f..eeb2b5a50 100644 --- a/Source/Utils/ERF_PoissonSolve.cpp +++ b/Source/Utils/ERF_PoissonSolve.cpp @@ -132,7 +132,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& vmf, MultiFa if (need_adjust_rhs) { Real offset = volWgtSumMF(lev, rhs[0], 0, *mapfac_m[lev], false, false); - amrex::Print() << "Poisson solvability offset = " << offset << std::endl; + // amrex::Print() << "Poisson solvability offset = " << offset << std::endl; rhs[0].plus(-offset, 0, 1); } @@ -162,7 +162,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& vmf, MultiFa MultiFab::Add(vmf[Vars::yvel], fluxes[0][1], 0,0,1,0); MultiFab::Add(vmf[Vars::zvel], fluxes[0][2], 0,0,1,0); -#if 1 +#if 0 // Confirm that the velocity is now divergence free u[0] = &(vmf[Vars::xvel]); u[1] = &(vmf[Vars::yvel]); diff --git a/Source/Utils/ERF_PoissonSolve_tb.cpp b/Source/Utils/ERF_PoissonSolve_tb.cpp index cdeb207fc..2daa04013 100644 --- a/Source/Utils/ERF_PoissonSolve_tb.cpp +++ b/Source/Utils/ERF_PoissonSolve_tb.cpp @@ -20,7 +20,7 @@ void ERF::project_velocities_tb (int lev, Real l_dt, Vector& vmf, Mult Vector ba_tmp; ba_tmp.push_back(vmf[Vars::cons].boxArray()); Vector dm_tmp; dm_tmp.push_back(vmf[Vars::cons].DistributionMap()); Vector geom_tmp; geom_tmp.push_back(geom[lev]); - amrex::Print() << "AT LEVEL " << lev << " BA FOR POISSON SOLVE " << vmf[Vars::cons].boxArray() << std::endl; + // amrex::Print() << "AT LEVEL " << lev << " BA FOR POISSON SOLVE " << vmf[Vars::cons].boxArray() << std::endl; // Use the default settings LPInfo info; @@ -76,12 +76,14 @@ void ERF::project_velocities_tb (int lev, Real l_dt, Vector& vmf, Mult deltaf[0][idim].setVal(0.0); // start with f^* == f^{n-1} } +#if 0 // DEBUG u[0] = &(vmf[Vars::xvel]); u[1] = &(vmf[Vars::yvel]); u[2] = &(vmf[Vars::zvel]); computeDivergence(rhs[0], u, geom[0]); Print() << "Max norm of divergence before solve at level 0 : " << rhs[0].norm0() << std::endl; +#endif for (int itp = 0; itp < solverChoice.ncorr; ++itp) { @@ -99,6 +101,7 @@ void ERF::project_velocities_tb (int lev, Real l_dt, Vector& vmf, Mult u[2] = &(u_plus_dtdf[0][2]); computeDivergence(rhs[0], u, geom_tmp[0]); +#if 0 // DEBUG if (itp==0) { for (MFIter mfi(rhs[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) @@ -120,11 +123,12 @@ void ERF::project_velocities_tb (int lev, Real l_dt, Vector& vmf, Mult }); } } +#endif // If all Neumann BCs, adjust RHS to make sure we can converge if (need_adjust_rhs) { Real offset = volWgtSumMF(lev, rhs[0], 0, *mapfac_m[lev], false, false); - amrex::Print() << "Poisson solvability offset = " << offset << std::endl; + // amrex::Print() << "Poisson solvability offset = " << offset << std::endl; rhs[0].plus(-offset, 0, 1); } @@ -212,7 +216,7 @@ void ERF::project_velocities_tb (int lev, Real l_dt, Vector& vmf, Mult // } // } -#if 1 +#if 0 // Confirm that the velocity is now divergence free u[0] = &(vmf[Vars::xvel]); u[1] = &(vmf[Vars::yvel]); @@ -220,6 +224,7 @@ void ERF::project_velocities_tb (int lev, Real l_dt, Vector& vmf, Mult computeDivergence(rhs[0], u, geom_tmp[0]); Print() << "Max norm of divergence after solve at level " << lev << " : " << rhs[0].norm0() << std::endl; +#if 0 for (MFIter mfi(rhs[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); @@ -242,5 +247,6 @@ void ERF::project_velocities_tb (int lev, Real l_dt, Vector& vmf, Mult }); } #endif +#endif } #endif From 9349d664f9e8bfdd55984653ee53281ab1210d0a Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Fri, 10 May 2024 13:40:13 -0700 Subject: [PATCH 05/20] Moisture average prof (#1616) * Wrong indexing for different moisture models. * missed one other case. --- Source/IO/ERF_Write1DProfiles.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Source/IO/ERF_Write1DProfiles.cpp b/Source/IO/ERF_Write1DProfiles.cpp index 2c5ac1681..4b6d9db3e 100644 --- a/Source/IO/ERF_Write1DProfiles.cpp +++ b/Source/IO/ERF_Write1DProfiles.cpp @@ -327,18 +327,18 @@ void ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVecto fab_arr(i, j, k,19) = p * w_cc_arr(i,j,k); // p'*w fab_arr(i, j, k,20) = cons_arr(i,j,k,RhoQ1_comp) / cons_arr(i,j,k,Rho_comp); // qv fab_arr(i, j, k,21) = cons_arr(i,j,k,RhoQ2_comp) / cons_arr(i,j,k,Rho_comp); // qc - fab_arr(i, j, k,22) = cons_arr(i,j,k,RhoQ4_comp) / cons_arr(i,j,k,Rho_comp); // qr + fab_arr(i, j, k,22) = cons_arr(i,j,k,RhoQr_comp) / cons_arr(i,j,k,Rho_comp); // qr fab_arr(i, j, k,23) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQ1_comp) / cons_arr(i,j,k,Rho_comp); // w*qv fab_arr(i, j, k,24) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQ2_comp) / cons_arr(i,j,k,Rho_comp); // w*qc - fab_arr(i, j, k,25) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQ4_comp) / cons_arr(i,j,k,Rho_comp); // w*qr - if (n_qstate > 3) { - fab_arr(i, j, k,26) = cons_arr(i,j,k,RhoQ3_comp) / cons_arr(i,j,k,Rho_comp); // qi - fab_arr(i, j, k,27) = cons_arr(i,j,k,RhoQ5_comp) / cons_arr(i,j,k,Rho_comp); // qs - fab_arr(i, j, k,28) = cons_arr(i,j,k,RhoQ6_comp) / cons_arr(i,j,k,Rho_comp); // qg + fab_arr(i, j, k,25) = w_cc_arr(i,j,k) * cons_arr(i,j,k,RhoQr_comp) / cons_arr(i,j,k,Rho_comp); // w*qr + if (n_qstate > 3) { + fab_arr(i, j, k,26) = cons_arr(i,j,k,RhoQ3_comp) / cons_arr(i,j,k,Rho_comp); // qi + fab_arr(i, j, k,27) = cons_arr(i,j,k,RhoQ5_comp) / cons_arr(i,j,k,Rho_comp); // qs + fab_arr(i, j, k,28) = cons_arr(i,j,k,RhoQ6_comp) / cons_arr(i,j,k,Rho_comp); // qg } else { - fab_arr(i, j, k,26) = 0.0; // qi - fab_arr(i, j, k,27) = 0.0; // qs - fab_arr(i, j, k,28) = 0.0; // qg + fab_arr(i, j, k,26) = 0.0; // qi + fab_arr(i, j, k,27) = 0.0; // qs + fab_arr(i, j, k,28) = 0.0; // qg } }); } // mfi From 96958c13a4fe401fee80f57b2515ffa01446cba3 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 11 May 2024 12:52:50 -0700 Subject: [PATCH 06/20] clean up to remove code duplication (#1617) --- Source/Advection/Advection.H | 7 +- Source/Advection/AdvectionSrcForState.cpp | 8 +- Source/TimeIntegration/ERF_slow_rhs_inc.cpp | 500 -------------------- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 89 +++- Source/TimeIntegration/Make.package | 4 - Source/TimeIntegration/TI_slow_headers.H | 5 +- Source/TimeIntegration/TI_slow_rhs_fun.H | 38 +- 7 files changed, 107 insertions(+), 544 deletions(-) delete mode 100644 Source/TimeIntegration/ERF_slow_rhs_inc.cpp diff --git a/Source/Advection/Advection.H b/Source/Advection/Advection.H index 1a91c2eeb..9ac750cd1 100644 --- a/Source/Advection/Advection.H +++ b/Source/Advection/Advection.H @@ -27,11 +27,8 @@ void AdvectionSrcForRho (const amrex::Box& bx, const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, - const amrex::GpuArray, AMREX_SPACEDIM>& flx_arr -#ifdef ERF_USE_POISSON_SOLVE - ,const bool const_rho = false -#endif - ); + const amrex::GpuArray, AMREX_SPACEDIM>& flx_arr, + const bool const_rho); /** Compute advection tendency for all scalars other than density and potential temperature */ void AdvectionSrcForScalars (const amrex::Box& bx, diff --git a/Source/Advection/AdvectionSrcForState.cpp b/Source/Advection/AdvectionSrcForState.cpp index bf117927c..59a01fe78 100644 --- a/Source/Advection/AdvectionSrcForState.cpp +++ b/Source/Advection/AdvectionSrcForState.cpp @@ -43,10 +43,8 @@ AdvectionSrcForRho (const Box& bx, const Array4& mf_m, const Array4& mf_u, const Array4& mf_v, - const GpuArray, AMREX_SPACEDIM>& flx_arr -#ifdef ERF_USE_POISSON_SOLVE - ,const bool const_rho -#endif + const GpuArray, AMREX_SPACEDIM>& flx_arr, + const bool const_rho ) { BL_PROFILE_VAR("AdvectionSrcForRho", AdvectionSrcForRho); @@ -73,14 +71,12 @@ AdvectionSrcForRho (const Box& bx, avg_zmom(i,j,k ) = (flx_arr[2])(i,j,k,0); }); -#ifdef ERF_USE_POISSON_SOLVE if (const_rho) { ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { advectionSrc(i,j,k,0) = 0.0; }); } else -#endif { ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp deleted file mode 100644 index 076eef8d9..000000000 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ /dev/null @@ -1,500 +0,0 @@ -#include -#include -#include -#include - -#include -#include - -using namespace amrex; - -/** - * Function for computing the slow RHS for the evolution equations for the density, potential temperature and momentum. - * - * @param[in] level level of resolution - * @param[in] nrk which RK stage - * @param[in] dt slow time step - * @param[out] S_rhs RHS computed here - * @param[in] S_old old-time solution - * @param[in] S_data current solution - * @param[in] S_prim primitive variables (i.e. conserved variables divided by density) - * @param[in] S_scratch scratch space - * @param[in] xvel x-component of velocity - * @param[in] yvel y-component of velocity - * @param[in] zvel z-component of velocity - * @param[in] qv water vapor - * @param[in] Omega component of the momentum normal to the z-coordinate surface - * @param[in] cc_src source terms for conserved variables - * @param[in] xmom_src source terms for x-momentum - * @param[in] ymom_src source terms for y-momentum - * @param[in] zmom_src source terms for z-momentum - * @param[in] Tau11 tau_11 component of stress tensor - * @param[in] Tau22 tau_22 component of stress tensor - * @param[in] Tau33 tau_33 component of stress tensor - * @param[in] Tau12 tau_12 component of stress tensor - * @param[in] Tau12 tau_13 component of stress tensor - * @param[in] Tau21 tau_21 component of stress tensor - * @param[in] Tau23 tau_23 component of stress tensor - * @param[in] Tau31 tau_31 component of stress tensor - * @param[in] Tau32 tau_32 component of stress tensor - * @param[in] SmnSmn strain rate magnitude - * @param[in] eddyDiffs diffusion coefficients for LES turbulence models - * @param[in] Hfx3 heat flux in z-dir - * @param[in] Diss dissipation of turbulent kinetic energy - * @param[in] geom Container for geometric information - * @param[in] solverChoice Container for solver parameters - * @param[in] most Pointer to MOST class for Monin-Obukhov Similarity Theory boundary condition - * @param[in] domain_bcs_type_d device vector for domain boundary conditions - * @param[in] domain_bcs_type_h host vector for domain boundary conditions - * @param[in] z_phys_nd height coordinate at nodes - * @param[in] ax area fractions on x-faces - * @param[in] ay area fractions on y-faces - * @param[in] az area fractions on z-faces - * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false) - * @param[in] p0 Reference (hydrostatically stratified) pressure - * @param[in] mapfac_m map factor at cell centers - * @param[in] mapfac_u map factor at x-faces - * @param[in] mapfac_v map factor at y-faces - * @param[in] pp_inc perturbational pressure for incompressible flows - */ - -void erf_slow_rhs_inc (int level, int nrk, - Real dt, - Vector& S_rhs, - Vector& S_old, - Vector& S_data, - const MultiFab& S_prim, - Vector& S_scratch, - const MultiFab& xvel, - const MultiFab& yvel, - const MultiFab& zvel, - MultiFab& Omega, - const MultiFab& cc_src, - const MultiFab& xmom_src, - const MultiFab& ymom_src, - const MultiFab& zmom_src, - MultiFab* Tau11, MultiFab* Tau22, MultiFab* Tau33, - MultiFab* Tau12, MultiFab* Tau13, MultiFab* Tau21, - MultiFab* Tau23, MultiFab* Tau31, MultiFab* Tau32, - MultiFab* SmnSmn, - MultiFab* eddyDiffs, - MultiFab* Hfx3, MultiFab* Diss, - const Geometry geom, - const SolverChoice& solverChoice, - std::unique_ptr& most, - const Gpu::DeviceVector& domain_bcs_type_d, - const Vector& domain_bcs_type_h, - std::unique_ptr& z_phys_nd, - std::unique_ptr& ax, - std::unique_ptr& ay, - std::unique_ptr& az, - std::unique_ptr& detJ, - const MultiFab* p0, - const MultiFab& pp_inc, - std::unique_ptr& mapfac_m, - std::unique_ptr& mapfac_u, - std::unique_ptr& mapfac_v) -{ - BL_PROFILE_REGION("erf_slow_rhs_pre_inc()"); - - const BCRec* bc_ptr_d = domain_bcs_type_d.data(); - const BCRec* bc_ptr_h = domain_bcs_type_h.data(); - - DiffChoice dc = solverChoice.diffChoice; - TurbChoice tc = solverChoice.turbChoice[level]; - - const MultiFab* t_mean_mf = nullptr; - if (most) t_mean_mf = most->get_mac_avg(0,2); - - int start_comp = 0; - int num_comp = 2; - int end_comp = start_comp + num_comp - 1; - - const bool l_const_rho = solverChoice.constant_density; - const AdvType l_horiz_adv_type = solverChoice.advChoice.dycore_horiz_adv_type; - const AdvType l_vert_adv_type = solverChoice.advChoice.dycore_vert_adv_type; - const Real l_horiz_upw_frac = solverChoice.advChoice.dycore_horiz_upw_frac; - const Real l_vert_upw_frac = solverChoice.advChoice.dycore_vert_upw_frac; - const bool l_use_terrain = solverChoice.use_terrain; - - AMREX_ALWAYS_ASSERT (!l_use_terrain); - - const bool l_use_diff = ( (dc.molec_diff_type != MolecDiffType::None) || - (tc.les_type != LESType::None) || - (tc.pbl_type != PBLType::None) ); - const bool l_use_turb = ( tc.les_type == LESType::Smagorinsky || - tc.les_type == LESType::Deardorff || - tc.pbl_type == PBLType::MYNN25 || - tc.pbl_type == PBLType::YSU ); - - const bool use_most = (most != nullptr); - const bool exp_most = (solverChoice.use_explicit_most); - - const Box& domain = geom.Domain(); - const int domhi_z = domain.bigEnd(2); - const int domlo_z = domain.smallEnd(2); - - const GpuArray dxInv = geom.InvCellSizeArray(); - - // ***************************************************************************** - // Combine external forcing terms - // ***************************************************************************** - const Array grav{0.0, 0.0, -solverChoice.gravity}; - const GpuArray grav_gpu{grav[0], grav[1], grav[2]}; - - // ***************************************************************************** - // Pre-computed quantities - // ***************************************************************************** - int nvars = S_data[IntVars::cons].nComp(); - const BoxArray& ba = S_data[IntVars::cons].boxArray(); - const DistributionMapping& dm = S_data[IntVars::cons].DistributionMap(); - - std::unique_ptr expr; - std::unique_ptr dflux_x; - std::unique_ptr dflux_y; - std::unique_ptr dflux_z; - - if (l_use_diff) { - erf_make_tau_terms(level,nrk,domain_bcs_type_d,domain_bcs_type_h,z_phys_nd, - S_data,xvel,yvel,zvel,Omega, - Tau11,Tau22,Tau33,Tau12,Tau13,Tau21,Tau23,Tau31,Tau32, - SmnSmn,eddyDiffs,Hfx3,Diss,geom,solverChoice,most, - detJ,mapfac_m,mapfac_u,mapfac_v); - - dflux_x = std::make_unique(convert(ba,IntVect(1,0,0)), dm, nvars, 0); - dflux_y = std::make_unique(convert(ba,IntVect(0,1,0)), dm, nvars, 0); - dflux_z = std::make_unique(convert(ba,IntVect(0,0,1)), dm, nvars, 0); - } // l_use_diff - - // ************************************************************************* - // Define updates and fluxes in the current RK stage - // ************************************************************************* - - // Open bc will be imposed upon all vars (we only access cons here for simplicity) - const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open); - const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open); - const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open); - const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open); - -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - { - std::array flux; - - for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) - { - Box bx = mfi.tilebox(); - Box tbx = mfi.nodaltilebox(0); - Box tby = mfi.nodaltilebox(1); - Box tbz = mfi.nodaltilebox(2); - - // If we are imposing open bc's then don't add rhs terms at the boundary locations - if ( xlo_open && (tbx.smallEnd(0) == domain.smallEnd(0)) ) {tbx.growLo(0,-1);} - if ( xhi_open && (tbx.bigEnd(0) == domain.bigEnd(0)+1) ) {tbx.growHi(0,-1);} - if ( ylo_open && (tby.smallEnd(1) == domain.smallEnd(1)) ) {tby.growLo(1,-1);} - if ( yhi_open && (tby.bigEnd(1) == domain.bigEnd(1)+1) ) {tby.growHi(1,-1);} - - // We don't compute a source term for z-momentum on the bottom or top boundary - tbz.growLo(2,-1); - tbz.growHi(2,-1); - - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4 & cell_prim = S_prim.array(mfi); - const Array4 & cell_rhs = S_rhs[IntVars::cons].array(mfi); - - const Array4 & cell_data_old = S_old[IntVars::cons].array(mfi); - - // We must initialize these to zero each RK step - S_scratch[IntVars::xmom][mfi].template setVal(0.,tbx); - S_scratch[IntVars::ymom][mfi].template setVal(0.,tby); - S_scratch[IntVars::zmom][mfi].template setVal(0.,tbz); - - Array4 avg_xmom = S_scratch[IntVars::xmom].array(mfi); - Array4 avg_ymom = S_scratch[IntVars::ymom].array(mfi); - Array4 avg_zmom = S_scratch[IntVars::zmom].array(mfi); - - const Array4 & u = xvel.array(mfi); - const Array4 & v = yvel.array(mfi); - const Array4 & w = zvel.array(mfi); - - const Array4& rho_u = S_data[IntVars::xmom].array(mfi); - const Array4& rho_v = S_data[IntVars::ymom].array(mfi); - const Array4& rho_w = S_data[IntVars::zmom].array(mfi); - - const Array4& rho_u_old = S_old[IntVars::xmom].array(mfi); - const Array4& rho_v_old = S_old[IntVars::ymom].array(mfi); - const Array4& rho_w_old = S_old[IntVars::zmom].array(mfi); - - const Array4& xmom_src_arr = xmom_src.const_array(mfi); - const Array4& ymom_src_arr = ymom_src.const_array(mfi); - const Array4& zmom_src_arr = zmom_src.const_array(mfi); - - // Perturbational pressure - const Array4& pp_arr = pp_inc.const_array(mfi); - - // Map factors - const Array4& mf_m = mapfac_m->const_array(mfi); - const Array4& mf_u = mapfac_u->const_array(mfi); - const Array4& mf_v = mapfac_v->const_array(mfi); - - const Array4< Real>& omega_arr = Omega.array(mfi); - - const Array4& rho_u_rhs = S_rhs[IntVars::xmom].array(mfi); - const Array4& rho_v_rhs = S_rhs[IntVars::ymom].array(mfi); - const Array4& rho_w_rhs = S_rhs[IntVars::zmom].array(mfi); - - const Array4& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : Array4{}; - - // Terrain metrics - const Array4& z_nd = l_use_terrain ? z_phys_nd->const_array(mfi) : Array4{}; - const Array4& detJ_arr = detJ->const_array(mfi); - - // Base state - // const Array4& p0_arr = p0->const_array(mfi); - - // ***************************************************************************** - // ***************************************************************************** - // Define flux arrays for use in advection - // ***************************************************************************** - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - flux[dir].resize(surroundingNodes(bx,dir),2); - flux[dir].setVal(0.); - } - const GpuArray, AMREX_SPACEDIM> - flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}}; - - // ***************************************************************************** - // Contravariant flux field - // ***************************************************************************** - { - BL_PROFILE("slow_rhs_making_omega"); - Box gbxo = surroundingNodes(bx,2); gbxo.grow(IntVect(1,1,0)); - // Now create Omega with momentum (not velocity) - ParallelFor(gbxo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - omega_arr(i,j,k) = rho_w(i,j,k); - }); - } // end profile - - - // ***************************************************************************** - // Diffusive terms (pre-computed above) - // ***************************************************************************** - // Expansion - Array4 er_arr; - if (expr) { - er_arr = expr->array(mfi); - } else { - er_arr = Array4{}; - } - - // No terrain diffusion - Array4 tau11,tau22,tau33; - Array4 tau12,tau13,tau23; - if (Tau11) { - tau11 = Tau11->array(mfi); tau22 = Tau22->array(mfi); tau33 = Tau33->array(mfi); - tau12 = Tau12->array(mfi); tau13 = Tau13->array(mfi); tau23 = Tau23->array(mfi); - } else { - tau11 = Array4{}; tau22 = Array4{}; tau33 = Array4{}; - tau12 = Array4{}; tau13 = Array4{}; tau23 = Array4{}; - } - // Terrain diffusion - Array4 tau21,tau31,tau32; - if (Tau21) { - tau21 = Tau21->array(mfi); tau31 = Tau31->array(mfi); tau32 = Tau32->array(mfi); - } else { - tau21 = Array4{}; tau31 = Array4{}; tau32 = Array4{}; - } - - // Strain magnitude - Array4 SmnSmn_a; - if (tc.les_type == LESType::Deardorff) { - SmnSmn_a = SmnSmn->array(mfi); - } else { - SmnSmn_a = Array4{}; - } - - // ************************************************************************** - // Define updates in the RHS of continuity and potential temperature equations - // ************************************************************************** - auto const& ax_arr = ax->const_array(mfi); - auto const& ay_arr = ay->const_array(mfi); - auto const& az_arr = az->const_array(mfi); - - AdvectionSrcForRho(bx, cell_rhs, - rho_u, rho_v, omega_arr, // these are being used to build the fluxes - avg_xmom, avg_ymom, avg_zmom, // these are being defined from the fluxes - ax_arr, ay_arr, az_arr, detJ_arr, - dxInv, mf_m, mf_u, mf_v, - flx_arr, l_const_rho); - - int icomp = RhoTheta_comp; int ncomp = 1; - AdvectionSrcForScalars(bx, icomp, ncomp, - avg_xmom, avg_ymom, avg_zmom, - cell_prim, cell_rhs, detJ_arr, - dxInv, mf_m, - l_horiz_adv_type, l_vert_adv_type, - l_horiz_upw_frac, l_vert_upw_frac, - flx_arr, domain, bc_ptr_h); - - if (l_use_diff) { - Array4 diffflux_x = dflux_x->array(mfi); - Array4 diffflux_y = dflux_y->array(mfi); - Array4 diffflux_z = dflux_z->array(mfi); - - Array4 hfx_z = Hfx3->array(mfi); - Array4 diss = Diss->array(mfi); - - const Array4 tm_arr = t_mean_mf ? t_mean_mf->const_array(mfi) : Array4{}; - - // NOTE: No diffusion for continuity, so n starts at 1. - int n_start = amrex::max(start_comp,RhoTheta_comp); - int n_comp = end_comp - n_start + 1; - - DiffusionSrcForState_N(bx, domain, n_start, n_comp, exp_most, u, v, - cell_data, cell_prim, cell_rhs, - diffflux_x, diffflux_y, diffflux_z, - dxInv, SmnSmn_a, mf_m, mf_u, mf_v, - hfx_z, diss, - mu_turb, dc, tc, - tm_arr, grav_gpu, bc_ptr_d, use_most); - } - - // Add source terms for (rho theta) - { - auto const& src_arr = cc_src.const_array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - cell_rhs(i,j,k,RhoTheta_comp) += src_arr(i,j,k,RhoTheta_comp); - }); - } - - // If in second RK stage, take average of old-time and new-time source - if (nrk == 1) - { - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - cell_rhs(i,j,k, Rho_comp) *= 0.5; - cell_rhs(i,j,k,RhoTheta_comp) *= 0.5; - - cell_rhs(i,j,k, Rho_comp) += 0.5 / dt * (cell_data(i,j,k, Rho_comp) - cell_data_old(i,j,k, Rho_comp)); - cell_rhs(i,j,k,RhoTheta_comp) += 0.5 / dt * (cell_data(i,j,k,RhoTheta_comp) - cell_data_old(i,j,k,RhoTheta_comp)); - }); - } - - // ***************************************************************************** - // Define updates in the RHS of {x, y, z}-momentum equations - // ***************************************************************************** - int lo_z_face; - int hi_z_face; - if (level == 0) { - lo_z_face = domain.smallEnd(2); - hi_z_face = domain.bigEnd(2)+1; - } else { - lo_z_face = mfi.validbox().smallEnd(2); - hi_z_face = mfi.validbox().bigEnd(2)+1; - } - AdvectionSrcForMom(bx, tbx, tby, tbz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - cell_data, u, v, w, - rho_u, rho_v, omega_arr, - z_nd, ax_arr, ay_arr, az_arr, detJ_arr, - dxInv, mf_m, mf_u, mf_v, - l_horiz_adv_type, l_vert_adv_type, - l_horiz_upw_frac, l_vert_upw_frac, - l_use_terrain, lo_z_face, hi_z_face, - domain, bc_ptr_h); - - if (l_use_diff) { - DiffusionSrcForMom_N(tbx, tby, tbz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - tau11, tau22, tau33, - tau12, tau13, tau23, - dxInv, - mf_m, mf_u, mf_v); - } - - ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { // x-momentum equation - - Real gpx = dxInv[0] * (pp_arr(i,j,k) - pp_arr(i-1,j,k)); - - rho_u_rhs(i, j, k) += xmom_src_arr(i,j,k) - gpx - solverChoice.abl_pressure_grad[0]; - - if (nrk == 1) { - rho_u_rhs(i,j,k) *= 0.5; - rho_u_rhs(i,j,k) += 0.5 / dt * (rho_u(i,j,k) - rho_u_old(i,j,k)); - } - }); - - ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { // y-momentum equation - - Real gpy = dxInv[1] * (pp_arr(i,j,k) - pp_arr(i,j-1,k)); - - rho_v_rhs(i, j, k) += ymom_src_arr(i,j,k) - gpy - solverChoice.abl_pressure_grad[1]; - - 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)); - } - }); - - // ***************************************************************************** - // Zero out source terms for x- and y- momenta if at walls or inflow - // We need to do this -- even though we call the boundary conditions later -- - // because the slow source is used to update the state in the fast interpolater. - // ***************************************************************************** - if ( (bx.smallEnd(0) == domain.smallEnd(0)) && - (bc_ptr_h[BCVars::xvel_bc].lo(0) == ERFBCType::ext_dir) ) { - Box lo_x_dom_face(bx); lo_x_dom_face.setBig(0,bx.smallEnd(0)); - ParallelFor(lo_x_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - rho_u_rhs(i,j,k) = 0.; - }); - } - if ( (bx.bigEnd(0) == domain.bigEnd(0)) && - (bc_ptr_h[BCVars::xvel_bc].hi(0) == ERFBCType::ext_dir) ) { - Box hi_x_dom_face(bx); hi_x_dom_face.setSmall(0,bx.bigEnd(0)+1); hi_x_dom_face.setBig(0,bx.bigEnd(0)+1); - ParallelFor(hi_x_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - rho_u_rhs(i,j,k) = 0.; - }); - } - if ( (bx.smallEnd(1) == domain.smallEnd(1)) && - (bc_ptr_h[BCVars::yvel_bc].lo(1) == ERFBCType::ext_dir) ) { - Box lo_y_dom_face(bx); lo_y_dom_face.setBig(1,bx.smallEnd(1)); - ParallelFor(lo_y_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - rho_v_rhs(i,j,k) = 0.; - }); - } - if ( (bx.bigEnd(1) == domain.bigEnd(1)) && - (bc_ptr_h[BCVars::yvel_bc].hi(1) == ERFBCType::ext_dir) ) { - Box hi_y_dom_face(bx); hi_y_dom_face.setSmall(1,bx.bigEnd(1)+1); hi_y_dom_face.setBig(1,bx.bigEnd(1)+1);; - ParallelFor(hi_y_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - rho_v_rhs(i,j,k) = 0.; - }); - } - - amrex::Box b2d = tbz; - b2d.setSmall(2,0); - b2d.setBig(2,0); - // Enforce no forcing term at top and bottom boundaries - ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) { - rho_w_rhs(i,j, 0) = 0.; - rho_w_rhs(i,j,domhi_z+1) = 0.; - }); - - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { // z-momentum equation - - Real gpz = dxInv[2] * (pp_arr(i,j,k) - pp_arr(i,j,k-1)); - - rho_w_rhs(i, j, k) += zmom_src_arr(i,j,k) - gpz - solverChoice.abl_pressure_grad[2]; - - 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)); - } - }); - } // mfi - } // omp -} diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 9fb3a964c..c21213fa8 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -17,6 +17,7 @@ using namespace amrex; * @param[in] nrk which RK stage * @param[in] dt slow time step * @param[out] S_rhs RHS computed here + * @param[in] S_old old-time solution -- used only for incompressible * @param[in] S_data current solution * @param[in] S_prim primitive variables (i.e. conserved variables divided by density) * @param[in] S_scratch scratch space @@ -66,6 +67,7 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk, Real dt, Vector& S_rhs, + Vector& S_old, Vector& S_data, const MultiFab& S_prim, Vector& S_scratch, @@ -148,6 +150,18 @@ void erf_slow_rhs_pre (int level, int finest_level, const bool use_most = (most != nullptr); const bool exp_most = (solverChoice.use_explicit_most); +#ifdef ERF_USE_POISSON_SOlVE + const bool l_incompressible = solverChoice.incompressible; + const bool l_const_rho = solverChoice.constant_density; + + // We cannot use incompressible with terrain or with moisture + AMREX_ALWAYS_ASSERT(!l_use_terrain || !incompressible); + AMREX_ALWAYS_ASSERT(!l_use_moisture || !incompressible); +#else + const bool l_incompressible = false; + const bool l_const_rho = false; +#endif + const Box& domain = geom.Domain(); const int domhi_z = domain.bigEnd(2); const int domlo_z = domain.smallEnd(2); @@ -222,10 +236,23 @@ void erf_slow_rhs_pre (int level, int finest_level, const Array4 & cell_prim = S_prim.array(mfi); const Array4 & cell_rhs = S_rhs[IntVars::cons].array(mfi); + const Array4 & cell_old = S_old[IntVars::cons].array(mfi); + const Array4& xmom_src_arr = xmom_src.const_array(mfi); const Array4& ymom_src_arr = ymom_src.const_array(mfi); const Array4& zmom_src_arr = zmom_src.const_array(mfi); + const Array4& rho_u_old = S_old[IntVars::xmom].array(mfi); + const Array4& rho_v_old = S_old[IntVars::ymom].array(mfi); + const Array4& rho_w_old = S_old[IntVars::zmom].array(mfi); + + if (l_incompressible) { + // When incompressible we must reset these to 0 each RK step + S_scratch[IntVars::xmom][mfi].template setVal(0.0,tbx); + S_scratch[IntVars::ymom][mfi].template setVal(0.0,tby); + S_scratch[IntVars::zmom][mfi].template setVal(0.0,tbz); + } + Array4 avg_xmom = S_scratch[IntVars::xmom].array(mfi); Array4 avg_ymom = S_scratch[IntVars::ymom].array(mfi); Array4 avg_zmom = S_scratch[IntVars::zmom].array(mfi); @@ -277,21 +304,26 @@ void erf_slow_rhs_pre (int level, int finest_level, // ***************************************************************************** // Perturbational pressure field // ***************************************************************************** - Box gbx = mfi.tilebox(); gbx.grow(IntVect(1,1,1)); - if (gbx.smallEnd(2) < 0) gbx.setSmall(2,0); - FArrayBox pprime; pprime.resize(gbx,1,The_Async_Arena()); - const Array4 & pp_arr = pprime.array(); - { - BL_PROFILE("slow_rhs_pre_pprime"); - ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - //if (cell_data(i,j,k,RhoTheta_comp) < 0.) printf("BAD THETA AT %d %d %d %e %e \n", - // i,j,k,cell_data(i,j,k,RhoTheta_comp),cell_data(i,j,k+1,RhoTheta_comp)); - AMREX_ASSERT(cell_data(i,j,k,RhoTheta_comp) > 0.); - Real qv_for_p = (use_moisture) ? cell_data(i,j,k,RhoQ1_comp)/cell_data(i,j,k,Rho_comp) : 0.0; - pp_arr(i,j,k) = getPgivenRTh(cell_data(i,j,k,RhoTheta_comp),qv_for_p) - p0_arr(i,j,k); - }); - } // end profile + FArrayBox pprime; + if (!l_incompressible) { + Box gbx = mfi.tilebox(); gbx.grow(IntVect(1,1,1)); + if (gbx.smallEnd(2) < 0) gbx.setSmall(2,0); + pprime.resize(gbx,1,The_Async_Arena()); + const Array4& pptemp_arr = pprime.array(); + ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + //if (cell_data(i,j,k,RhoTheta_comp) < 0.) printf("BAD THETA AT %d %d %d %e %e \n", + // i,j,k,cell_data(i,j,k,RhoTheta_comp),cell_data(i,j,k+1,RhoTheta_comp)); + AMREX_ASSERT(cell_data(i,j,k,RhoTheta_comp) > 0.); + Real qv_for_p = (use_moisture) ? cell_data(i,j,k,RhoQ1_comp)/cell_data(i,j,k,Rho_comp) : 0.0; + pptemp_arr(i,j,k) = getPgivenRTh(cell_data(i,j,k,RhoTheta_comp),qv_for_p) - p0_arr(i,j,k); + }); + } +#ifdef ERF_USE_POISSON_SOLVE + const Array4& pp_arr = (l_incompressible) ? pp_inc.const_array(mfi) : pprime.const_array(); +#else + const Array4& pp_arr = pprime.const_array(); +#endif // ***************************************************************************** // Contravariant flux field @@ -393,7 +425,7 @@ void erf_slow_rhs_pre (int level, int finest_level, avg_xmom, avg_ymom, avg_zmom, // these are being defined from the fluxes ax_arr, ay_arr, az_arr, detJ_arr, dxInv, mf_m, mf_u, mf_v, - flx_arr); + flx_arr, l_const_rho); int icomp = RhoTheta_comp; int ncomp = 1; AdvectionSrcForScalars(bx, icomp, ncomp, @@ -453,6 +485,19 @@ void erf_slow_rhs_pre (int level, int finest_level, }); } + // If incompressible and in second RK stage, take average of old-time and new-time source + if ( l_incompressible && (nrk == 1) ) + { + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + cell_rhs(i,j,k, Rho_comp) *= 0.5; + cell_rhs(i,j,k,RhoTheta_comp) *= 0.5; + + cell_rhs(i,j,k, Rho_comp) += 0.5 / dt * (cell_data(i,j,k, Rho_comp) - cell_old(i,j,k, Rho_comp)); + cell_rhs(i,j,k,RhoTheta_comp) += 0.5 / dt * (cell_data(i,j,k,RhoTheta_comp) - cell_old(i,j,k,RhoTheta_comp)); + }); + } + // ***************************************************************************** // Define updates in the RHS of {x, y, z}-momentum equations // ***************************************************************************** @@ -546,6 +591,11 @@ void erf_slow_rhs_pre (int level, int finest_level, Real h_zeta = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd); rho_u_rhs(i, j, k) *= h_zeta; } + + if ( l_incompressible && (nrk == 1) ) { + rho_u_rhs(i,j,k) *= 0.5; + rho_u_rhs(i,j,k) += 0.5 / dt * (rho_u(i,j,k) - rho_u_old(i,j,k)); + } }); ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -589,6 +639,11 @@ void erf_slow_rhs_pre (int level, int finest_level, Real h_zeta = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd); rho_v_rhs(i, j, k) *= h_zeta; } + + if ( l_incompressible && (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)); + } }); // ***************************************************************************** @@ -666,7 +721,7 @@ void erf_slow_rhs_pre (int level, int finest_level, // NOTE: for now we are only refluxing density not (rho theta) since the latter seems to introduce // a problem at top and bottom boundaries if (l_reflux && nrk == 2) { - int strt_comp_reflux = 0; + int strt_comp_reflux = (l_const_rho) ? 1 : 0; int num_comp_reflux = 1; if (level < finest_level) { fr_as_crse->CrseAdd(mfi, diff --git a/Source/TimeIntegration/Make.package b/Source/TimeIntegration/Make.package index 8a845a8cc..7f38f3eab 100644 --- a/Source/TimeIntegration/Make.package +++ b/Source/TimeIntegration/Make.package @@ -12,10 +12,6 @@ CEXE_sources += ERF_fast_rhs_N.cpp CEXE_sources += ERF_fast_rhs_T.cpp CEXE_sources += ERF_fast_rhs_MT.cpp -ifeq ($(USE_POISSON_SOLVE),TRUE) -CEXE_sources += ERF_slow_rhs_inc.cpp -endif - CEXE_headers += TI_fast_rhs_fun.H CEXE_headers += TI_slow_rhs_fun.H CEXE_headers += TI_no_substep_fun.H diff --git a/Source/TimeIntegration/TI_slow_headers.H b/Source/TimeIntegration/TI_slow_headers.H index c3d350429..dafbf5f1b 100644 --- a/Source/TimeIntegration/TI_slow_headers.H +++ b/Source/TimeIntegration/TI_slow_headers.H @@ -57,9 +57,10 @@ void erf_make_tau_terms (int level, int nrk, void erf_slow_rhs_pre (int level, int finest_level, int nrk, amrex::Real dt, amrex::Vector& S_rhs, + amrex::Vector& S_old, amrex::Vector& S_data, - const amrex::MultiFab& S_prim, - amrex::Vector& S_scratch, + const amrex::MultiFab & S_prim, + amrex::Vector& S_scratch, const amrex::MultiFab& xvel, const amrex::MultiFab& yvel, const amrex::MultiFab& zvel, diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 23694da64..54a709392 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -4,11 +4,7 @@ * Wrapper for calling the routine that creates the slow RHS */ auto slow_rhs_fun_pre = [&](Vector& S_rhs, -#ifdef ERF_USE_NETCDF Vector& S_old, -#else - Vector& /*S_old*/, -#endif Vector& S_data, Vector& S_scratch, const Real old_step_time, @@ -119,7 +115,7 @@ dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, n_qstate); - erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, + erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, z_t_rk[level], Omega, cc_src, xmom_src, ymom_src, zmom_src, Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(), @@ -222,7 +218,7 @@ dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, n_qstate); - erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, + erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, z_t_rk[level], Omega, cc_src, xmom_src, ymom_src, zmom_src, Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(), @@ -392,6 +388,21 @@ Real slow_dt = new_stage_time - old_step_time; + // ************************************************************************* + // Set up flux registers if using two_way coupling + // ************************************************************************* + YAFluxRegister* fr_as_crse = nullptr; + YAFluxRegister* fr_as_fine = nullptr; + if (solverChoice.coupling_type == CouplingType::TwoWay) { + if (level < finest_level) { + fr_as_crse = getAdvFluxReg(level+1); + fr_as_crse->reset(); + } + if (level > 0) { + fr_as_fine = getAdvFluxReg(level); + } + } + 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; @@ -413,16 +424,23 @@ dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, n_qstate); - erf_slow_rhs_inc(level, nrk, slow_dt, + erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, - Omega, cc_src, xmom_src, ymom_src, zmom_src, + z_t_rk[level], Omega, cc_src, xmom_src, ymom_src, zmom_src, Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(), Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(), Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx3, Diss, fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, - z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0, pp_inc[level], - mapfac_m[level], mapfac_u[level], mapfac_v[level]); + z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0, +#ifdef ERF_USE_POISSON_SOLVE + pp_inc[level], +#endif + mapfac_m[level], mapfac_u[level], mapfac_v[level], +#ifdef ERF_USE_EB + EBFactory(level), +#endif + fr_as_crse, fr_as_fine); add_thin_body_sources(xmom_src, ymom_src, zmom_src, xflux_imask[level], yflux_imask[level], zflux_imask[level], From 4bed605557d412ab0a94af14e71a9c8e6a25303c Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 11 May 2024 16:42:42 -0700 Subject: [PATCH 07/20] fix typo in ifdef for incompressible (#1618) --- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 28 ++++++++++----------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index c21213fa8..fce69e346 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -146,17 +146,17 @@ void erf_slow_rhs_pre (int level, int finest_level, tc.pbl_type == PBLType::MYNN25 || tc.pbl_type == PBLType::YSU ); - const bool use_moisture = (solverChoice.moisture_type != MoistureType::None); - const bool use_most = (most != nullptr); - const bool exp_most = (solverChoice.use_explicit_most); + const bool l_use_moisture = (solverChoice.moisture_type != MoistureType::None); + const bool l_use_most = (most != nullptr); + const bool l_exp_most = (solverChoice.use_explicit_most); -#ifdef ERF_USE_POISSON_SOlVE +#ifdef ERF_USE_POISSON_SOLVE const bool l_incompressible = solverChoice.incompressible; const bool l_const_rho = solverChoice.constant_density; // We cannot use incompressible with terrain or with moisture - AMREX_ALWAYS_ASSERT(!l_use_terrain || !incompressible); - AMREX_ALWAYS_ASSERT(!l_use_moisture || !incompressible); + AMREX_ALWAYS_ASSERT(!l_use_terrain || !l_incompressible); + AMREX_ALWAYS_ASSERT(!l_use_moisture || !l_incompressible); #else const bool l_incompressible = false; const bool l_const_rho = false; @@ -315,7 +315,7 @@ void erf_slow_rhs_pre (int level, int finest_level, //if (cell_data(i,j,k,RhoTheta_comp) < 0.) printf("BAD THETA AT %d %d %d %e %e \n", // i,j,k,cell_data(i,j,k,RhoTheta_comp),cell_data(i,j,k+1,RhoTheta_comp)); AMREX_ASSERT(cell_data(i,j,k,RhoTheta_comp) > 0.); - Real qv_for_p = (use_moisture) ? cell_data(i,j,k,RhoQ1_comp)/cell_data(i,j,k,Rho_comp) : 0.0; + Real qv_for_p = (l_use_moisture) ? cell_data(i,j,k,RhoQ1_comp)/cell_data(i,j,k,Rho_comp) : 0.0; pptemp_arr(i,j,k) = getPgivenRTh(cell_data(i,j,k,RhoTheta_comp),qv_for_p) - p0_arr(i,j,k); }); } @@ -451,21 +451,21 @@ void erf_slow_rhs_pre (int level, int finest_level, int n_comp = end_comp - n_start + 1; if (l_use_terrain) { - DiffusionSrcForState_T(bx, domain, n_start, n_comp, exp_most, u, v, + DiffusionSrcForState_T(bx, domain, n_start, n_comp, l_exp_most, u, v, cell_data, cell_prim, cell_rhs, diffflux_x, diffflux_y, diffflux_z, z_nd, ax_arr, ay_arr, az_arr, detJ_arr, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_z, diss, mu_turb, dc, tc, - tm_arr, grav_gpu, bc_ptr_d, use_most); + tm_arr, grav_gpu, bc_ptr_d, l_use_most); } else { - DiffusionSrcForState_N(bx, domain, n_start, n_comp, exp_most, u, v, + DiffusionSrcForState_N(bx, domain, n_start, n_comp, l_exp_most, u, v, cell_data, cell_prim, cell_rhs, diffflux_x, diffflux_y, diffflux_z, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_z, diss, mu_turb, dc, tc, - tm_arr, grav_gpu, bc_ptr_d, use_most); + tm_arr, grav_gpu, bc_ptr_d, l_use_most); } } @@ -579,7 +579,7 @@ void erf_slow_rhs_pre (int level, int finest_level, gpx *= mf_u(i,j,0); Real q = 0.0; - if (use_moisture) { + if (l_use_moisture) { q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i-1,j,k,PrimQ1_comp) +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i-1,j,k,PrimQ2_comp) ); } @@ -628,7 +628,7 @@ void erf_slow_rhs_pre (int level, int finest_level, gpy *= mf_v(i,j,0); Real q = 0.0; - if (use_moisture) { + if (l_use_moisture) { 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) ); } @@ -704,7 +704,7 @@ void erf_slow_rhs_pre (int level, int finest_level, Real gpz = dxInv[2] * ( pp_arr(i,j,k)-pp_arr(i,j,k-1) ) / met_h_zeta; Real q = 0.0; - if (use_moisture) { + if (l_use_moisture) { q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i,j,k-1,PrimQ1_comp) +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j,k-1,PrimQ2_comp) ); } From 8f95c95eb9825ce32f745fa9ece4822f19ccd24e Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 13 May 2024 17:19:07 -0700 Subject: [PATCH 08/20] update amrex submodule to fix HIP CI (#1619) --- Submodules/AMReX | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Submodules/AMReX b/Submodules/AMReX index a89b4651c..af6e77226 160000 --- a/Submodules/AMReX +++ b/Submodules/AMReX @@ -1 +1 @@ -Subproject commit a89b4651c99e43489185ea26ac18ae048bce4051 +Subproject commit af6e77226cdb1668a8936ee3b35cf0d441409869 From e8d31360566389c0a03c7929f093d975936ac97d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Tue, 14 May 2024 14:22:03 -0700 Subject: [PATCH 09/20] Add the viscous terms for open bc boundaries (#1615) Co-authored-by: Mahesh Natarajan Co-authored-by: Ann Almgren --- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index fce69e346..188e27dbc 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -222,12 +222,6 @@ void erf_slow_rhs_pre (int level, int finest_level, Box tby = mfi.nodaltilebox(1); Box tbz = mfi.nodaltilebox(2); - // If we are imposing open bc's then don't add rhs terms at the boundary locations - if ( xlo_open && (tbx.smallEnd(0) == domain.smallEnd(0)) ) {tbx.growLo(0,-1);} - if ( xhi_open && (tbx.bigEnd(0) == domain.bigEnd(0)+1) ) {tbx.growHi(0,-1);} - if ( ylo_open && (tby.smallEnd(1) == domain.smallEnd(1)) ) {tby.growLo(1,-1);} - if ( yhi_open && (tby.bigEnd(1) == domain.bigEnd(1)+1) ) {tby.growHi(1,-1);} - // We don't compute a source term for z-momentum on the bottom or top boundary tbz.growLo(2,-1); tbz.growHi(2,-1); From 6fe430bbf493fdec4d8bba7d3279ec18d3ecf772 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 14 May 2024 18:59:18 -0700 Subject: [PATCH 10/20] move windfarm details into Source/WindFarm (#1620) --- Exec/Make.ERF | 8 +++- Source/DataStructs/DataStruct.H | 46 ++++++++++++------- Source/ERF.H | 4 +- Source/ERF.cpp | 8 ++-- Source/ERF_make_new_arrays.cpp | 12 ++--- Source/TimeIntegration/ERF_Advance.cpp | 14 ++---- Source/TimeIntegration/ERF_make_tau_terms.cpp | 21 --------- Source/TimeIntegration/ERF_slow_rhs_post.cpp | 2 +- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 8 ++-- Source/TimeIntegration/TI_slow_headers.H | 3 -- .../EWP/AdvanceEWP.cpp | 12 ++--- Source/WindFarmParametrization/Make.package | 2 + Source/WindFarmParametrization/WindFarm.H | 16 +++++++ Source/WindFarmParametrization/WindFarm.cpp | 22 +++++++++ 14 files changed, 100 insertions(+), 78 deletions(-) create mode 100644 Source/WindFarmParametrization/Make.package create mode 100644 Source/WindFarmParametrization/WindFarm.H create mode 100644 Source/WindFarmParametrization/WindFarm.cpp diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 1b52fbb37..7efb8b1fa 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -141,10 +141,14 @@ INCLUDE_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) # at runtime from the inputs ifeq ($(USE_WINDFARM), TRUE) DEFINES += -DERF_USE_WINDFARM - ERF_WINDFARM_FITCH_DIR = $(ERF_SOURCE_DIR)/WindFarmParametrization/Fitch - ERF_WINDFARM_EWP_DIR = $(ERF_SOURCE_DIR)/WindFarmParametrization/EWP + ERF_WINDFARM_DIR = $(ERF_SOURCE_DIR)/WindFarmParametrization + ERF_WINDFARM_FITCH_DIR = $(ERF_WINDFARM_DIR)/Fitch + ERF_WINDFARM_EWP_DIR = $(ERF_WINDFARM_DIR)/EWP + include $(ERF_WINDFARM_DIR)/Make.package include $(ERF_WINDFARM_FITCH_DIR)/Make.package include $(ERF_WINDFARM_EWP_DIR)/Make.package + VPATH_LOCATIONS += $(ERF_WINDFARM_DIR) + INCLUDE_LOCATIONS += $(ERF_WINDFARM_DIR) VPATH_LOCATIONS += $(ERF_WINDFARM_FITCH_DIR) INCLUDE_LOCATIONS += $(ERF_WINDFARM_FITCH_DIR) VPATH_LOCATIONS += $(ERF_WINDFARM_EWP_DIR) diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 5a092b237..b738ec678 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -37,7 +37,7 @@ enum struct MoistureType { }; enum struct WindFarmType { - Fitch, EWP + Fitch, EWP, None }; enum struct LandSurfaceType { @@ -161,12 +161,24 @@ struct SolverChoice { // Should we project the initial velocity field to make it divergence-free? pp.query("project_initial_velocity", project_initial_velocity); - pp.query("incompressible", incompressible); + int nvals_inc = pp.countval("incompressible"); + AMREX_ALWAYS_ASSERT(nvals_inc == 1 || nvals_inc >= max_level+1); + amrex::Vector inc_in; inc_in.resize(nvals_inc); + pp.getarr("incompressible",inc_in); + if (nvals_inc == 1) { + for (int i = 0; i <= max_level; ++i) incompressible.push_back(inc_in[0]); + } else { + for (int i = 0; i <= max_level; ++i) incompressible.push_back(inc_in[0]); + } + pp.query("constant_density", constant_density); pp.query("project_every_stage", project_every_stage); pp.query("ncorr", ncorr); pp.query("poisson_abstol", poisson_abstol); pp.query("poisson_reltol", poisson_reltol); +#else + incompressible.resize(max_level+1); + for (int i = 0; i <= max_level; ++i) incompressible[i] = 0; #endif // Turn off acoustic substepping? @@ -174,10 +186,11 @@ struct SolverChoice { pp.query("force_stage1_single_substep", force_stage1_single_substep); #if defined(ERF_USE_POISSON_SOLVE) - // If this is set, it must be even - if (incompressible != 0 && no_substepping == 0) - { - amrex::Abort("If you specify incompressible, you must specific no_substepping"); + for (int lev = 0; lev <= max_level; lev++) { + if (incompressible[lev] != 0 && no_substepping == 0) + { + amrex::Abort("If you specify incompressible, you must specific no_substepping"); + } } #endif @@ -274,19 +287,19 @@ struct SolverChoice { amrex::Abort("Dont know this coupling_type"); } - // Which type of windfarm model - static std::string windfarm_type_string = "Fitch"; - - pp.query("windfarm_type", windfarm_type_string); pp.query("latitude_lo", latitude_lo); pp.query("longitude_lo", longitude_lo); + + // Which type of windfarm model + static std::string windfarm_type_string = "None"; + pp.query("windfarm_type", windfarm_type_string); if (windfarm_type_string == "Fitch") { windfarm_type = WindFarmType::Fitch; } - else if(windfarm_type_string == "EWP"){ + else if (windfarm_type_string == "EWP") { windfarm_type = WindFarmType::EWP; } - else { + else if (windfarm_type_string != "None") { amrex::Abort("Dont know this windfarm_type. Currently only Fitch and EWP models are supported."); } @@ -294,12 +307,14 @@ struct SolverChoice { pp.query("time_avg_vel",time_avg_vel); } - void display() + void display(int max_level) { amrex::Print() << "SOLVER CHOICE: " << std::endl; amrex::Print() << "no_substepping : " << no_substepping << std::endl; amrex::Print() << "force_stage1_single_substep : " << force_stage1_single_substep << std::endl; - amrex::Print() << "incompressible : " << incompressible << std::endl; + for (int lev = 0; lev <= max_level; lev++) { + amrex::Print() << "incompressible at level : " << lev << " is " << incompressible[lev] << std::endl; + } amrex::Print() << "use_coriolis : " << use_coriolis << std::endl; amrex::Print() << "use_gravity : " << use_gravity << std::endl; @@ -336,7 +351,6 @@ struct SolverChoice { diffChoice.display(); spongeChoice.display(); - int max_level = turbChoice.size()-1; for (int lev = 0; lev <= max_level; lev++) { turbChoice[lev].display(lev); } @@ -387,7 +401,7 @@ struct SolverChoice { int no_substepping = 0; int force_stage1_single_substep = 1; - int incompressible = 0; + amrex::Vector incompressible; int constant_density = 0; int project_every_stage = 1; int ncorr = 1; diff --git a/Source/ERF.H b/Source/ERF.H index 5beb392c3..e299cb61c 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -631,8 +631,8 @@ private: amrex::Vector> qmoist; // (lev,ncomp) This has up to 8 components: qt, qv, qc, qi, qp, qr, qs, qg amrex::Vector Nturb; - amrex::Vector vars_fitch; // Vabs, Vabsdt, dudt, dvdt, dTKEdt - amrex::Vector vars_ewp; // dudt, dvdt, dTKEdt + amrex::Vector vars_windfarm; // Fitch: Vabs, Vabsdt, dudt, dvdt, dTKEdt + // EWP: dudt, dvdt, dTKEdt amrex::Real hub_height, rotor_dia, thrust_coeff_standing, nominal_power; amrex::Vector wind_speed, thrust_coeff, power; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index c6215973f..c0f1effb5 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -119,8 +119,7 @@ ERF::ERF () #ifdef ERF_USE_WINDFARM Nturb.resize(nlevs_max); - vars_fitch.resize(nlevs_max); - vars_ewp.resize(nlevs_max); + vars_windfarm.resize(nlevs_max); #endif #if defined(ERF_USE_RRTMGP) @@ -1456,7 +1455,7 @@ ERF::ReadParameters () } if (verbose > 0) { - solverChoice.display(); + solverChoice.display(max_level); } if (solverChoice.coupling_type == CouplingType::TwoWay && cf_width > 0) { @@ -1803,8 +1802,7 @@ ERF::ERF (const RealBox& rb, int max_level_in, #ifdef ERF_USE_WINDFARM Nturb.resize(nlevs_max); - vars_fitch.resize(nlevs_max); - vars_ewp.resize(nlevs_max); + vars_windfarm.resize(nlevs_max); #endif #if defined(ERF_USE_RRTMGP) diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 8dbda2bff..393d1c446 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -236,14 +236,12 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, // Variables for Ftich model for windfarm parametrization //********************************************************* if (solverChoice.windfarm_type == WindFarmType::Fitch){ - int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1; - vars_fitch[lev].define(ba, dm, 5, ngrow_state); // V, dVabsdt, dudt, dvdt, dTKEdt - Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell + vars_windfarm[lev].define(ba, dm, 5, ngrow_state); // V, dVabsdt, dudt, dvdt, dTKEdt + Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell } if (solverChoice.windfarm_type == WindFarmType::EWP){ - int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1; - vars_ewp[lev].define(ba, dm, 3, ngrow_state); // dudt, dvdt, dTKEdt - Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell + vars_windfarm[lev].define(ba, dm, 3, ngrow_state); // dudt, dvdt, dTKEdt + Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell } #endif @@ -411,7 +409,7 @@ ERF::initialize_integrator (int lev, MultiFab& cons_mf, MultiFab& vel_mf) mri_integrator_mem[lev] = std::make_unique > >(int_state); mri_integrator_mem[lev]->setNoSubstepping(solverChoice.no_substepping); - mri_integrator_mem[lev]->setIncompressible(solverChoice.incompressible); + mri_integrator_mem[lev]->setIncompressible(solverChoice.incompressible[lev]); mri_integrator_mem[lev]->setNcompCons(ncomp_cons); mri_integrator_mem[lev]->setForceFirstStageSingleSubstep(solverChoice.force_stage1_single_substep); } diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index eff2ea0a8..60fd99139 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -2,8 +2,7 @@ #include #ifdef ERF_USE_WINDFARM -#include -#include +#include #endif using namespace amrex; @@ -70,14 +69,9 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) } #if defined(ERF_USE_WINDFARM) - // Update with the Fitch source terms - if (solverChoice.windfarm_type == WindFarmType::Fitch) { - fitch_advance(lev, Geom(lev), dt_lev, S_old, - U_old, V_old, W_old, vars_fitch[lev], Nturb[lev]); - } - else if (solverChoice.windfarm_type == WindFarmType::EWP) { - ewp_advance(lev, Geom(lev), dt_lev, S_old, - U_old, V_old, W_old, vars_ewp[lev], Nturb[lev]); + if (solverChoice.windfarm_type != WindFarmType::None) { + advance_windfarm(lev, Geom(lev), dt_lev, S_old, + U_old, V_old, W_old, vars_windfarm[lev], Nturb[lev], solverChoice); } #endif diff --git a/Source/TimeIntegration/ERF_make_tau_terms.cpp b/Source/TimeIntegration/ERF_make_tau_terms.cpp index 456697719..6d3300193 100644 --- a/Source/TimeIntegration/ERF_make_tau_terms.cpp +++ b/Source/TimeIntegration/ERF_make_tau_terms.cpp @@ -10,7 +10,6 @@ using namespace amrex; void erf_make_tau_terms (int level, int nrk, - const Gpu::DeviceVector& domain_bcs_type_d, const Vector& domain_bcs_type_h, std::unique_ptr& z_phys_nd, Vector& S_data, @@ -23,7 +22,6 @@ void erf_make_tau_terms (int level, int nrk, MultiFab* Tau23, MultiFab* Tau31, MultiFab* Tau32, MultiFab* SmnSmn, MultiFab* eddyDiffs, - MultiFab* Hfx3, MultiFab* Diss, const Geometry geom, const SolverChoice& solverChoice, std::unique_ptr& most, @@ -34,25 +32,15 @@ void erf_make_tau_terms (int level, int nrk, { BL_PROFILE_REGION("erf_make_tau_terms()"); - const BCRec* bc_ptr_d = domain_bcs_type_d.data(); const BCRec* bc_ptr_h = domain_bcs_type_h.data(); DiffChoice dc = solverChoice.diffChoice; TurbChoice tc = solverChoice.turbChoice[level]; - int start_comp = 0; - int num_comp = 2; - int end_comp = start_comp + num_comp - 1; - - const AdvType l_horiz_adv_type = solverChoice.advChoice.dycore_horiz_adv_type; - const AdvType l_vert_adv_type = solverChoice.advChoice.dycore_vert_adv_type; - const Real l_horiz_upw_frac = solverChoice.advChoice.dycore_horiz_upw_frac; - const Real l_vert_upw_frac = solverChoice.advChoice.dycore_vert_upw_frac; const bool l_use_terrain = solverChoice.use_terrain; const bool l_moving_terrain = (solverChoice.terrain_type == TerrainType::Moving); if (l_moving_terrain) AMREX_ALWAYS_ASSERT (l_use_terrain); - const bool l_reflux = (solverChoice.coupling_type == CouplingType::TwoWay); const bool l_use_diff = ( (dc.molec_diff_type != MolecDiffType::None) || (tc.les_type != LESType::None) || @@ -63,22 +51,13 @@ void erf_make_tau_terms (int level, int nrk, tc.pbl_type == PBLType::MYNN25 || tc.pbl_type == PBLType::YSU ); - const bool use_moisture = (solverChoice.moisture_type != MoistureType::None); const bool use_most = (most != nullptr); const bool exp_most = (solverChoice.use_explicit_most); const Box& domain = geom.Domain(); - const int domhi_z = domain.bigEnd(2); const int domlo_z = domain.smallEnd(2); const GpuArray dxInv = geom.InvCellSizeArray(); - const Real* dx = geom.CellSize(); - - // ***************************************************************************** - // Combine external forcing terms - // ***************************************************************************** - const Array grav{0.0, 0.0, -solverChoice.gravity}; - const GpuArray grav_gpu{grav[0], grav[1], grav[2]}; // ***************************************************************************** // Pre-computed quantities diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index 22e4067c2..b985380d0 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -257,7 +257,7 @@ void erf_slow_rhs_post (int level, int finest_level, // We have projected the velocities stored in S_data but we will use // the velocities stored in S_scratch to update the scalars, so // we need to copy from S_data (projected) into S_scratch - if (solverChoice.incompressible) { + if (solverChoice.incompressible[level]) { Box tbx_inc = mfi.nodaltilebox(0); Box tby_inc = mfi.nodaltilebox(1); Box tbz_inc = mfi.nodaltilebox(2); diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 188e27dbc..6b33f0e9f 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -151,7 +151,7 @@ void erf_slow_rhs_pre (int level, int finest_level, const bool l_exp_most = (solverChoice.use_explicit_most); #ifdef ERF_USE_POISSON_SOLVE - const bool l_incompressible = solverChoice.incompressible; + const bool l_incompressible = solverChoice.incompressible[level]; const bool l_const_rho = solverChoice.constant_density; // We cannot use incompressible with terrain or with moisture @@ -164,7 +164,6 @@ void erf_slow_rhs_pre (int level, int finest_level, const Box& domain = geom.Domain(); const int domhi_z = domain.bigEnd(2); - const int domlo_z = domain.smallEnd(2); const GpuArray dxInv = geom.InvCellSizeArray(); const Real* dx = geom.CellSize(); @@ -188,10 +187,10 @@ void erf_slow_rhs_pre (int level, int finest_level, std::unique_ptr dflux_z; if (l_use_diff) { - erf_make_tau_terms(level,nrk,domain_bcs_type_d,domain_bcs_type_h,z_phys_nd, + erf_make_tau_terms(level,nrk,domain_bcs_type_h,z_phys_nd, S_data,xvel,yvel,zvel,Omega, Tau11,Tau22,Tau33,Tau12,Tau13,Tau21,Tau23,Tau31,Tau32, - SmnSmn,eddyDiffs,Hfx3,Diss,geom,solverChoice,most, + SmnSmn,eddyDiffs,geom,solverChoice,most, detJ,mapfac_m,mapfac_u,mapfac_v); dflux_x = std::make_unique(convert(ba,IntVect(1,0,0)), dm, nvars, 0); @@ -238,7 +237,6 @@ void erf_slow_rhs_pre (int level, int finest_level, const Array4& rho_u_old = S_old[IntVars::xmom].array(mfi); const Array4& rho_v_old = S_old[IntVars::ymom].array(mfi); - const Array4& rho_w_old = S_old[IntVars::zmom].array(mfi); if (l_incompressible) { // When incompressible we must reset these to 0 each RK step diff --git a/Source/TimeIntegration/TI_slow_headers.H b/Source/TimeIntegration/TI_slow_headers.H index dafbf5f1b..f3da39587 100644 --- a/Source/TimeIntegration/TI_slow_headers.H +++ b/Source/TimeIntegration/TI_slow_headers.H @@ -21,7 +21,6 @@ #endif void erf_make_tau_terms (int level, int nrk, - const amrex::Gpu::DeviceVector& domain_bcs_type_d, const amrex::Vector& domain_bcs_type, std::unique_ptr& z_phys_nd, amrex::Vector& S_data, @@ -40,8 +39,6 @@ void erf_make_tau_terms (int level, int nrk, amrex::MultiFab* Tau32, amrex::MultiFab* SmnSmn, amrex::MultiFab* eddyDiffs, - amrex::MultiFab* Hfx3, - amrex::MultiFab* Diss, const amrex::Geometry geom, const SolverChoice& solverChoice, std::unique_ptr& most, diff --git a/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp b/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp index 13c360363..fbe9f506e 100644 --- a/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp +++ b/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp @@ -7,12 +7,12 @@ Real z_c_ewp = 100.0; Real C_T_ewp = 0.8, C_TKE_ewp = 0.0; void ewp_advance (int lev, - const Geometry& geom, - const Real& dt_advance, - MultiFab& cons_in, - MultiFab& U_old, MultiFab& V_old, MultiFab& W_old, - MultiFab& mf_vars_ewp, const amrex::MultiFab& mf_Nturb) -{ + const Geometry& geom, + const Real& dt_advance, + MultiFab& cons_in, + MultiFab& U_old, MultiFab& V_old, MultiFab& W_old, + MultiFab& mf_vars_ewp, const amrex::MultiFab& mf_Nturb) + { ewp_source_terms_cellcentered(geom, cons_in, U_old, V_old, W_old, mf_vars_ewp, mf_Nturb); ewp_update(dt_advance, cons_in, U_old, V_old, mf_vars_ewp); } diff --git a/Source/WindFarmParametrization/Make.package b/Source/WindFarmParametrization/Make.package new file mode 100644 index 000000000..1918ac501 --- /dev/null +++ b/Source/WindFarmParametrization/Make.package @@ -0,0 +1,2 @@ +CEXE_sources += WindFarm.cpp +CEXE_headers += WindFarm.H diff --git a/Source/WindFarmParametrization/WindFarm.H b/Source/WindFarmParametrization/WindFarm.H new file mode 100644 index 000000000..132e8706c --- /dev/null +++ b/Source/WindFarmParametrization/WindFarm.H @@ -0,0 +1,16 @@ +#ifndef WINDFARM_H +#define WINDFARM_H + +#include +#include +#include + +void advance_windfarm (int lev, + const amrex::Geometry& geom, + const amrex::Real& dt_advance, + amrex::MultiFab& cons_in, + amrex::MultiFab& U_old, amrex::MultiFab& V_old, amrex::MultiFab& W_old, + amrex::MultiFab& mf_vars_windfarm, const amrex::MultiFab& mf_Nturb, + SolverChoice& solver_choice); +#endif + diff --git a/Source/WindFarmParametrization/WindFarm.cpp b/Source/WindFarmParametrization/WindFarm.cpp new file mode 100644 index 000000000..4162e7b73 --- /dev/null +++ b/Source/WindFarmParametrization/WindFarm.cpp @@ -0,0 +1,22 @@ +#include +#include +#include +using namespace amrex; + + +void advance_windfarm (int lev, + const Geometry& geom, + const Real& dt_advance, + MultiFab& cons_in, + MultiFab& U_old, MultiFab& V_old, MultiFab& W_old, + MultiFab& mf_vars_windfarm, const amrex::MultiFab& mf_Nturb, + SolverChoice& solver_choice) +{ + if (solver_choice.windfarm_type == WindFarmType::Fitch) { + fitch_advance(lev, geom, dt_advance, cons_in, U_old, V_old, W_old, + mf_vars_windfarm, mf_Nturb); + } else if (solver_choice.windfarm_type == WindFarmType::EWP) { + ewp_advance(lev, geom, dt_advance, cons_in, U_old, V_old, W_old, + mf_vars_windfarm, mf_Nturb); + } +} From 9b33ecc981f2006935d1ffb27692714a32bd559e Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Wed, 15 May 2024 10:34:17 -0700 Subject: [PATCH 11/20] fix incompressible build (#1622) --- Source/TimeIntegration/ERF_advance_dycore.cpp | 2 +- Source/TimeIntegration/TI_no_substep_fun.H | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index ae74c83d9..d4691bf38 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -272,7 +272,7 @@ void ERF::advance_dycore(int level, mri_integrator.set_post_update(post_update_fun); #ifdef ERF_USE_POISSON_SOLVE - if (solverChoice.incompressible) { + if (solverChoice.incompressible[level]) { mri_integrator.set_slow_rhs_inc(slow_rhs_fun_inc); } #endif diff --git a/Source/TimeIntegration/TI_no_substep_fun.H b/Source/TimeIntegration/TI_no_substep_fun.H index bebba6845..57b18d878 100644 --- a/Source/TimeIntegration/TI_no_substep_fun.H +++ b/Source/TimeIntegration/TI_no_substep_fun.H @@ -135,7 +135,7 @@ apply_bcs(S_sum, time_for_fp, ng_cons, ng_vel, fast_only=true, vel_and_mom_synced=false); #ifdef ERF_USE_POISSON_SOLVE - if (solverChoice.incompressible) { + if (solverChoice.incompressible[level]) { bool have_tb = (thin_xforce[0] || thin_yforce[0] || thin_zforce[0]); if (solverChoice.project_every_stage || (nrk==1)) { if (have_tb) { From 80f78f69f35ffe3db35c98507b4cf3dfb9149e66 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Thu, 16 May 2024 12:53:44 -0700 Subject: [PATCH 12/20] Fix custom srcs for bomex. (#1624) --- Exec/DevTests/Bomex/prob.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Exec/DevTests/Bomex/prob.cpp b/Exec/DevTests/Bomex/prob.cpp index 9470eb58d..4090dba3f 100644 --- a/Exec/DevTests/Bomex/prob.cpp +++ b/Exec/DevTests/Bomex/prob.cpp @@ -227,7 +227,8 @@ Problem::update_rhotheta_sources (const Real& /*time*/, if (z_cc < parms.cutoff) { src[k] = parms.advection_heating_rate; } else if (z_cc < parms.cutoff+parms.cutoff_transition) { - src[k] = parms.advection_heating_rate * (z_cc-parms.cutoff)/parms.cutoff_transition; + Real slope = -parms.advection_heating_rate / parms.cutoff_transition; + src[k] = (z_cc-parms.cutoff) * slope + parms.advection_heating_rate; } else { src[k] = 0.0; } @@ -269,7 +270,8 @@ Problem::update_rhoqt_sources (const Real& /*time*/, if (z_cc < parms.moisture_cutoff) { qsrc[k] = parms.advection_moisture_rate; } else if (z_cc < parms.moisture_cutoff+parms.moisture_cutoff_transition) { - qsrc[k] = parms.advection_moisture_rate * (z_cc-parms.moisture_cutoff)/parms.moisture_cutoff_transition; + Real slope = -parms.advection_moisture_rate / parms.moisture_cutoff_transition; + qsrc[k] = (z_cc-parms.moisture_cutoff) * slope + parms.advection_moisture_rate; } else { qsrc[k] = 0.0; } From 0c7d8955ede15eae3d335fbacdfe7cf69898f975 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Wed, 22 May 2024 15:45:23 -0700 Subject: [PATCH 13/20] Bdy Plane Offset Fix (#1626) * Need to use the delta distance from the problem low to comput i,j,k. * Forgot const. --- Source/IO/ERF_WriteBndryPlanes.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/Source/IO/ERF_WriteBndryPlanes.cpp b/Source/IO/ERF_WriteBndryPlanes.cpp index fc8d493e0..eda3353c9 100644 --- a/Source/IO/ERF_WriteBndryPlanes.cpp +++ b/Source/IO/ERF_WriteBndryPlanes.cpp @@ -64,15 +64,16 @@ WriteBndryPlanes::WriteBndryPlanes (Vector& grids, // If the target area is contained at a finer level, use the finest data possible for (int ilev = 0; ilev < grids.size(); ilev++) { - auto const dxi = geom[ilev].InvCellSizeArray(); + const Real* xLo = m_geom[ilev].ProbLo(); + auto const dxi = geom[ilev].InvCellSizeArray(); const Box& domain = m_geom[ilev].Domain(); // We create the smallest box that contains all of the cell centers // in the physical region specified - int ilo = static_cast(Math::floor(box_lo[0] * dxi[0])+.5); - int jlo = static_cast(Math::floor(box_lo[1] * dxi[1])+.5); - int ihi = static_cast(Math::floor(box_hi[0] * dxi[0])+.5)-1; - int jhi = static_cast(Math::floor(box_hi[1] * dxi[1])+.5)-1; + int ilo = static_cast(Math::floor((box_lo[0] - xLo[0]) * dxi[0])+.5); + int jlo = static_cast(Math::floor((box_lo[1] - xLo[1]) * dxi[1])+.5); + int ihi = static_cast(Math::floor((box_hi[0] - xLo[0]) * dxi[0])+.5)-1; + int jhi = static_cast(Math::floor((box_hi[1] - xLo[1]) * dxi[1])+.5)-1; // Map this to index space -- for now we do no interpolation target_box.setSmall(IntVect(ilo,jlo,0)); From 517e56c77eee49c3715f3730a4d0d08f35c1bb57 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Fri, 24 May 2024 11:53:25 -0700 Subject: [PATCH 14/20] SAM model without precipitation and ice (#1628) * fixes to thermal sourcing. * Adding SAM_NoPrecip_NoIce * Removing unwanted print * Correcting GPU errors * Correcting GPU errors * Checking GPU compilation * Checking GPU compilation * Correcting GPU compilation --------- Co-authored-by: Aaron Lattanzi Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> Co-authored-by: Mahesh Natarajan Co-authored-by: Ann Almgren --- ...nputs_BF02_moist_bubble_SAM_NoPrecip_NoIce | 81 +++++++++++ Exec/RegTests/Bubble/prob.cpp | 15 +- Source/DataStructs/DataStruct.H | 7 +- Source/Microphysics/EulerianMicrophysics.H | 3 +- Source/Microphysics/Microphysics.H | 1 + Source/Microphysics/SAM/Cloud_SAM.cpp | 137 +++++++++++------- Source/Microphysics/SAM/Init_SAM.cpp | 2 +- Source/Microphysics/SAM/Precip.cpp | 16 +- Source/Microphysics/SAM/SAM.H | 14 +- Source/SourceTerms/ERF_make_buoyancy.cpp | 2 +- 10 files changed, 204 insertions(+), 74 deletions(-) create mode 100644 Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoPrecip_NoIce diff --git a/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoPrecip_NoIce b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoPrecip_NoIce new file mode 100644 index 000000000..fc18f8f23 --- /dev/null +++ b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoPrecip_NoIce @@ -0,0 +1,81 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 2000 +stop_time = 3600.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 20000.0 400.0 10000.0 +amr.n_cell = 200 4 100 +geometry.is_periodic = 0 1 0 +xlo.type = "SlipWall" +xhi.type = "SlipWall" +zlo.type = "SlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 0.5 +erf.fixed_mri_dt_ratio = 4 +#erf.no_substepping = 1 +#erf.fixed_dt = 0.1 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 100 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoadv_0 x_velocity y_velocity z_velocity pressure theta scalar temp pres_hse dens_hse pert_pres pert_dens eq_pot_temp qt qv qc qi qrain qsnow qgraup + +# SOLVER CHOICES +erf.use_gravity = true +erf.use_coriolis = false + +erf.dycore_horiz_adv_type = "Upwind_3rd" +erf.dycore_vert_adv_type = "Upwind_3rd" +erf.dryscal_horiz_adv_type = "Upwind_3rd" +erf.dryscal_vert_adv_type = "Upwind_3rd" +erf.moistscal_horiz_adv_type = "Upwind_3rd" +erf.moistscal_vert_adv_type = "Upwind_3rd" + +# PHYSICS OPTIONS +erf.les_type = "None" +erf.pbl_type = "None" +erf.moisture_model = "SAM_NoPrecip_NoIce" +erf.buoyancy_type = 1 +erf.use_moist_background = true + +erf.molec_diff_type = "ConstantAlpha" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 0.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 0.0 # [m^2/s] +erf.alpha_C = 0.0 + +# INITIAL CONDITIONS +#erf.init_type = "input_sounding" +#erf.input_sounding_file = "BF02_moist_sounding" +#erf.init_sounding_ideal = true + +# PROBLEM PARAMETERS (optional) +# warm bubble input +prob.x_c = 10000.0 +prob.z_c = 2000.0 +prob.x_r = 2000.0 +prob.z_r = 2000.0 +prob.T_0 = 300.0 + +prob.do_moist_bubble = true +prob.theta_pert = 2.0 +prob.qt_init = 0.02 +prob.eq_pot_temp = 320.0 diff --git a/Exec/RegTests/Bubble/prob.cpp b/Exec/RegTests/Bubble/prob.cpp index 8957ab3f1..b6382081b 100644 --- a/Exec/RegTests/Bubble/prob.cpp +++ b/Exec/RegTests/Bubble/prob.cpp @@ -360,6 +360,14 @@ Problem::init_custom_pert( const Real x_r = parms.x_r, y_r = parms.y_r, z_r = parms.z_r; const Real theta_pert = parms.theta_pert; + int moisture_type = 1; + + if(sc.moisture_type == MoistureType::SAM) { + moisture_type = 1; + } else if(sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) { + moisture_type = 2; + } + amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) { // Geometry (note we must include these here to get the data on device) @@ -410,7 +418,12 @@ Problem::init_custom_pert( // Cold microphysics are present int nstate = state_pert.ncomp; if (nstate == NVAR_max) { - Real omn = std::max(0.0,std::min(1.0,(T-tbgmin)*a_bg)); + Real omn; + if(moisture_type == 1) { + omn = std::max(0.0,std::min(1.0,(T-tbgmin)*a_bg)); + } else if(moisture_type == 2) { + omn = 1.0; + } Real qn = state_pert(i, j, k, RhoQ2_comp); state_pert(i, j, k, RhoQ2_comp) = qn * omn; state_pert(i, j, k, RhoQ3_comp) = qn * (1.0-omn); diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index b738ec678..0ab5ea5cc 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -33,7 +33,7 @@ enum struct MoistureModelType{ }; enum struct MoistureType { - Kessler, SAM, Kessler_NoRain, None + Kessler, SAM, SAM_NoPrecip_NoIce, Kessler_NoRain, None }; enum struct WindFarmType { @@ -95,6 +95,8 @@ struct SolverChoice { pp.query("moisture_model", moisture_model_string); if (moisture_model_string == "SAM") { moisture_type = MoistureType::SAM; + } else if (moisture_model_string == "SAM_NoPrecip_NoIce") { + moisture_type = MoistureType::SAM_NoPrecip_NoIce; } else if (moisture_model_string == "Kessler") { moisture_type = MoistureType::Kessler; }else if (moisture_model_string == "Kessler_NoRain") { @@ -107,7 +109,8 @@ struct SolverChoice { // Set a different default for moist vs dry if (moisture_type != MoistureType::None) { if (moisture_type == MoistureType::Kessler_NoRain || - moisture_type == MoistureType::SAM) { + moisture_type == MoistureType::SAM || + moisture_type == MoistureType::SAM_NoPrecip_NoIce) { buoyancy_type = 1; // asserted in make buoyancy } else { buoyancy_type = 2; // uses Tprime diff --git a/Source/Microphysics/EulerianMicrophysics.H b/Source/Microphysics/EulerianMicrophysics.H index 4df1af178..4e1ac1af2 100644 --- a/Source/Microphysics/EulerianMicrophysics.H +++ b/Source/Microphysics/EulerianMicrophysics.H @@ -27,7 +27,8 @@ public: { AMREX_ASSERT( Microphysics::modelType(a_model_type) == MoistureModelType::Eulerian ); m_moist_model.resize(nlev); - if (a_model_type == MoistureType::SAM) { + if (a_model_type == MoistureType::SAM or + a_model_type == MoistureType::SAM_NoPrecip_NoIce) { SetModel(); amrex::Print() << "SAM moisture model!\n"; } else if (a_model_type == MoistureType::Kessler or diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index efd601460..ea96f0e79 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -58,6 +58,7 @@ public: static MoistureModelType modelType (const MoistureType a_moisture_type) { if ( (a_moisture_type == MoistureType::SAM) + || (a_moisture_type == MoistureType::SAM_NoPrecip_NoIce) || (a_moisture_type == MoistureType::Kessler) || (a_moisture_type == MoistureType::Kessler_NoRain) || (a_moisture_type == MoistureType::None) ) { diff --git a/Source/Microphysics/SAM/Cloud_SAM.cpp b/Source/Microphysics/SAM/Cloud_SAM.cpp index 8c9264f3e..89483d887 100644 --- a/Source/Microphysics/SAM/Cloud_SAM.cpp +++ b/Source/Microphysics/SAM/Cloud_SAM.cpp @@ -8,7 +8,7 @@ using namespace amrex; /** * Split cloud components according to saturation pressures; source theta from latent heat. */ -void SAM::Cloud () { +void SAM::Cloud (const SolverChoice& solverChoice) { constexpr Real an = 1.0/(tbgmax-tbgmin); constexpr Real bn = tbgmin*an; @@ -20,6 +20,16 @@ void SAM::Cloud () { Real tol = 1.0e-4; + SolverChoice sc = solverChoice; + + int moisture_type = 1; + + if(solverChoice.moisture_type == MoistureType::SAM) { + moisture_type = 1; + } else if(solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) { + moisture_type = 2; + } + for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]), TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi); auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi); @@ -54,42 +64,57 @@ void SAM::Cloud () { // This ensures the omn splitting is enforced // before the Newton iteration, which assumes it is. - // Cloud ice not permitted (melt to form water) - if(tabs_array(i,j,k) >= tbgmax) { - omn = 1.0; - delta_qi = qci_array(i,j,k); - qci_array(i,j,k) = 0.0; - qcl_array(i,j,k) += delta_qi; - tabs_array(i,j,k) -= fac_fus * delta_qi; - pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - * (1.0 + R_v/R_d * qv_array(i,j,k)); - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); - pres_array(i,j,k) /= 100.0; + omn = 1.0; + if(moisture_type == 1){ + // Cloud ice not permitted (melt to form water) + if(tabs_array(i,j,k) >= tbgmax) { + omn = 1.0; + delta_qi = qci_array(i,j,k); + qci_array(i,j,k) = 0.0; + qcl_array(i,j,k) += delta_qi; + tabs_array(i,j,k) -= fac_fus * delta_qi; + pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) + * (1.0 + R_v/R_d * qv_array(i,j,k)); + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); + pres_array(i,j,k) *= 0.01; + } + // Cloud water not permitted (freeze to form ice) + else if(tabs_array(i,j,k) <= tbgmin) { + omn = 0.0; + delta_qc = qcl_array(i,j,k); + qcl_array(i,j,k) = 0.0; + qci_array(i,j,k) += delta_qc; + tabs_array(i,j,k) += fac_fus * delta_qc; + pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) + * (1.0 + R_v/R_d * qv_array(i,j,k)); + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); + pres_array(i,j,k) *= 0.01; + } + // Mixed cloud phase (split according to omn) + else { + omn = an*tabs_array(i,j,k)-bn; + delta_qc = qcl_array(i,j,k) - qn_array(i,j,k) * omn; + delta_qi = qci_array(i,j,k) - qn_array(i,j,k) * (1.0 - omn); + qcl_array(i,j,k) = qn_array(i,j,k) * omn; + qci_array(i,j,k) = qn_array(i,j,k) * (1.0 - omn); + tabs_array(i,j,k) += fac_fus * delta_qc; + pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) + * (1.0 + R_v/R_d * qv_array(i,j,k)); + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); + pres_array(i,j,k) *= 0.01; + } } - // Cloud water not permitted (freeze to form ice) - else if(tabs_array(i,j,k) <= tbgmin) { - omn = 0.0; - delta_qc = qcl_array(i,j,k); - qcl_array(i,j,k) = 0.0; - qci_array(i,j,k) += delta_qc; - tabs_array(i,j,k) += fac_fus * delta_qc; - pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - * (1.0 + R_v/R_d * qv_array(i,j,k)); - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); - pres_array(i,j,k) /= 100.0; - } - // Mixed cloud phase (split according to omn) - else { - omn = an*tabs_array(i,j,k)-bn; - delta_qc = qcl_array(i,j,k) - qn_array(i,j,k) * omn; - delta_qi = qci_array(i,j,k) - qn_array(i,j,k) * (1.0 - omn); - qcl_array(i,j,k) = qn_array(i,j,k) * omn; - qci_array(i,j,k) = qn_array(i,j,k) * (1.0 - omn); - tabs_array(i,j,k) += fac_fus * delta_qc; + else if(moisture_type == 2){ + // No ice. ie omn = 1.0 + delta_qc = qcl_array(i,j,k) - qn_array(i,j,k); + delta_qi = 0.0; + qcl_array(i,j,k) = qn_array(i,j,k); + qci_array(i,j,k) = 0.0; + tabs_array(i,j,k) += fac_cond * delta_qc; pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) * (1.0 + R_v/R_d * qv_array(i,j,k)); theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); - pres_array(i,j,k) /= 100.0; + pres_array(i,j,k) *= 0.01; } // Initial guess for temperature & pressure @@ -120,19 +145,24 @@ void SAM::Cloud () { erf_dtqsatw(tabs, pres, dqsatw); erf_dtqsati(tabs, pres, dqsati); - // Cloud ice not permitted (condensation & fusion) - if(tabs >= tbgmax) { - omn = 1.0; - } - // Cloud water not permitted (sublimation & fusion) - else if(tabs <= tbgmin) { - omn = 0.0; - lstarw = fac_sub; - } - // Mixed cloud phase (condensation & fusion) - else { - omn = an*tabs-bn; - domn = an; + if(moisture_type == 1) { + // Cloud ice not permitted (condensation & fusion) + if(tabs >= tbgmax) { + omn = 1.0; + } + // Cloud water not permitted (sublimation & fusion) + else if(tabs <= tbgmin) { + omn = 0.0; + lstarw = fac_sub; + } + // Mixed cloud phase (condensation & fusion) + else { + omn = an*tabs-bn; + domn = an; + } + } else if(moisture_type == 2) { + omn = 1.0; + domn = 0.0; } // Linear combination of each component @@ -155,8 +185,7 @@ void SAM::Cloud () { // Update the pressure pres = rho_array(i,j,k) * R_d * tabs - * (1.0 + R_v/R_d * qsat); - pres /= 100.0; + * (1.0 + R_v/R_d * qsat) * 0.01; // Update iteration niter = niter+1; @@ -188,15 +217,15 @@ void SAM::Cloud () { theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); // Pressure unit conversion - pres_array(i,j,k) /= 100.0; + pres_array(i,j,k) *= 0.01; } // We cannot blindly relax to qsat, but we can convert qc/qi -> qv else { // Changes in each component - delta_qv = qcl_array(i,j,k) + qci_array(i,j,k); - delta_qc = -qcl_array(i,j,k); - delta_qi = -qci_array(i,j,k); + delta_qv = qcl_array(i,j,k) + qci_array(i,j,k); + delta_qc = qcl_array(i,j,k); + delta_qi = qci_array(i,j,k); // Partition the change in non-precipitating q qv_array(i,j,k) += delta_qv; @@ -206,7 +235,7 @@ void SAM::Cloud () { qt_array(i,j,k) = qv_array(i,j,k); // Update temperature (endothermic since we evap/sublime) - tabs_array(i,j,k) -= fac_fus * delta_qc + fac_sub * delta_qi; + tabs_array(i,j,k) -= fac_cond * delta_qc + fac_sub * delta_qi; // Update pressure pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) @@ -216,7 +245,7 @@ void SAM::Cloud () { theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); // Pressure unit conversion - pres_array(i,j,k) /= 100.0; + pres_array(i,j,k) *= 0.01; } }); } // mfi diff --git a/Source/Microphysics/SAM/Init_SAM.cpp b/Source/Microphysics/SAM/Init_SAM.cpp index 0fc219432..1d593044f 100644 --- a/Source/Microphysics/SAM/Init_SAM.cpp +++ b/Source/Microphysics/SAM/Init_SAM.cpp @@ -129,7 +129,7 @@ void SAM::Copy_State_to_Micro (const MultiFab& cons_in) tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp), states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); - pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k))/100.; + pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)) * 0.01; }); } } diff --git a/Source/Microphysics/SAM/Precip.cpp b/Source/Microphysics/SAM/Precip.cpp index e28aadfae..296cc776c 100644 --- a/Source/Microphysics/SAM/Precip.cpp +++ b/Source/Microphysics/SAM/Precip.cpp @@ -167,12 +167,12 @@ void SAM::Precip () { qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k); - // Latent heat source for theta + // Latent heat source for theta (endothermic fusion & exothermic melting) tabs_array(i,j,k) += fac_fus * ( dqca * (1.0 - omp) - dqia * omp ); pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) * (1.0 + R_v/R_d * qv_array(i,j,k)); theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); - pres_array(i,j,k) /= 100.0; + pres_array(i,j,k) *= 0.01; } //================================================== @@ -188,12 +188,12 @@ void SAM::Precip () { dqpg = evapg1_t(k)*sqrt(qpg) + evapg2_t(k)*pow(qpg,powg2); // NOTE: This is always a sink for precipitating comps - // since qv0. If we are // in a super-saturated state (qv>qsat) the Newton // iterations in Cloud() will have handled condensation. - dqpr *= -dtn * (qv_array(i,j,k)/qsat - 1.0); - dqps *= -dtn * (qv_array(i,j,k)/qsat - 1.0); - dqpg *= -dtn * (qv_array(i,j,k)/qsat - 1.0); + dqpr *= dtn * (1.0 - qv_array(i,j,k)/qsat); + dqps *= dtn * (1.0 - qv_array(i,j,k)/qsat); + dqpg *= dtn * (1.0 - qv_array(i,j,k)/qsat); // Limit to avoid negative moisture fractions dqpr = std::min(qpr_array(i,j,k),dqpr); @@ -211,8 +211,8 @@ void SAM::Precip () { qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k); - // Latent heat source for theta - tabs_array(i,j,k) += fac_cond * dqpr + fac_sub * (dqps + dqpg); + // Latent heat source for theta (endothermic) + tabs_array(i,j,k) -= fac_cond * dqpr + fac_sub * (dqps + dqpg); pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) * (1.0 + R_v/R_d * qv_array(i,j,k)); theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index 998712d89..3b4d9f280 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -63,7 +63,7 @@ public: virtual ~SAM () = default; // cloud physics - void Cloud (); + void Cloud (const SolverChoice& sc); // ice physics void IceFall (); @@ -121,14 +121,16 @@ public: // wrapper to do all the updating void Advance (const amrex::Real& dt_advance, - const SolverChoice& /*solverChoice*/) override + const SolverChoice& sc) override { dt = dt_advance; - this->Cloud(); - this->IceFall(); - this->Precip(); - this->PrecipFall(); + this->Cloud(sc); + if(sc.moisture_type == MoistureType::SAM) { + this->IceFall(); + this->Precip(); + this->PrecipFall(); + } } amrex::MultiFab* diff --git a/Source/SourceTerms/ERF_make_buoyancy.cpp b/Source/SourceTerms/ERF_make_buoyancy.cpp index 5dabd32dd..65533381e 100644 --- a/Source/SourceTerms/ERF_make_buoyancy.cpp +++ b/Source/SourceTerms/ERF_make_buoyancy.cpp @@ -136,7 +136,7 @@ void make_buoyancy (Vector& S_data, AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); } - if (solverChoice.moisture_type == MoistureType::SAM) { + if (solverChoice.moisture_type == MoistureType::SAM or solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) { AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); } From 120f71f19bcfc617f83923d8413c2e581ed59303 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Fri, 24 May 2024 15:59:20 -0700 Subject: [PATCH 15/20] Fix bomex MOST. (#1627) Co-authored-by: Aaron Lattanzi Co-authored-by: Ann Almgren --- Source/BoundaryConditions/MOSTStress.H | 9 +++------ Source/Diffusion/ComputeTurbulentViscosity.cpp | 5 +++-- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index 8d39e12b8..8cc6f27cc 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -1198,7 +1198,7 @@ struct custom_flux je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s] eta = amrex::max(eta,eta_eps); - dest_arr(i,j,k,icomp+n) = rho*(qv - moflux*rho/eta*deltaz); + dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) - moflux*rho/eta*deltaz; return moflux; } @@ -1238,7 +1238,7 @@ struct custom_flux amrex::Real ustar = u_star_arr(ic,jc,zlo); amrex::Real tstar = t_star_arr(ic,jc,zlo); - amrex::Real moflux = (std::abs(tstar) > eps) ? -tstar*ustar : 0.0; + amrex::Real moflux = (std::abs(tstar) > eps) ? -rho * tstar : 0.0; amrex::Real deltaz = dz * (zlo - k); if (exp_most) { @@ -1254,10 +1254,7 @@ struct custom_flux je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] eta = amrex::max(eta,eta_eps); - // Note: Kh = eta/rho - // hfx = -Kh dT/dz ==> +ve hfx corresponds to heating from the surface - // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel - dest_arr(i,j,k,icomp+n) = rho*(theta + moflux*rho/eta*deltaz); + dest_arr(i,j,k,icomp+n) = dest_arr(i,j,zlo,icomp+n) - moflux*rho/eta*deltaz; } return moflux; diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index 44c90a505..db7742162 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -191,10 +191,11 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, // Calculate SFS quantities // - dissipation Real Ce; - if ((l_C_e_wall > 0) && (k==0)) + if ((l_C_e_wall > 0) && (k==0)) { Ce = l_C_e_wall; - else + } else { Ce = 1.9*l_C_k + Ce_lcoeff*length / DeltaMsf; + } diss(i,j,k) = cell_data(i,j,k,Rho_comp) * Ce * std::pow(E,1.5) / length; // - heat flux // (Note: If using ERF_EXPLICIT_MOST_STRESS, the value at k=0 will From fb65da652b890ebc4deaa4487c4a7c098801353c Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 24 May 2024 18:12:57 -0700 Subject: [PATCH 16/20] fix SAM for moist bubble case, add qsat as a derived quantity and print buoyancy_type at the start (#1629) --- Source/DataStructs/DataStruct.H | 2 + Source/ERF.H | 7 +++- Source/IO/Plotfile.cpp | 48 ++++++++++++++++++++-- Source/Microphysics/Kessler/Kessler.cpp | 43 +++++++++----------- Source/Microphysics/SAM/Cloud_SAM.cpp | 54 ++++++++++++++++--------- 5 files changed, 105 insertions(+), 49 deletions(-) diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 0ab5ea5cc..69a5c1b02 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -350,6 +350,8 @@ struct SolverChoice { amrex::Print() << ")" << std::endl; } + amrex::Print() << "Buoyancy_type : " << buoyancy_type << std::endl; + advChoice.display(); diffChoice.display(); spongeChoice.display(); diff --git a/Source/ERF.H b/Source/ERF.H index e299cb61c..d3712f66d 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -799,7 +799,8 @@ private: const amrex::Vector derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar", "vorticity_x","vorticity_y","vorticity_z", "magvel", "divU", - "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "eq_pot_temp", "num_turb", + "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", + "eq_pot_temp", "num_turb", "dpdx", "dpdy", "pres_hse_x", "pres_hse_y", "z_phys", "detJ" , "mapfac", "lat_m", "lon_m", // Time averaged velocity @@ -811,8 +812,10 @@ private: // mynn pbl lengthscale "Lpbl", // moisture vars - "qt", "qv", "qc", "qi", "qp", "qrain", "qsnow", "qgraup", "rain_accum", "snow_accum", "graup_accum" + "qt", "qv", "qc", "qi", "qp", "qrain", "qsnow", "qgraup", "qsat", + "rain_accum", "snow_accum", "graup_accum" #ifdef ERF_COMPUTE_ERROR + // error vars ,"xvel_err", "yvel_err", "zvel_err", "pp_err" #endif }; diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 1ca5db770..ada718f94 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -51,9 +51,15 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector for (int i = 0; i < ncomp_cons; ++i) { if ( containerHasElement(plot_var_names, cons_names[i]) ) { - tmp_plot_names.push_back(cons_names[i]); + if ( (solverChoice.moisture_type == MoistureType::SAM) || (derived_names[i] != "rhoQ4" && + derived_names[i] != "rhoQ5" && + derived_names[i] != "rhoQ6") ) + { + tmp_plot_names.push_back(cons_names[i]); + } // moisture_type } } + // check for velocity since it's not in cons_names // if we are asked for any velocity component, we will need them all if (containerHasElement(plot_var_names, "x_velocity") || @@ -63,13 +69,26 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector tmp_plot_names.push_back("y_velocity"); tmp_plot_names.push_back("z_velocity"); } + + // + // If the model we are running doesn't have the variable listed in the inputs file, + // just ignore it rather than aborting + // for (int i = 0; i < derived_names.size(); ++i) { if ( containerHasElement(plot_var_names, derived_names[i]) ) { if (solverChoice.use_terrain || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) { - tmp_plot_names.push_back(derived_names[i]); - } - } + if ( (solverChoice.moisture_type == MoistureType::SAM) || (derived_names[i] != "qi" && + derived_names[i] != "qsnow" && + derived_names[i] != "qgraup" && + derived_names[i] != "snow_accum" && + derived_names[i] != "graup_accum") ) + { + tmp_plot_names.push_back(derived_names[i]); + } // moisture_type + } // use_terrain? + } // hasElement } + #ifdef ERF_USE_PARTICLES const auto& particles_namelist( particleData.getNamesUnalloc() ); for (auto it = particles_namelist.cbegin(); it != particles_namelist.cend(); ++it) { @@ -967,6 +986,27 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp += 1; } + if (containerHasElement(plot_var_names, "qsat")) + { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Real qv = S_arr(i,j,k,RhoQ1_comp) / S_arr(i,j,k,Rho_comp); + Real T = getTgivenRandRTh(S_arr(i,j,k,Rho_comp), S_arr(i,j,k,RhoTheta_comp), qv); + Real pressure = getPgivenRTh(S_arr(i,j,k,RhoTheta_comp), qv); + erf_qsatw(T, pressure, derdat(i,j,k,mf_comp)); + }); + } + mf_comp ++; + } + if(solverChoice.moisture_type == MoistureType::Kessler){ if (containerHasElement(plot_var_names, "rain_accum")) { diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index 15d8fa5a1..7e8abbf96 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -69,7 +69,6 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) }); } - for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi); auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi); @@ -118,30 +117,28 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) //Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k)); // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature - if(qv_array(i,j,k) > qsat){ + if (qv_array(i,j,k) > qsat) { dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac)); } - - // If water vapor is less than the saturated value, then the cloud water can evaporate, leading to evaporative cooling and - // reducing temperature - if(qv_array(i,j,k) < qsat && qc_array(i,j,k) > 0.0){ + // If water vapor is less than the saturated value, then the cloud water can evaporate, + // leading to evaporative cooling and reducing temperature + if (qv_array(i,j,k) < qsat && qc_array(i,j,k) > 0.0) { dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac)); - } + } - if(qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) { + if (qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) { Real C = 1.6 + 124.9*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.2046); dq_rain_to_vapor = 1.0/(0.001*rho_array(i,j,k))*(1.0 - qv_array(i,j,k)/qsat)*C*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.525)/ (5.4e5 + 2.55e6/(pressure*qsat))*dtn; // The negative sign is to make this variable (vapor formed from evaporation) - // a poistive quantity (as qv/qs < 1) + // a positive quantity (as qv/qs < 1) dq_rain_to_vapor = std::min({qp_array(i,j,k), dq_rain_to_vapor}); // Removing latent heat due to evaporation from rain water to water vapor, reduces the (potential) temperature } // If there is cloud water present then do accretion and autoconversion to rain - if (qc_array(i,j,k) > 0.0) { qcc = qc_array(i,j,k); @@ -163,12 +160,12 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) Real dq_sed = dtn * dJinv * (1.0/rho_array(i,j,k)) * (fz_array(i,j,k+1) - fz_array(i,j,k))/dz; if(std::fabs(dq_sed) < 1e-14) dq_sed = 0.0; - qv_array(i,j,k) = qv_array(i,j,k) - dq_vapor_to_clwater + dq_clwater_to_vapor + dq_rain_to_vapor; - qc_array(i,j,k) = qc_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain; - qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor; + qv_array(i,j,k) += -dq_vapor_to_clwater + dq_clwater_to_vapor + dq_rain_to_vapor; + qc_array(i,j,k) += dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain; + qp_array(i,j,k) += dq_sed + dq_clwater_to_rain - dq_rain_to_vapor; - theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond * - (dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor); + Real theta_over_T = theta_array(i,j,k)/tabs_array(i,j,k); + theta_array(i,j,k) += theta_over_T * d_fac_cond * (dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor); qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k)); qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); @@ -179,8 +176,7 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) } } - - if(solverChoice.moisture_type == MoistureType::Kessler_NoRain){ + if (solverChoice.moisture_type == MoistureType::Kessler_NoRain){ // get the temperature, dentisy, theta, qt and qc from input for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -217,20 +213,21 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k)); // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature - if(qv_array(i,j,k) > qsat){ + if (qv_array(i,j,k) > qsat){ dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac)); } // If water vapor is less than the saturated value, then the cloud water can evaporate, leading to evaporative cooling and // reducing temperature - if(qv_array(i,j,k) < qsat && qc_array(i,j,k) > 0.0){ + if (qv_array(i,j,k) < qsat && qc_array(i,j,k) > 0.0){ dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac)); } - qv_array(i,j,k) = qv_array(i,j,k) - dq_vapor_to_clwater + dq_clwater_to_vapor; - qc_array(i,j,k) = qc_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor; + qv_array(i,j,k) += -dq_vapor_to_clwater + dq_clwater_to_vapor; + qc_array(i,j,k) += dq_vapor_to_clwater - dq_clwater_to_vapor; + + Real theta_over_T = theta_array(i,j,k)/tabs_array(i,j,k); - theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond * - (dq_vapor_to_clwater - dq_clwater_to_vapor); + theta_array(i,j,k) += theta_over_T * d_fac_cond * (dq_vapor_to_clwater - dq_clwater_to_vapor); qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k)); qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); diff --git a/Source/Microphysics/SAM/Cloud_SAM.cpp b/Source/Microphysics/SAM/Cloud_SAM.cpp index 89483d887..7783cb5bf 100644 --- a/Source/Microphysics/SAM/Cloud_SAM.cpp +++ b/Source/Microphysics/SAM/Cloud_SAM.cpp @@ -22,12 +22,12 @@ void SAM::Cloud (const SolverChoice& solverChoice) { SolverChoice sc = solverChoice; - int moisture_type = 1; + int SAM_moisture_type = 1; if(solverChoice.moisture_type == MoistureType::SAM) { - moisture_type = 1; + SAM_moisture_type = 1; } else if(solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) { - moisture_type = 2; + SAM_moisture_type = 2; } for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]), TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -65,9 +65,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) { // before the Newton iteration, which assumes it is. omn = 1.0; - if(moisture_type == 1){ + if (SAM_moisture_type == 1){ // Cloud ice not permitted (melt to form water) - if(tabs_array(i,j,k) >= tbgmax) { + if (tabs_array(i,j,k) >= tbgmax) { omn = 1.0; delta_qi = qci_array(i,j,k); qci_array(i,j,k) = 0.0; @@ -79,7 +79,7 @@ void SAM::Cloud (const SolverChoice& solverChoice) { pres_array(i,j,k) *= 0.01; } // Cloud water not permitted (freeze to form ice) - else if(tabs_array(i,j,k) <= tbgmin) { + else if (tabs_array(i,j,k) <= tbgmin) { omn = 0.0; delta_qc = qcl_array(i,j,k); qcl_array(i,j,k) = 0.0; @@ -104,7 +104,8 @@ void SAM::Cloud (const SolverChoice& solverChoice) { pres_array(i,j,k) *= 0.01; } } - else if(moisture_type == 2){ + else if (SAM_moisture_type == 2) + { // No ice. ie omn = 1.0 delta_qc = qcl_array(i,j,k) - qn_array(i,j,k); delta_qi = 0.0; @@ -145,7 +146,7 @@ void SAM::Cloud (const SolverChoice& solverChoice) { erf_dtqsatw(tabs, pres, dqsatw); erf_dtqsati(tabs, pres, dqsati); - if(moisture_type == 1) { + if (SAM_moisture_type == 1) { // Cloud ice not permitted (condensation & fusion) if(tabs >= tbgmax) { omn = 1.0; @@ -160,7 +161,7 @@ void SAM::Cloud (const SolverChoice& solverChoice) { omn = an*tabs-bn; domn = an; } - } else if(moisture_type == 2) { + } else if (SAM_moisture_type == 2) { omn = 1.0; domn = 0.0; } @@ -183,16 +184,18 @@ void SAM::Cloud (const SolverChoice& solverChoice) { dtabs = -fff/dfff; tabs = tabs+dtabs; - // Update the pressure - pres = rho_array(i,j,k) * R_d * tabs - * (1.0 + R_v/R_d * qsat) * 0.01; + // For now at least we perform this iteration at constant pressure + // For the moist bubble case, the results are indistinguisable + // between running with this used vs commented out + // pres = rho_array(i,j,k) * R_d * tabs + // * (1.0 + R_v/R_d * qsat) * 0.01; // Update iteration niter = niter+1; } while(std::abs(dtabs) > tol && niter < 20); // Update qsat from last iteration (dq = dq/dt * dt) - qsat = qsat + dqsat*dtabs; + qsat += dqsat*dtabs; // Changes in each component delta_qv = qv_array(i,j,k) - qsat; @@ -209,19 +212,30 @@ void SAM::Cloud (const SolverChoice& solverChoice) { // Update temperature tabs_array(i,j,k) = tabs; - // Update pressure - pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - * (1.0 + R_v/R_d * qv_array(i,j,k)); + // Update theta from temperature (it is essential to do this BEFORE the pressure is updated) + // This would be inconsistent with updating the pressure as part of the iteration above. + // Empirically based on the moist bubble rise case, getting the correct theta here + // depends on using the old (unchanged by the phase changes) pressure. - // Update theta from temperature - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); + + // Update pressure to be consistent with updated theta_array + pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k)); + + // This was used in the earlier implmentation when we updated theta using this new pressure + // pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) + // * (1.0 + R_v/R_d * qv_array(i,j,k)); // Pressure unit conversion pres_array(i,j,k) *= 0.01; - } + // // We cannot blindly relax to qsat, but we can convert qc/qi -> qv - else { + // The concept here is that if we assume the temperature change due to this conversion + // doesn't change qsat, we can safely put all the moisture into qv without reaching qv > qsat + // because we are in the "else" part of the test on (qt_array(i,j,k) > qsat) + // + } else { // Changes in each component delta_qv = qcl_array(i,j,k) + qci_array(i,j,k); delta_qc = qcl_array(i,j,k); From cef8354a520758525f635cc04da8c70958257f05 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 25 May 2024 20:38:59 -0700 Subject: [PATCH 17/20] =?UTF-8?q?Enable=20refined=20region=20to=20be=20set?= =?UTF-8?q?=20by=20specifing=20the=20bounding=20box=20of=20the=20=E2=80=A6?= =?UTF-8?q?=20(#1631)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Enable refined region to be set by specifing the bounding box of the refined region using the fine level index space rather than the real physical positions * remove unused inputs_file --- Docs/sphinx_doc/MeshRefinement.rst | 57 ++++++++++----- .../Couette_Poiseuille/inputs_neumann_test | 72 ------------------- Source/ERF_Tagging.cpp | 55 +++++++++++++- 3 files changed, 92 insertions(+), 92 deletions(-) delete mode 100644 Exec/RegTests/Couette_Poiseuille/inputs_neumann_test diff --git a/Docs/sphinx_doc/MeshRefinement.rst b/Docs/sphinx_doc/MeshRefinement.rst index b6da7e410..efe1da1e5 100644 --- a/Docs/sphinx_doc/MeshRefinement.rst +++ b/Docs/sphinx_doc/MeshRefinement.rst @@ -23,14 +23,13 @@ See the `Gridding`_ section of the AMReX documentation for details of how indivi Static Mesh Refinement ---------------------- -For static refinement, we control the placement of grids by specifying -the low and high extents (in physical space) of each box in the lateral -directions. ERF enforces that all refinement spans the entire vertical direction. +For static refinement, we can control the placement of grids by specifying +the low and high extents (in physical space or index space) of each box. The following example demonstrates how to tag regions for static refinement. -In this first example, all cells in the region ((.15,.25,prob_lo_z)(.35,.45,prob_hi_z)) -and in the region ((.65,.75,prob_lo_z)(.85,.95,prob_hi_z)) are tagged for -one level of refinement, where prob_lo_z and prob_hi_z are the vertical extents of the domain: +In this first example, all cells in the region ((.15,.25,0.)(.35,.45,1.)) +and in the region ((.65,.75,0.0)(.85,.95,1.0)) are tagged for +one level of refinement. :: @@ -39,13 +38,13 @@ one level of refinement, where prob_lo_z and prob_hi_z are the vertical extents erf.refinement_indicators = box1 box2 - erf.box1.in_box_lo = .15 .25 - erf.box1.in_box_hi = .35 .45 + erf.box1.in_box_lo = .15 .25 0.0 + erf.box1.in_box_hi = .35 .45 1.0 - erf.box2.in_box_lo = .65 .75 - erf.box2.in_box_hi = .85 .95 + erf.box2.in_box_lo = .65 .75 0.0 + erf.box2.in_box_hi = .85 .95 1.0 -In the example below, we refine the region ((.15,.25,prob_lo_z)(.35,.45,prob_hi_z)) +In the example below, we refine the region ((.15,.25,0.)(.35,.45,.5)) by two levels of factor 3 refinement. In this case, the refined region at level 1 will be sufficient to enclose the refined region at level 2. @@ -56,11 +55,11 @@ be sufficient to enclose the refined region at level 2. erf.refinement_indicators = box1 - erf.box1.in_box_lo = .15 .25 - erf.box1.in_box_hi = .35 .45 + erf.box1.in_box_lo = .15 .25 0.0 + erf.box1.in_box_hi = .35 .45 1.0 -And in this final example, the region ((.15,.25,prob_lo_z)(.35,.45,prob_hi_z)) -will be refined by two levels of factor 3, but the larger region, ((.05,.05,prob_lo_z)(.75,.75,prob_hi_z)) +And in this final example, the region ((.15,.25,0.)(.35,.45,1.)) +will be refined by two levels of factor 3, but the larger region, ((.05,.05,0.)(.75,.75,1.)) will be refined by a single factor 3 refinement. :: @@ -70,12 +69,30 @@ will be refined by a single factor 3 refinement. erf.refinement_indicators = box1 box2 - erf.box1.in_box_lo = .15 .25 - erf.box1.in_box_hi = .35 .45 + erf.box1.in_box_lo = .15 .25 0.0 + erf.box1.in_box_hi = .35 .45 1.0 - erf.box2.in_box_lo = .05 .05 - erf.box2.in_box_hi = .75 .75 erf.box2.max_level = 1 + erf.box2.in_box_lo = .05 .05 0.0 + erf.box2.in_box_hi = .75 .75 1.0 + + +We note that instead of specifying the physical extent enclosed, we can instead specify the indices of +the bounding box of the refined region in the index space of that fine level. +To do this we use +``in_box_lo_indices`` and ``in_box_hi_indices`` instead of ``in_box_lo`` and ``in_box_hi``. +If we want to refine the inner region (spanning half the width in each direction) by one level of +factor 2 refinement, and the domain has 32x64x8 cells at level 0 covering the domain, then we would set + +:: + + amr.max_level = 1 + amr.ref_ratio = 2 + + erf.refinement_indicators = box1 + + erf.box1.in_box_lo_indices = 16 32 4 + erf.box1.in_box_hi_indices = 47 95 11 Dynamic Mesh Refinement @@ -84,6 +101,8 @@ Dynamic Mesh Refinement Dynamically created tagging functions are based on runtime data specified in the inputs file. These dynamically generated functions test on either state variables or derived variables defined in ERF_derive.cpp and included in the derive_lst in Setup.cpp. +(We note that static refinement can also be achieved by using the refinement criteria as specified below +but setting ``erf.regrid_int`` to a number greater than the total number of steps that will be taken.) Available tests include diff --git a/Exec/RegTests/Couette_Poiseuille/inputs_neumann_test b/Exec/RegTests/Couette_Poiseuille/inputs_neumann_test deleted file mode 100644 index 507aa12d5..000000000 --- a/Exec/RegTests/Couette_Poiseuille/inputs_neumann_test +++ /dev/null @@ -1,72 +0,0 @@ -# --------------------------------------------- -# --------------------------------------------- -# NOTE: to run this test effectively, you must -# hack Source/EOS.H so that getPgivenRTh -# returns just p0; this allows theta to change -# without changing the pressure and generating -# a velocity field -# --------------------------------------------- -# --------------------------------------------- - -max_step = 2000 - -amrex.fpe_trap_invalid = 1 - -fabarray.mfiter_tile_size = 1024 1024 1024 - -# PROBLEM SIZE & GEOMETRY -geometry.prob_lo = 0 0. -1. -geometry.prob_hi = 4. 1. 1. -amr.n_cell = 32 4 16 - -geometry.is_periodic = 1 1 0 - -zlo.type = "NoSlipWall" -zhi.type = "NoSlipWall" - -zlo.density = 2.0 -zlo.theta = 300.0 - -zhi.theta_grad = 3.0 - -# TIME STEP CONTROL -erf.no_substepping = 1 -erf.fixed_dt = 0.005 - -# DIAGNOSTICS & VERBOSITY -erf.sum_interval = 1 # timesteps between computing mass -erf.v = 1 # verbosity in ERF.cpp -amr.v = 1 # verbosity in Amr.cpp - -# REFINEMENT / REGRIDDING -amr.max_level = 0 # maximum level number allowed - -# CHECKPOINT FILES -erf.check_file = chk # root name of checkpoint file -erf.check_int = 1000 # number of timesteps between checkpoints - -# PLOTFILES -erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 1000 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhotheta x_velocity y_velocity z_velocity theta - -# SOLVER CHOICE -erf.use_gravity = false - -erf.alpha_T = 1.0 -erf.alpha_C = 0.0 - -erf.les_type = "None" -erf.rho0_trans = 2.0 -erf.molec_diff_type = "Constant" -erf.dynamicViscosity = 0.1 - -erf.use_coriolis = false - -erf.init_type = "uniform" - -# PROBLEM PARAMETERS -prob.prob_type = 10 - -prob.rho_0 = 2.0 -prob.T_0 = 300.0 diff --git a/Source/ERF_Tagging.cpp b/Source/ERF_Tagging.cpp index fd12b90fe..e8af57f23 100644 --- a/Source/ERF_Tagging.cpp +++ b/Source/ERF_Tagging.cpp @@ -119,11 +119,26 @@ ERF::refinement_criteria_setup () ParmParse ppr(ref_prefix); RealBox realbox; int lev_for_box; - if (ppr.countval("in_box_lo")) { + + int num_real_lo = ppr.countval("in_box_lo"); + int num_indx_lo = ppr.countval("in_box_lo_indices"); + + if ( !((num_real_lo == AMREX_SPACEDIM && num_indx_lo == 0) || + (num_indx_lo == AMREX_SPACEDIM && num_real_lo == 0) || + (num_indx_lo == 0 && num_real_lo == 0)) ) + { + amrex::Abort("Must only specify box for refinement using real OR index space"); + } + + if (num_real_lo > 0) { std::vector box_lo(3), box_hi(3); ppr.get("max_level",lev_for_box); if (lev_for_box <= max_level) { + if (n_error_buf[0] != IntVect::TheZeroVector()) { + amrex::Abort("Don't use n_error_buf > 0 when setting the box explicitly"); + } + ppr.getarr("in_box_lo",box_lo,0,AMREX_SPACEDIM); ppr.getarr("in_box_hi",box_hi,0,AMREX_SPACEDIM); realbox = RealBox(&(box_lo[0]),&(box_hi[0])); @@ -152,6 +167,44 @@ ERF::refinement_criteria_setup () } } + + } else if (num_indx_lo > 0) { + + std::vector box_lo(3), box_hi(3); + ppr.get("max_level",lev_for_box); + if (lev_for_box <= max_level) + { + if (n_error_buf[0] != IntVect::TheZeroVector()) { + amrex::Abort("Don't use n_error_buf > 0 when setting the box explicitly"); + } + + ppr.getarr("in_box_lo_indices",box_lo,0,AMREX_SPACEDIM); + ppr.getarr("in_box_hi_indices",box_hi,0,AMREX_SPACEDIM); + + Box bx(IntVect(box_lo[0],box_lo[1],box_lo[2]),IntVect(box_hi[0],box_hi[1],box_hi[2])); + amrex::Print() << "BOX " << bx << std::endl; + + const auto* dx = geom[lev_for_box].CellSize(); + const Real* plo = geom[lev_for_box].ProbLo(); + realbox = RealBox(plo[0]+ box_lo[0] *dx[0],plo[1] +box_lo[1] *dx[1],plo[2] +box_lo[2] *dx[2], + plo[0]+(box_hi[0]+1)*dx[0],plo[1]+(box_hi[1]+1)*dx[1],plo[2]+(box_hi[2]+1)*dx[2]); + + Print() << "Reading " << bx << " at level " << lev_for_box << std::endl; + num_boxes_at_level[lev_for_box] += 1; + + if ( (box_lo[0]%ref_ratio[lev_for_box-1][0] != 0) || ((box_hi[0]+1)%ref_ratio[lev_for_box-1][0] != 0) || + (box_lo[1]%ref_ratio[lev_for_box-1][1] != 0) || ((box_hi[1]+1)%ref_ratio[lev_for_box-1][1] != 0) || + (box_lo[2]%ref_ratio[lev_for_box-1][2] != 0) || ((box_hi[2]+1)%ref_ratio[lev_for_box-1][2] != 0) ) + amrex::Error("Fine box is not legit with this ref_ratio"); + boxes_at_level[lev_for_box].push_back(bx); + Print() << "Saving in 'boxes at level' as " << bx << std::endl; + } // lev + if (init_type == "real" || init_type == "metgrid") { + if (num_boxes_at_level[lev_for_box] != num_files_at_level[lev_for_box]) { + amrex::Error("Number of boxes doesnt match number of input files"); + + } + } } AMRErrorTagInfo info; From d80cad7478d483778b9b17d86f5d43b5a647bab8 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Tue, 28 May 2024 14:56:40 -0700 Subject: [PATCH 18/20] Cloud SAM (#1632) * Extend changes in SAM cloud. * Make inline function static to avoid implicit capture issue. --- Source/BoundaryConditions/MOSTStress.H | 9 +- Source/Microphysics/SAM/Cloud_SAM.cpp | 162 +++++++------------- Source/Microphysics/SAM/SAM.H | 116 ++++++++++++++ Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 6 - 4 files changed, 171 insertions(+), 122 deletions(-) diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index 8cc6f27cc..2e4a45649 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -1174,7 +1174,7 @@ struct custom_flux const amrex::Array4& /*t_surf_arr*/, const amrex::Array4& dest_arr) const { - amrex::Real rho, eta, qv; + amrex::Real rho, eta; int ic, jc; ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i; @@ -1183,7 +1183,6 @@ struct custom_flux jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc; rho = cons_arr(ic,jc,zlo,Rho_comp); - qv = cons_arr(ic,jc,zlo,RhoQ1_comp) / rho; amrex::Real qstar = q_star_arr(ic,jc,zlo); amrex::Real moflux = (std::abs(qstar) > eps) ? -rho * qstar : 0.0; @@ -1220,12 +1219,12 @@ struct custom_flux const amrex::Array4& /*vely_arr*/, const amrex::Array4& /*umm_arr*/, const amrex::Array4& /*tm_arr*/, - const amrex::Array4& u_star_arr, + const amrex::Array4& /*u_star_arr*/, const amrex::Array4& t_star_arr, const amrex::Array4& /*t_surf_arr*/, const amrex::Array4& dest_arr) const { - amrex::Real rho, eta, theta; + amrex::Real rho, eta; int ic, jc; ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i; @@ -1234,9 +1233,7 @@ struct custom_flux jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc; rho = cons_arr(ic,jc,zlo,Rho_comp); - theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho; - amrex::Real ustar = u_star_arr(ic,jc,zlo); amrex::Real tstar = t_star_arr(ic,jc,zlo); amrex::Real moflux = (std::abs(tstar) > eps) ? -rho * tstar : 0.0; amrex::Real deltaz = dz * (zlo - k); diff --git a/Source/Microphysics/SAM/Cloud_SAM.cpp b/Source/Microphysics/SAM/Cloud_SAM.cpp index 7783cb5bf..84d51549f 100644 --- a/Source/Microphysics/SAM/Cloud_SAM.cpp +++ b/Source/Microphysics/SAM/Cloud_SAM.cpp @@ -18,8 +18,6 @@ void SAM::Cloud (const SolverChoice& solverChoice) { Real fac_fus = m_fac_fus; Real rdOcp = m_rdOcp; - Real tol = 1.0e-4; - SolverChoice sc = solverChoice; int SAM_moisture_type = 1; @@ -47,16 +45,12 @@ void SAM::Cloud (const SolverChoice& solverChoice) { ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // Saturation moisture fractions - Real omn, domn; - Real qsat, dqsat; - Real qsatw, dqsatw; - Real qsati, dqsati; + Real omn; + Real qsat; + Real qsatw; + Real qsati; // Newton iteration vars - int niter; - Real fff, dfff, dtabs; - Real lstar, dlstar; - Real lstarw, lstari; Real delta_qv, delta_qc, delta_qi; // NOTE: Conversion before iterations is necessary to @@ -118,122 +112,40 @@ void SAM::Cloud (const SolverChoice& solverChoice) { pres_array(i,j,k) *= 0.01; } - // Initial guess for temperature & pressure - Real tabs = tabs_array(i,j,k); - Real pres = pres_array(i,j,k); - // Saturation moisture fractions - erf_qsatw(tabs, pres, qsatw); - erf_qsati(tabs, pres, qsati); + erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw); + erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati); qsat = omn * qsatw + (1.0-omn) * qsati; // We have enough total moisture to relax to equilibrium if (qt_array(i,j,k) > qsat) { - niter = 0; - dtabs = 1; - //================================================== - // Newton iteration to qv=qsat (cloud phase only) - //================================================== - do { - // Latent heats and their derivatives wrt to T - lstarw = fac_cond; - lstari = fac_fus; - domn = 0.0; - - // Saturation moisture fractions - erf_qsatw(tabs, pres, qsatw); - erf_qsati(tabs, pres, qsati); - erf_dtqsatw(tabs, pres, dqsatw); - erf_dtqsati(tabs, pres, dqsati); - - if (SAM_moisture_type == 1) { - // Cloud ice not permitted (condensation & fusion) - if(tabs >= tbgmax) { - omn = 1.0; - } - // Cloud water not permitted (sublimation & fusion) - else if(tabs <= tbgmin) { - omn = 0.0; - lstarw = fac_sub; - } - // Mixed cloud phase (condensation & fusion) - else { - omn = an*tabs-bn; - domn = an; - } - } else if (SAM_moisture_type == 2) { - omn = 1.0; - domn = 0.0; - } - - // Linear combination of each component - qsat = omn * qsatw + (1.0-omn) * qsati; - dqsat = omn * dqsatw + (1.0-omn) * dqsati - + domn * qsatw - domn * qsati; - lstar = omn * lstarw + (1.0-omn) * lstari; - dlstar = domn * lstarw - domn * lstari; - - // Function for root finding: - // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat) - fff = -tabs + tabs_array(i,j,k) + lstar*(qv_array(i,j,k) - qsat); - - // Derivative of function (T_new iterated on) - dfff = -1.0 + dlstar*(qv_array(i,j,k) - qsat) - lstar*dqsat; - - // Update the temperature - dtabs = -fff/dfff; - tabs = tabs+dtabs; - - // For now at least we perform this iteration at constant pressure - // For the moist bubble case, the results are indistinguisable - // between running with this used vs commented out - // pres = rho_array(i,j,k) * R_d * tabs - // * (1.0 + R_v/R_d * qsat) * 0.01; - - // Update iteration - niter = niter+1; - } while(std::abs(dtabs) > tol && niter < 20); - - // Update qsat from last iteration (dq = dq/dt * dt) - qsat += dqsat*dtabs; - - // Changes in each component - delta_qv = qv_array(i,j,k) - qsat; - delta_qc = std::max(-qcl_array(i,j,k), delta_qv * omn); - delta_qi = std::max(-qci_array(i,j,k), delta_qv * (1.0-omn)); - - // Partition the change in non-precipitating q - qv_array(i,j,k) = qsat; - qcl_array(i,j,k) += delta_qc; - qci_array(i,j,k) += delta_qi; - qn_array(i,j,k) = qcl_array(i,j,k) + qci_array(i,j,k); - qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); // Update temperature - tabs_array(i,j,k) = tabs; + tabs_array(i,j,k) = NewtonIterSat(i, j, k , SAM_moisture_type , + fac_cond , fac_fus , fac_sub , + an , bn , + tabs_array, pres_array, + qv_array , qcl_array , qci_array, + qn_array , qt_array); // Update theta from temperature (it is essential to do this BEFORE the pressure is updated) // This would be inconsistent with updating the pressure as part of the iteration above. // Empirically based on the moist bubble rise case, getting the correct theta here // depends on using the old (unchanged by the phase changes) pressure. - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); // Update pressure to be consistent with updated theta_array pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k)); - // This was used in the earlier implmentation when we updated theta using this new pressure - // pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - // * (1.0 + R_v/R_d * qv_array(i,j,k)); - // Pressure unit conversion pres_array(i,j,k) *= 0.01; // - // We cannot blindly relax to qsat, but we can convert qc/qi -> qv - // The concept here is that if we assume the temperature change due to this conversion - // doesn't change qsat, we can safely put all the moisture into qv without reaching qv > qsat - // because we are in the "else" part of the test on (qt_array(i,j,k) > qsat) + // We cannot blindly relax to qsat, but we can convert qc/qi -> qv. + // The concept here is that if we put all the moisture into qv and modify + // the temperature, we can then check if qv > qsat occurs (for final T/P/qv). + // If the reduction in T/qsat and increase in qv does trigger the + // aforementioned condition, we can do Newton iteration to drive qv = qsat. // } else { // Changes in each component @@ -251,15 +163,45 @@ void SAM::Cloud (const SolverChoice& solverChoice) { // Update temperature (endothermic since we evap/sublime) tabs_array(i,j,k) -= fac_cond * delta_qc + fac_sub * delta_qi; - // Update pressure - pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - * (1.0 + R_v/R_d * qv_array(i,j,k)); + // Update theta from temperature (it is essential to do this BEFORE the pressure is updated) + // This would be inconsistent with updating the pressure as part of the iteration above. + // Empirically based on the moist bubble rise case, getting the correct theta here + // depends on using the old (unchanged by the phase changes) pressure. + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); - // Update theta from temperature - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); + // Update pressure to be consistent with updated theta_array + pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k)); // Pressure unit conversion pres_array(i,j,k) *= 0.01; + + // Verify assumption that qv > qsat does not occur + erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw); + erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati); + qsat = omn * qsatw + (1.0-omn) * qsati; + if (qt_array(i,j,k) > qsat) { + + // Update temperature + tabs_array(i,j,k) = NewtonIterSat(i, j, k , SAM_moisture_type , + fac_cond , fac_fus , fac_sub , + an , bn , + tabs_array, pres_array, + qv_array , qcl_array , qci_array, + qn_array , qt_array); + + // Update theta from temperature (it is essential to do this BEFORE the pressure is updated) + // This would be inconsistent with updating the pressure as part of the iteration above. + // Empirically based on the moist bubble rise case, getting the correct theta here + // depends on using the old (unchanged by the phase changes) pressure. + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); + + // Update pressure to be consistent with updated theta_array + pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k)); + + // Pressure unit conversion + pres_array(i,j,k) *= 0.01; + + } } }); } // mfi diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index 3b4d9f280..8a88866a6 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -149,6 +149,122 @@ public: int Qstate_Size () override { return SAM::m_qstate_size; } + AMREX_GPU_HOST_DEVICE + AMREX_FORCE_INLINE + static amrex::Real + NewtonIterSat (int& i, int& j, int& k, + const int& SAM_moisture_type, + const amrex::Real& fac_cond, + const amrex::Real& fac_fus, + const amrex::Real& fac_sub, + const amrex::Real& an, + const amrex::Real& bn, + const amrex::Array4& tabs_array, + const amrex::Array4& pres_array, + const amrex::Array4& qv_array, + const amrex::Array4& qc_array, + const amrex::Array4& qi_array, + const amrex::Array4& qn_array, + const amrex::Array4& qt_array) + { + // Solution tolerance + amrex::Real tol = 1.0e-4; + + // Saturation moisture fractions + amrex::Real omn, domn; + amrex::Real qsat, dqsat; + amrex::Real qsatw, dqsatw; + amrex::Real qsati, dqsati; + + // Newton iteration vars + int niter; + amrex::Real fff, dfff, dtabs; + amrex::Real lstar, dlstar; + amrex::Real lstarw, lstari; + amrex::Real delta_qv, delta_qc, delta_qi; + + // Initial guess for temperature & pressure + amrex::Real tabs = tabs_array(i,j,k); + amrex::Real pres = pres_array(i,j,k); + + niter = 0; + dtabs = 1; + //================================================== + // Newton iteration to qv=qsat (cloud phase only) + //================================================== + do { + // Latent heats and their derivatives wrt to T + lstarw = fac_cond; + lstari = fac_fus; + domn = 0.0; + + // Saturation moisture fractions + erf_qsatw(tabs, pres, qsatw); + erf_qsati(tabs, pres, qsati); + erf_dtqsatw(tabs, pres, dqsatw); + erf_dtqsati(tabs, pres, dqsati); + + if (SAM_moisture_type == 1) { + // Cloud ice not permitted (condensation & fusion) + if(tabs >= tbgmax) { + omn = 1.0; + } + // Cloud water not permitted (sublimation & fusion) + else if(tabs <= tbgmin) { + omn = 0.0; + lstarw = fac_sub; + } + // Mixed cloud phase (condensation & fusion) + else { + omn = an*tabs-bn; + domn = an; + } + } else if (SAM_moisture_type == 2) { + omn = 1.0; + domn = 0.0; + } + + // Linear combination of each component + qsat = omn * qsatw + (1.0-omn) * qsati; + dqsat = omn * dqsatw + (1.0-omn) * dqsati + + domn * qsatw - domn * qsati; + lstar = omn * lstarw + (1.0-omn) * lstari; + dlstar = domn * lstarw - domn * lstari; + + // Function for root finding: + // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat) + fff = -tabs + tabs_array(i,j,k) + lstar*(qv_array(i,j,k) - qsat); + + // Derivative of function (T_new iterated on) + dfff = -1.0 + dlstar*(qv_array(i,j,k) - qsat) - lstar*dqsat; + + // Update the temperature + dtabs = -fff/dfff; + tabs += dtabs; + + // Update iteration + niter = niter+1; + } while(std::abs(dtabs) > tol && niter < 20); + + // Update qsat from last iteration (dq = dq/dt * dt) + qsat += dqsat*dtabs; + + // Changes in each component + delta_qv = qv_array(i,j,k) - qsat; + delta_qc = std::max(-qc_array(i,j,k), delta_qv * omn); + delta_qi = std::max(-qi_array(i,j,k), delta_qv * (1.0-omn)); + + // Partition the change in non-precipitating q + qv_array(i,j,k) = qsat; + qc_array(i,j,k) += delta_qc; + qi_array(i,j,k) += delta_qi; + qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); + qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); + + // Return to temperature + return tabs; + } + private: // Number of qmoist variables (qt, qv, qcl, qci, qp, qpr, qps, qpg) int m_qmoist_size = 11; diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 6b33f0e9f..059923b56 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -202,12 +202,6 @@ void erf_slow_rhs_pre (int level, int finest_level, // Define updates and fluxes in the current RK stage // ***************************************************************************** - // Open bc will be imposed upon all vars (we only access cons here for simplicity) - const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open); - const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open); - const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open); - const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open); - #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif From 3d68389156cdb3b9d758c751f53a5fdbd56b3635 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Wed, 29 May 2024 12:10:33 -0700 Subject: [PATCH 19/20] This compared well with Kessler. (#1634) Co-authored-by: Aaron Lattanzi --- .../Bubble/inputs_BF02_moist_bubble_Kessler | 81 +++++++++++++++++++ .../Bubble/inputs_BF02_moist_bubble_SAM_NoIce | 81 +++++++++++++++++++ Exec/RegTests/Bubble/prob.cpp | 3 +- Source/DataStructs/DataStruct.H | 5 +- Source/ERF.cpp | 2 +- Source/IO/Plotfile.cpp | 12 +-- Source/Microphysics/EulerianMicrophysics.H | 5 +- Source/Microphysics/Kessler/Kessler.cpp | 4 +- Source/Microphysics/Microphysics.H | 1 + Source/Microphysics/SAM/Cloud_SAM.cpp | 45 ++--------- Source/Microphysics/SAM/IceFall.cpp | 6 +- Source/Microphysics/SAM/Init_SAM.cpp | 6 +- Source/Microphysics/SAM/Precip.cpp | 41 ++++++---- Source/Microphysics/SAM/PrecipFall.cpp | 40 +++++++-- Source/Microphysics/SAM/SAM.H | 14 ++-- Source/Microphysics/SAM/Update_SAM.cpp | 3 +- 16 files changed, 267 insertions(+), 82 deletions(-) create mode 100644 Exec/RegTests/Bubble/inputs_BF02_moist_bubble_Kessler create mode 100644 Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoIce diff --git a/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_Kessler b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_Kessler new file mode 100644 index 000000000..56e38b8fc --- /dev/null +++ b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_Kessler @@ -0,0 +1,81 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 2000 +stop_time = 3600.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 20000.0 400.0 10000.0 +amr.n_cell = 200 4 100 +geometry.is_periodic = 0 1 0 +xlo.type = "SlipWall" +xhi.type = "SlipWall" +zlo.type = "SlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 0.5 +erf.fixed_mri_dt_ratio = 4 +#erf.no_substepping = 1 +#erf.fixed_dt = 0.1 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 100 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoadv_0 x_velocity y_velocity z_velocity pressure theta scalar temp pres_hse dens_hse pert_pres pert_dens eq_pot_temp qt qv qc qrain + +# SOLVER CHOICES +erf.use_gravity = true +erf.use_coriolis = false + +erf.dycore_horiz_adv_type = "Upwind_3rd" +erf.dycore_vert_adv_type = "Upwind_3rd" +erf.dryscal_horiz_adv_type = "Upwind_3rd" +erf.dryscal_vert_adv_type = "Upwind_3rd" +erf.moistscal_horiz_adv_type = "Upwind_3rd" +erf.moistscal_vert_adv_type = "Upwind_3rd" + +# PHYSICS OPTIONS +erf.les_type = "None" +erf.pbl_type = "None" +erf.moisture_model = "Kessler" +erf.buoyancy_type = 1 +erf.use_moist_background = true + +erf.molec_diff_type = "ConstantAlpha" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 0.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 0.0 # [m^2/s] +erf.alpha_C = 0.0 + +# INITIAL CONDITIONS +#erf.init_type = "input_sounding" +#erf.input_sounding_file = "BF02_moist_sounding" +#erf.init_sounding_ideal = true + +# PROBLEM PARAMETERS (optional) +# warm bubble input +prob.x_c = 10000.0 +prob.z_c = 2000.0 +prob.x_r = 2000.0 +prob.z_r = 2000.0 +prob.T_0 = 300.0 + +prob.do_moist_bubble = true +prob.theta_pert = 2.0 +prob.qt_init = 0.02 +prob.eq_pot_temp = 320.0 diff --git a/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoIce b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoIce new file mode 100644 index 000000000..7faf745ed --- /dev/null +++ b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoIce @@ -0,0 +1,81 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 2000 +stop_time = 3600.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 20000.0 400.0 10000.0 +amr.n_cell = 200 4 100 +geometry.is_periodic = 0 1 0 +xlo.type = "SlipWall" +xhi.type = "SlipWall" +zlo.type = "SlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 0.5 +erf.fixed_mri_dt_ratio = 4 +#erf.no_substepping = 1 +#erf.fixed_dt = 0.1 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 100 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoadv_0 x_velocity y_velocity z_velocity pressure theta scalar temp pres_hse dens_hse pert_pres pert_dens eq_pot_temp qt qv qc qi qp qrain qsnow qgraup + +# SOLVER CHOICES +erf.use_gravity = true +erf.use_coriolis = false + +erf.dycore_horiz_adv_type = "Upwind_3rd" +erf.dycore_vert_adv_type = "Upwind_3rd" +erf.dryscal_horiz_adv_type = "Upwind_3rd" +erf.dryscal_vert_adv_type = "Upwind_3rd" +erf.moistscal_horiz_adv_type = "Upwind_3rd" +erf.moistscal_vert_adv_type = "Upwind_3rd" + +# PHYSICS OPTIONS +erf.les_type = "None" +erf.pbl_type = "None" +erf.moisture_model = "SAM_NoIce" +erf.buoyancy_type = 1 +erf.use_moist_background = true + +erf.molec_diff_type = "ConstantAlpha" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 0.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 0.0 # [m^2/s] +erf.alpha_C = 0.0 + +# INITIAL CONDITIONS +#erf.init_type = "input_sounding" +#erf.input_sounding_file = "BF02_moist_sounding" +#erf.init_sounding_ideal = true + +# PROBLEM PARAMETERS (optional) +# warm bubble input +prob.x_c = 10000.0 +prob.z_c = 2000.0 +prob.x_r = 2000.0 +prob.z_r = 2000.0 +prob.T_0 = 300.0 + +prob.do_moist_bubble = true +prob.theta_pert = 2.0 +prob.qt_init = 0.02 +prob.eq_pot_temp = 320.0 diff --git a/Exec/RegTests/Bubble/prob.cpp b/Exec/RegTests/Bubble/prob.cpp index b6382081b..911f18aa8 100644 --- a/Exec/RegTests/Bubble/prob.cpp +++ b/Exec/RegTests/Bubble/prob.cpp @@ -364,7 +364,8 @@ Problem::init_custom_pert( if(sc.moisture_type == MoistureType::SAM) { moisture_type = 1; - } else if(sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) { + } else if(sc.moisture_type == MoistureType::SAM_NoIce || + sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) { moisture_type = 2; } diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 69a5c1b02..a3351c756 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -33,7 +33,7 @@ enum struct MoistureModelType{ }; enum struct MoistureType { - Kessler, SAM, SAM_NoPrecip_NoIce, Kessler_NoRain, None + Kessler, SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler_NoRain, None }; enum struct WindFarmType { @@ -95,6 +95,8 @@ struct SolverChoice { pp.query("moisture_model", moisture_model_string); if (moisture_model_string == "SAM") { moisture_type = MoistureType::SAM; + } else if (moisture_model_string == "SAM_NoIce") { + moisture_type = MoistureType::SAM_NoIce; } else if (moisture_model_string == "SAM_NoPrecip_NoIce") { moisture_type = MoistureType::SAM_NoPrecip_NoIce; } else if (moisture_model_string == "Kessler") { @@ -110,6 +112,7 @@ struct SolverChoice { if (moisture_type != MoistureType::None) { if (moisture_type == MoistureType::Kessler_NoRain || moisture_type == MoistureType::SAM || + moisture_type == MoistureType::SAM_NoIce || moisture_type == MoistureType::SAM_NoPrecip_NoIce) { buoyancy_type = 1; // asserted in make buoyancy } else { diff --git a/Source/ERF.cpp b/Source/ERF.cpp index c0f1effb5..68c267430 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1451,7 +1451,7 @@ ERF::ReadParameters () lsm.SetModel(); Print() << "Null land surface model!\n"; } else { - Abort("Dont know this moisture_type!") ; + Abort("Dont know this LandSurfaceType!") ; } if (verbose > 0) { diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index ada718f94..c9e0ff65e 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -77,11 +77,13 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector for (int i = 0; i < derived_names.size(); ++i) { if ( containerHasElement(plot_var_names, derived_names[i]) ) { if (solverChoice.use_terrain || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) { - if ( (solverChoice.moisture_type == MoistureType::SAM) || (derived_names[i] != "qi" && - derived_names[i] != "qsnow" && - derived_names[i] != "qgraup" && - derived_names[i] != "snow_accum" && - derived_names[i] != "graup_accum") ) + if ( (solverChoice.moisture_type == MoistureType::SAM || + solverChoice.moisture_type == MoistureType::SAM_NoIce) || + (derived_names[i] != "qi" && + derived_names[i] != "qsnow" && + derived_names[i] != "qgraup" && + derived_names[i] != "snow_accum" && + derived_names[i] != "graup_accum") ) { tmp_plot_names.push_back(derived_names[i]); } // moisture_type diff --git a/Source/Microphysics/EulerianMicrophysics.H b/Source/Microphysics/EulerianMicrophysics.H index 4e1ac1af2..ba5f72735 100644 --- a/Source/Microphysics/EulerianMicrophysics.H +++ b/Source/Microphysics/EulerianMicrophysics.H @@ -27,11 +27,12 @@ public: { AMREX_ASSERT( Microphysics::modelType(a_model_type) == MoistureModelType::Eulerian ); m_moist_model.resize(nlev); - if (a_model_type == MoistureType::SAM or + if (a_model_type == MoistureType::SAM || + a_model_type == MoistureType::SAM_NoIce || a_model_type == MoistureType::SAM_NoPrecip_NoIce) { SetModel(); amrex::Print() << "SAM moisture model!\n"; - } else if (a_model_type == MoistureType::Kessler or + } else if (a_model_type == MoistureType::Kessler || a_model_type == MoistureType::Kessler_NoRain) { SetModel(); amrex::Print() << "Kessler moisture model!\n"; diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index 7e8abbf96..f00d8759b 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -144,12 +144,12 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) autor = 0.0; if (qcc > qcw0) { - autor = 0.001; + autor = alphaelq; } accrr = 0.0; accrr = 2.2 * std::pow(qp_array(i,j,k) , 0.875); - dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - 0.001)); + dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - qcw0)); // If the amount of change is more than the amount of qc present, then dq = qc dq_clwater_to_rain = std::min(dq_clwater_to_rain, qc_array(i,j,k)); diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index ea96f0e79..dd3bc87b4 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -58,6 +58,7 @@ public: static MoistureModelType modelType (const MoistureType a_moisture_type) { if ( (a_moisture_type == MoistureType::SAM) + || (a_moisture_type == MoistureType::SAM_NoIce) || (a_moisture_type == MoistureType::SAM_NoPrecip_NoIce) || (a_moisture_type == MoistureType::Kessler) || (a_moisture_type == MoistureType::Kessler_NoRain) diff --git a/Source/Microphysics/SAM/Cloud_SAM.cpp b/Source/Microphysics/SAM/Cloud_SAM.cpp index 84d51549f..8f5ec071d 100644 --- a/Source/Microphysics/SAM/Cloud_SAM.cpp +++ b/Source/Microphysics/SAM/Cloud_SAM.cpp @@ -8,7 +8,9 @@ using namespace amrex; /** * Split cloud components according to saturation pressures; source theta from latent heat. */ -void SAM::Cloud (const SolverChoice& solverChoice) { +void +SAM::Cloud (const SolverChoice& sc) +{ constexpr Real an = 1.0/(tbgmax-tbgmin); constexpr Real bn = tbgmin*an; @@ -18,13 +20,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) { Real fac_fus = m_fac_fus; Real rdOcp = m_rdOcp; - SolverChoice sc = solverChoice; - int SAM_moisture_type = 1; - - if(solverChoice.moisture_type == MoistureType::SAM) { - SAM_moisture_type = 1; - } else if(solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) { + if (sc.moisture_type == MoistureType::SAM_NoIce || + sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) { SAM_moisture_type = 2; } @@ -128,18 +126,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) { qv_array , qcl_array , qci_array, qn_array , qt_array); - // Update theta from temperature (it is essential to do this BEFORE the pressure is updated) - // This would be inconsistent with updating the pressure as part of the iteration above. - // Empirically based on the moist bubble rise case, getting the correct theta here - // depends on using the old (unchanged by the phase changes) pressure. + // Update theta theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); - // Update pressure to be consistent with updated theta_array - pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k)); - - // Pressure unit conversion - pres_array(i,j,k) *= 0.01; - // // We cannot blindly relax to qsat, but we can convert qc/qi -> qv. // The concept here is that if we put all the moisture into qv and modify @@ -163,18 +152,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) { // Update temperature (endothermic since we evap/sublime) tabs_array(i,j,k) -= fac_cond * delta_qc + fac_sub * delta_qi; - // Update theta from temperature (it is essential to do this BEFORE the pressure is updated) - // This would be inconsistent with updating the pressure as part of the iteration above. - // Empirically based on the moist bubble rise case, getting the correct theta here - // depends on using the old (unchanged by the phase changes) pressure. + // Update theta theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); - // Update pressure to be consistent with updated theta_array - pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k)); - - // Pressure unit conversion - pres_array(i,j,k) *= 0.01; - // Verify assumption that qv > qsat does not occur erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw); erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati); @@ -189,18 +169,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) { qv_array , qcl_array , qci_array, qn_array , qt_array); - // Update theta from temperature (it is essential to do this BEFORE the pressure is updated) - // This would be inconsistent with updating the pressure as part of the iteration above. - // Empirically based on the moist bubble rise case, getting the correct theta here - // depends on using the old (unchanged by the phase changes) pressure. + // Update theta theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); - // Update pressure to be consistent with updated theta_array - pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k)); - - // Pressure unit conversion - pres_array(i,j,k) *= 0.01; - } } }); diff --git a/Source/Microphysics/SAM/IceFall.cpp b/Source/Microphysics/SAM/IceFall.cpp index 120b7c725..4951c70f9 100644 --- a/Source/Microphysics/SAM/IceFall.cpp +++ b/Source/Microphysics/SAM/IceFall.cpp @@ -7,7 +7,11 @@ using namespace amrex; /** * Sedimentation of cloud ice (A32) */ -void SAM::IceFall () { +void SAM::IceFall (const SolverChoice& sc) { + + if(sc.moisture_type == MoistureType::SAM_NoIce || + sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) + return; Real dz = m_geom.CellSize(2); Real dtn = dt; diff --git a/Source/Microphysics/SAM/Init_SAM.cpp b/Source/Microphysics/SAM/Init_SAM.cpp index 1d593044f..69abdf31f 100644 --- a/Source/Microphysics/SAM/Init_SAM.cpp +++ b/Source/Microphysics/SAM/Init_SAM.cpp @@ -18,7 +18,8 @@ using namespace amrex; * @param[in] geom Geometry associated with these MultiFabs and grids * @param[in] dt_advance Timestep for the advance */ -void SAM::Init (const MultiFab& cons_in, +void +SAM::Init (const MultiFab& cons_in, const BoxArray& grids, const Geometry& geom, const Real& dt_advance, @@ -83,7 +84,8 @@ void SAM::Init (const MultiFab& cons_in, * * @param[in] cons_in Conserved variables input */ -void SAM::Copy_State_to_Micro (const MultiFab& cons_in) +void +SAM::Copy_State_to_Micro (const MultiFab& cons_in) { // Get the temperature, density, theta, qt and qp from input for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) { diff --git a/Source/Microphysics/SAM/Precip.cpp b/Source/Microphysics/SAM/Precip.cpp index 296cc776c..6582c5290 100644 --- a/Source/Microphysics/SAM/Precip.cpp +++ b/Source/Microphysics/SAM/Precip.cpp @@ -6,7 +6,11 @@ using namespace amrex; /** * Autoconversion (A30), Accretion (A28), Evaporation (A24) */ -void SAM::Precip () { +void +SAM::Precip (const SolverChoice& sc) +{ + + if (sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) return; Real powr1 = (3.0 + b_rain) / 4.0; Real powr2 = (5.0 + b_rain) / 8.0; @@ -37,6 +41,11 @@ void SAM::Precip () { Real dtn = dt; + int SAM_moisture_type = 1; + if (sc.moisture_type == MoistureType::SAM_NoIce) { + SAM_moisture_type = 2; + } + // get the temperature, dentisy, theta, qt and qp from input for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]),TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto rho_array = mic_fab_vars[MicVar::rho]->array(mfi); @@ -77,9 +86,15 @@ void SAM::Precip () { // Work to be done for autoc/accr or evap if (qn_array(i,j,k)+qp_array(i,j,k) > 0.0) { - omn = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg)); - omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); - omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr)); + if (SAM_moisture_type == 2) { + omn = 1.0; + omp = 1.0; + omg = 0.0; + } else { + omn = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg)); + omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); + omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr)); + } qcc = qcl_array(i,j,k); qii = qci_array(i,j,k); @@ -167,12 +182,11 @@ void SAM::Precip () { qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k); - // Latent heat source for theta (endothermic fusion & exothermic melting) + // Update temperature tabs_array(i,j,k) += fac_fus * ( dqca * (1.0 - omp) - dqia * omp ); - pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - * (1.0 + R_v/R_d * qv_array(i,j,k)); - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); - pres_array(i,j,k) *= 0.01; + + // Update theta + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); } //================================================== @@ -211,12 +225,11 @@ void SAM::Precip () { qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k); - // Latent heat source for theta (endothermic) + // Update temperature tabs_array(i,j,k) -= fac_cond * dqpr + fac_sub * (dqps + dqpg); - pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - * (1.0 + R_v/R_d * qv_array(i,j,k)); - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); - pres_array(i,j,k) /= 100.0; + + // Update theta + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); } } }); diff --git a/Source/Microphysics/SAM/PrecipFall.cpp b/Source/Microphysics/SAM/PrecipFall.cpp index facf45966..d72277cc7 100644 --- a/Source/Microphysics/SAM/PrecipFall.cpp +++ b/Source/Microphysics/SAM/PrecipFall.cpp @@ -11,8 +11,11 @@ using namespace amrex; * * @param[in] hydro_type Type selection for the precipitation advection hydrodynamics scheme (0-3) */ -void SAM::PrecipFall () +void +SAM::PrecipFall (const SolverChoice& sc) { + if(sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) return; + Real rho_0 = 1.29; Real gamr3 = erf_gammafff(4.0+b_rain); @@ -49,6 +52,11 @@ void SAM::PrecipFall () MultiFab fz; fz.define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow); + int SAM_moisture_type = 1; + if (sc.moisture_type == MoistureType::SAM_NoIce) { + SAM_moisture_type = 2; + } + // Add sedimentation of precipitation field to the vert. vel. for (MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto qp_array = qp->array(mfi); @@ -80,8 +88,14 @@ void SAM::PrecipFall () Real Pprecip = 0.0; if(qp_avg > qp_threshold) { - Real omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr)); - Real omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr)); + Real omp, omg; + if (SAM_moisture_type == 2) { + omp = 1.0; + omg = 0.0; + } else { + omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr)); + omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr)); + } Real qrr = omp*qp_avg; Real qss = (1.0-omp)*(1.0-omg)*qp_avg; Real qgg = (1.0-omp)*(omg)*qp_avg; @@ -99,8 +113,14 @@ void SAM::PrecipFall () fz_array(i,j,k) = Pprecip * std::sqrt(rho_0/rho_avg); if(k==k_lo){ - Real omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr)); - Real omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr)); + Real omp, omg; + if (SAM_moisture_type == 2) { + omp = 1.0; + omg = 0.0; + } else { + omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr)); + omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr)); + } rain_accum_array(i,j,k) = rain_accum_array(i,j,k) + rho_avg*(omp*qp_avg)*vrain*dtn/rhor*1000.0; // Divide by rho_water and convert to mm snow_accum_array(i,j,k) = snow_accum_array(i,j,k) + rho_avg*(1.0-omp)*(1.0-omg)*qp_avg*vrain*dtn/rhos*1000.0; // Divide by rho_snow and convert to mm graup_accum_array(i,j,k) = graup_accum_array(i,j,k) + rho_avg*(1.0-omp)*(omg)*qp_avg*vrain*dtn/rhog*1000.0; // Divide by rho_graupel and convert to mm @@ -133,8 +153,14 @@ void SAM::PrecipFall () // Precipitating sedimentation (A19) //================================================== Real dqp = dJinv * (1.0/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef; - Real omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); - Real omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr)); + Real omp, omg; + if (SAM_moisture_type == 2) { + omp = 1.0; + omg = 0.0; + } else { + omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); + omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr)); + } qpr_array(i,j,k) = std::max(0.0, qpr_array(i,j,k) + dqp*omp); qps_array(i,j,k) = std::max(0.0, qps_array(i,j,k) + dqp*(1.0-omp)*(1.0-omg)); diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index 8a88866a6..0ba1a59c5 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -66,13 +66,13 @@ public: void Cloud (const SolverChoice& sc); // ice physics - void IceFall (); + void IceFall (const SolverChoice& sc); // precip - void Precip (); + void Precip (const SolverChoice& sc); // precip fall - void PrecipFall (); + void PrecipFall (const SolverChoice& sc); // Set up for first time void @@ -126,11 +126,9 @@ public: dt = dt_advance; this->Cloud(sc); - if(sc.moisture_type == MoistureType::SAM) { - this->IceFall(); - this->Precip(); - this->PrecipFall(); - } + this->IceFall(sc); + this->Precip(sc); + this->PrecipFall(sc); } amrex::MultiFab* diff --git a/Source/Microphysics/SAM/Update_SAM.cpp b/Source/Microphysics/SAM/Update_SAM.cpp index 3ddf97a29..ce4bb29ba 100644 --- a/Source/Microphysics/SAM/Update_SAM.cpp +++ b/Source/Microphysics/SAM/Update_SAM.cpp @@ -11,7 +11,8 @@ using namespace amrex; * @param[out] cons Conserved variables * @param[out] qmoist: qv, qc, qi, qr, qs, qg */ -void SAM::Copy_Micro_to_State (MultiFab& cons) +void +SAM::Copy_Micro_to_State (MultiFab& cons) { // Get the temperature, density, theta, qt and qp from input for ( MFIter mfi(cons,TilingIfNotGPU()); mfi.isValid(); ++mfi) { From 51a8eb0e97b8ebeb1a64b132ced9a83864bd49d2 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Fri, 31 May 2024 13:13:55 -0700 Subject: [PATCH 20/20] only call with terrain. (#1636) --- Source/BoundaryConditions/BoundaryConditions_zvel.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp index 3bd9e9930..ef04f7ee5 100644 --- a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp @@ -126,7 +126,9 @@ void ERFPhysBCFunct_w::impose_lateral_zvel_bcs (const Array4& dest_a int jflip = 2*dom_hi.y + 1 - j; if (bc_ptr_w[n].hi(1) == ERFBCType::ext_dir) { dest_arr(i,j,k) = l_bc_extdir_vals_d[n][4]; - dest_arr(i,j,k) = WFromOmega(i,j,k,dest_arr(i,j,k),xvel_arr,yvel_arr,z_phys_nd,dxInv); + if (l_use_terrain) { + dest_arr(i,j,k) = WFromOmega(i,j,k,dest_arr(i,j,k),xvel_arr,yvel_arr,z_phys_nd,dxInv); + } } else if (bc_ptr_w[n].hi(1) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,dom_hi.y,k); } else if (bc_ptr_w[n].hi(1) == ERFBCType::open) {