From 7ee67e76aeeda2bf16f5b2a24d24be5d0af8420b Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Thu, 12 Dec 2024 08:24:43 -0700 Subject: [PATCH 1/9] Update compressible dt calc to enable implicit substepping with nx=ny=1 (SCM) (#2019) --- Source/TimeIntegration/ERF_ComputeTimestep.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/Source/TimeIntegration/ERF_ComputeTimestep.cpp b/Source/TimeIntegration/ERF_ComputeTimestep.cpp index 96f37c1e4..6e91a1661 100644 --- a/Source/TimeIntegration/ERF_ComputeTimestep.cpp +++ b/Source/TimeIntegration/ERF_ComputeTimestep.cpp @@ -115,15 +115,16 @@ ERF::estTimeStep (int level, long& dt_fast_ratio) const // If we are doing implicit acoustic substepping, then the z-direction does not contribute // to the computation of the time step if (l_implicit_substepping) { - if (nxc > 1 && nyc > 1) { - new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), - ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt); - } else if (nxc > 1) { + if ((nxc > 1) && (nyc==1)) { + // 2-D in x-z new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), new_comp_dt); - } else if (nyc > 1) { + } else if ((nyc > 1) && (nxc==1)) { + // 2-D in y-z new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt); } else { - amrex::Abort("Not sure how to compute dt for this case"); + // 3-D or SCM + new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), + ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt); } // If we are not doing implicit acoustic substepping, then the z-direction contributes From 99a66902f02f674bdf7ecbcf1361292d71e17d38 Mon Sep 17 00:00:00 2001 From: Soonpil Kang <109235650+skang67@users.noreply.github.com> Date: Thu, 12 Dec 2024 07:28:56 -0800 Subject: [PATCH 2/9] Setup EB Scalar DevTest. (#2011) * Eliminated ERF_USE_EB preprocessor in AdvectionSrcForMom, because it directed to use_terrain part. * Added a scalar bubble to EB_Test. * Added ERF_redistribute to ERF class. * Fixed small things in EB_Test. * Remove whitespaces. --------- Co-authored-by: Soonpil Kang --- Exec/DevTests/EB_Test/ERF_Prob.H | 3 + Exec/DevTests/EB_Test/ERF_Prob.cpp | 66 +++++++++++++++++-- Exec/DevTests/EB_Test/inputs | 19 ++++-- ..._redistribute.cpp => ERF_Redistribute.cpp} | 36 +++++----- Source/EB/Make.package | 1 - Source/ERF.H | 12 ++++ .../TimeIntegration/ERF_TI_no_substep_fun.H | 8 ++- 7 files changed, 115 insertions(+), 30 deletions(-) rename Source/EB/{ERF_redistribute.cpp => ERF_Redistribute.cpp} (77%) diff --git a/Exec/DevTests/EB_Test/ERF_Prob.H b/Exec/DevTests/EB_Test/ERF_Prob.H index 8d670810a..9dc736590 100644 --- a/Exec/DevTests/EB_Test/ERF_Prob.H +++ b/Exec/DevTests/EB_Test/ERF_Prob.H @@ -23,6 +23,9 @@ struct ProbParm : ProbParmDefaults { amrex::Real yc_frac = 0.5; // Location of "center" of scalar (multiplies domain length) amrex::Real zc_frac = 0.5; // Location of "center" of scalar (multiplies domain length) + amrex::Real xradius = 10.0; // x-radius of scalar bubble + amrex::Real zradius = 10.0; // z-radius of scalar bubble + int prob_type = -1; }; // namespace ProbParm diff --git a/Exec/DevTests/EB_Test/ERF_Prob.cpp b/Exec/DevTests/EB_Test/ERF_Prob.cpp index 1c4ec28cc..e15bd4222 100644 --- a/Exec/DevTests/EB_Test/ERF_Prob.cpp +++ b/Exec/DevTests/EB_Test/ERF_Prob.cpp @@ -32,6 +32,9 @@ Problem::Problem() pp.query("prob_type", parms.prob_type); + pp.query("xradius", parms.xradius); + pp.query("zradius", parms.zradius); + init_base_parms(parms.rho_0, parms.T_0); } @@ -81,16 +84,71 @@ Problem::init_custom_pert( AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + // ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + // { + // // Set scalar = 0 everywhere + // state_pert(i, j, k, RhoScalar_comp) = 0.0; + + // if (use_moisture) { + // state_pert(i, j, k, RhoQ1_comp) = 0.0; + // state_pert(i, j, k, RhoQ2_comp) = 0.0; + // } + // }); + + // Set the state_pert + ParallelFor(bx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - // Set scalar = 0 everywhere - state_pert(i, j, k, RhoScalar_comp) = 0.0; + // Geometry + const Real* prob_lo = geomdata.ProbLo(); + const Real* prob_hi = geomdata.ProbHi(); + const Real* dx = geomdata.CellSize(); + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Define a point (xc,yc,zc) at the center of the domain + const Real xc = parms_d.xc_frac * (prob_lo[0] + prob_hi[0]); + const Real yc = parms_d.yc_frac * (prob_lo[1] + prob_hi[1]); + const Real zc = parms_d.zc_frac * (prob_lo[2] + prob_hi[2]); + + // Define ellipse parameters + const Real r0 = parms_d.rad_0 * (prob_hi[0] - prob_lo[0]); + const Real r3d = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); + const Real r2d_xy = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc)); + const Real r2d_xz = std::sqrt((x-xc)*(x-xc) + (z-zc)*(z-zc)); + const Real r2d_xz_nd = std::sqrt((x-xc)*(x-xc)/parms_d.xradius/parms_d.xradius + + (z-zc)*(z-zc)/parms_d.zradius/parms_d.zradius); + + if (parms_d.prob_type == 10) + { + // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain, + // + B_0*sin(x) + // state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0 * exp(-10.*r3d*r3d) + parms_d.B_0*sin(x); + state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0 * exp(-0.1*r2d_xz*r2d_xz) + parms_d.B_0*sin(x); + + } else if (parms_d.prob_type == 11) { + if (r2d_xz_nd < 1.0) + { + state_pert(i, j, k, RhoScalar_comp) = 0.5 * parms_d.A_0 * (1.0 + std::cos(PI*r2d_xz_nd)); + } else { + state_pert(i, j, k, RhoScalar_comp) = 0.0; + } + } else { + // Set scalar = A_0 in a ball of radius r0 and 0 elsewhere + if (r3d < r0) { + state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0; + } else { + state_pert(i, j, k, RhoScalar_comp) = 0.0; + } + } + + state_pert(i, j, k, RhoScalar_comp) *= parms_d.rho_0; if (use_moisture) { state_pert(i, j, k, RhoQ1_comp) = 0.0; state_pert(i, j, k, RhoQ2_comp) = 0.0; } - }); + }); // Set the x-velocity ParallelFor(xbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept diff --git a/Exec/DevTests/EB_Test/inputs b/Exec/DevTests/EB_Test/inputs index 6cf5daf48..75e3e23f7 100644 --- a/Exec/DevTests/EB_Test/inputs +++ b/Exec/DevTests/EB_Test/inputs @@ -1,6 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 100 -max_step = 10 +max_step = 800 amr.max_grid_size = 256 256 256 @@ -17,16 +16,16 @@ fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY geometry.prob_lo = 0. 0. 0. -geometry.prob_hi = 16. 1. 16. +geometry.prob_hi = 16. 4. 8. -amr.n_cell = 64 4 64 +amr.n_cell = 64 4 32 geometry.is_periodic = 0 1 0 xlo.type = Inflow xhi.type = Outflow -xlo.velocity = 1. 0. 0. +xlo.velocity = 10. 0. 0. xlo.density = 1.0 xlo.theta = 1.0 xlo.scalar = 0. @@ -36,7 +35,7 @@ zhi.type = SlipWall # TIME STEP CONTROL erf.substepping_type = None -erf.fixed_dt = 1.e-5 +erf.fixed_dt = 1.e-3 # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass @@ -72,4 +71,10 @@ erf.init_type = "uniform" # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.T_0 = 1.0 -prob.u_0 = 1.0 +prob.u_0 = 10.0 +prob.A_0 = 1.0 +prob.prob_type = 11 +prob.xradius = 3.0 +prob.zradius = 1.5 +prob.xc_frac = 0.25 +prob.zc_frac = 0.15 diff --git a/Source/EB/ERF_redistribute.cpp b/Source/EB/ERF_Redistribute.cpp similarity index 77% rename from Source/EB/ERF_redistribute.cpp rename to Source/EB/ERF_Redistribute.cpp index 184d73ad1..aa98941c1 100644 --- a/Source/EB/ERF_redistribute.cpp +++ b/Source/EB/ERF_Redistribute.cpp @@ -9,10 +9,11 @@ using namespace amrex; void ERF::redistribute_term ( int lev, - MultiFab& result, - MultiFab& result_tmp, // Saves doing a MF::copy. does this matter??? - MultiFab const& state, - BCRec const* bc) // this is bc for the state (needed for SRD slopes) + MultiFab& result, + MultiFab& result_tmp, // Saves doing a MF::copy. does this matter??? + MultiFab const& state, + BCRec const* bc, // this is bc for the state (needed for SRD slopes) + Real const dt) { // ************************************************************************ // Redistribute result_tmp and pass out result @@ -26,16 +27,17 @@ ERF::redistribute_term ( int lev, #endif for (MFIter mfi(state,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - redistribute_term(mfi, result, result_tmp, state, bc, lev); + redistribute_term(mfi, lev, result, result_tmp, state, bc, dt); } } void ERF::redistribute_term ( MFIter const& mfi, int lev, - MultiFab& result, - MultiFab& result_tmp, - MultiFab const& state, - BCRec const* bc) // this is bc for the state (needed for SRD slopes) + MultiFab& result, + MultiFab& result_tmp, + MultiFab const& state, + BCRec const* bc, // this is bc for the state (needed for SRD slopes) + Real const dt) { AMREX_ASSERT(result.nComp() == state.nComp()); @@ -57,13 +59,13 @@ ERF::redistribute_term ( MFIter const& mfi, int lev, auto const& vfrac = ebfact.getVolFrac().const_array(mfi); auto const& ccc = ebfact.getCentroid().const_array(mfi); - auto const& apx = ebfact.getAreaFrac()[0]->const_array(mfi);, - auto const& apy = ebfact.getAreaFrac()[1]->const_array(mfi);, - auto const& apz = ebfact.getAreaFrac()[2]->const_array(mfi);); + auto const& apx = ebfact.getAreaFrac()[0]->const_array(mfi); + auto const& apy = ebfact.getAreaFrac()[1]->const_array(mfi); + auto const& apz = ebfact.getAreaFrac()[2]->const_array(mfi); - auto const& fcx = ebfact.getFaceCent()[0]->const_array(mfi);, - auto const& fcy = ebfact.getFaceCent()[1]->const_array(mfi);, - auto const& fcz = ebfact.getFaceCent()[2]->const_array(mfi);); + auto const& fcx = ebfact.getFaceCent()[0]->const_array(mfi); + auto const& fcy = ebfact.getFaceCent()[1]->const_array(mfi); + auto const& fcz = ebfact.getFaceCent()[2]->const_array(mfi); Box gbx = bx; gbx.grow(3); @@ -77,7 +79,7 @@ ERF::redistribute_term ( MFIter const& mfi, int lev, // scratch(i,j,k) = 1.; // }); - std::string redistribution_type = "StateRedistribution"; + std::string redistribution_type = "StateRedist"; // State redist acts on the state. Array4 state_arr = state.const_array(mfi); @@ -85,7 +87,7 @@ ERF::redistribute_term ( MFIter const& mfi, int lev, scratch, flag, apx, apy, apz, vfrac, fcx, fcy, fcz, ccc, - bc, geom[lev], m_dt, edistribution_type); + bc, geom[lev], dt, redistribution_type); } else { diff --git a/Source/EB/Make.package b/Source/EB/Make.package index a9e8662af..8bae36093 100644 --- a/Source/EB/Make.package +++ b/Source/EB/Make.package @@ -6,6 +6,5 @@ CEXE_sources += ERF_EBCylinder.cpp CEXE_sources += ERF_EBRegular.cpp CEXE_sources += ERF_Redistribute.cpp CEXE_sources += ERF_WriteEBSurface.cpp - CEXE_headers += ERF_EBIF.H CEXE_headers += ERF_TerrainIF.H diff --git a/Source/ERF.H b/Source/ERF.H index d49fffe32..89163e184 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -475,6 +475,18 @@ public: void make_eb_box (); void make_eb_cylinder (); void make_eb_regular (); + void redistribute_term ( int lev, + amrex::MultiFab& result, + amrex::MultiFab& result_tmp, + amrex::MultiFab const& state, + amrex::BCRec const* bc, + amrex::Real const dt); + void redistribute_term ( amrex::MFIter const& mfi, int lev, + amrex::MultiFab& result, + amrex::MultiFab& result_tmp, + amrex::MultiFab const& state, + amrex::BCRec const* bc, + amrex::Real const dt); #endif // more flexible version of AverageDown() that lets you average down across multiple levels diff --git a/Source/TimeIntegration/ERF_TI_no_substep_fun.H b/Source/TimeIntegration/ERF_TI_no_substep_fun.H index 13f6e2e3f..074c0c2c7 100644 --- a/Source/TimeIntegration/ERF_TI_no_substep_fun.H +++ b/Source/TimeIntegration/ERF_TI_no_substep_fun.H @@ -151,7 +151,13 @@ } } // mfi #ifdef ERF_USE_EB - // call to redistribute_term -- pass in ssum[IntVars::cons] which is MF + // // call to redistribute_term -- pass in ssum[IntVars::cons] which is MF + // MultiFab dUdt_tmp (ba, dm, n_data, 3, MFInfo(), Factory(level)); + // dUdt_tmp.setVal(0.); + // dUdt_tmp.FillBoundary(fine_geom.periodicity()); + // const BCRec* bc_ptr_h = domain_bcs_type.data(); // Should be host or device or both? + // redistribute_term ( level, F_slow[IntVars::cons], dUdt_tmp, + // S_sum[IntVars::cons], bc_ptr_h, slow_dt); #endif } // omp From 6997ffc5c62cf83b107d3c6cb4beb7f99b8ebf19 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Thu, 12 Dec 2024 13:06:33 -0800 Subject: [PATCH 3/9] start to replace terrain_type by mesh_type (#2018) * start to replace terrain_type by mesh_type * if we read in terrain_type = Static, set mesh_type to VariableDz * update for FFT --- Source/DataStructs/ERF_DataStruct.H | 35 ++++++++++++++++- Source/ERF.H | 5 ++- Source/ERF.cpp | 39 +++++++++---------- Source/ERF_MakeNewArrays.cpp | 13 ++++--- Source/ERF_MakeNewLevel.cpp | 4 +- Source/ERF_Tagging.cpp | 2 +- Source/IO/ERF_Checkpoint.cpp | 4 +- Source/IO/ERF_Plotfile.cpp | 19 ++++----- Source/IO/ERF_WriteScalarProfiles.cpp | 2 +- Source/Initialization/ERF_Init1D.cpp | 4 +- .../LinearSolvers/ERF_ComputeDivergence.cpp | 17 ++++---- Source/LinearSolvers/ERF_PoissonSolve.cpp | 13 +++---- Source/LinearSolvers/ERF_PoissonSolve_tb.cpp | 2 +- Source/LinearSolvers/ERF_SolveWithFFT.cpp | 10 ++--- Source/TimeIntegration/ERF_MakeFastCoeffs.cpp | 8 ++-- Source/TimeIntegration/ERF_SlowRhsPost.cpp | 2 +- Source/TimeIntegration/ERF_TI_fast_headers.H | 2 +- Source/TimeIntegration/ERF_TI_substep_fun.H | 8 ++-- Source/Utils/ERF_AverageDown.cpp | 4 +- 19 files changed, 112 insertions(+), 81 deletions(-) diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 3e7349b88..865db5f64 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -29,6 +29,10 @@ AMREX_ENUM(SubsteppingType, None, Explicit, Implicit ); +AMREX_ENUM(MeshType, + ConstantDz, StretchedDz, VariableDz +); + AMREX_ENUM(TerrainType, None, Static, Moving ); @@ -87,6 +91,9 @@ struct SolverChoice { "The grid stretching ratio must be greater than 1"); } if (grid_stretching_ratio >= 1) { + if (mesh_type == MeshType::ConstantDz) { + mesh_type = MeshType::StretchedDz; + } if (terrain_type != TerrainType::Static) { amrex::Print() << "Turning terrain on to enable grid stretching" << std::endl; terrain_type = TerrainType::Static; @@ -148,10 +155,20 @@ struct SolverChoice { // Is the terrain none, static or moving? pp.query_enum_case_insensitive("terrain_type",terrain_type); + + if (terrain_type != TerrainType::None) { + mesh_type = MeshType::VariableDz; + } + int n_zlevels = pp.countval("terrain_z_levels"); - if (n_zlevels > 0 and terrain_type == TerrainType::None) + if (n_zlevels > 0) { - terrain_type = TerrainType::Static; + if (terrain_type == TerrainType::None) { + terrain_type = TerrainType::Static; + } + if (mesh_type == MeshType::ConstantDz) { + mesh_type = MeshType::StretchedDz; + } } // Use lagged_delta_rt in the fast integrator? @@ -434,6 +451,17 @@ struct SolverChoice { amrex::Print() << " None" << std::endl; } + amrex::Print() << " Mesh Type: " << std::endl; + if (mesh_type == MeshType::ConstantDz) { + amrex::Print() << " ConstantDz" << std::endl; + } else if (mesh_type == MeshType::StretchedDz) { + amrex::Print() << " StretchedDz" << std::endl; + } else if (mesh_type == MeshType::VariableDz) { + amrex::Print() << " VariableDz" << std::endl; + } else { + amrex::Abort("No mesh_type set!"); + } + amrex::Print() << "ABL Driver Type: " << std::endl; if (abl_driver_type == ABLDriverType::None) { amrex::Print() << " None" << std::endl; @@ -533,6 +561,9 @@ struct SolverChoice { inline static TerrainType terrain_type = TerrainType::None; + inline static + MeshType mesh_type = MeshType::ConstantDz; + static void set_flat_terrain_flag () { diff --git a/Source/ERF.H b/Source/ERF.H index 89163e184..2523005a5 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -826,8 +826,10 @@ private: amrex::Vector> SFS_q1fx1_lev, SFS_q1fx2_lev, SFS_q1fx3_lev; amrex::Vector> SFS_q2fx3_lev; - // Terrain / grid stretching + // Grid stretching amrex::Vector> zlevels_stag; // nominal height levels + + // Terrain data structures amrex::Vector> z_phys_nd; amrex::Vector> z_phys_cc; @@ -850,6 +852,7 @@ private: amrex::Vector> z_t_rk; + // Map scale factors amrex::Vector> mapfac_m; amrex::Vector> mapfac_u; amrex::Vector> mapfac_v; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 3c759adcf..38c87fea4 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -143,7 +143,7 @@ ERF::ERF_shared () const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1); const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2); - // This is only used when we have flat terrain and stretched z levels. + // This is only used when we have mesh_type == MeshType::StretchedDz stretched_dz_h.resize(nlevs_max); stretched_dz_d.resize(nlevs_max); @@ -159,7 +159,8 @@ ERF::ERF_shared () solverChoice.zsurf, solverChoice.dz0); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type == MeshType::StretchedDz || + SolverChoice::mesh_type == MeshType::VariableDz) { int nz = geom[0].Domain().length(2) + 1; // staggered if (std::fabs(zlevels_stag[0][nz-1]-geom[0].ProbHi(2)) > 1.0e-4) { Print() << "Note: prob_hi[2]=" << geom[0].ProbHi(2) @@ -448,7 +449,6 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) if (solverChoice.coupling_type == CouplingType::TwoWay) { - bool use_terrain = (SolverChoice::terrain_type != TerrainType::None); int ncomp = vars_new[0][Vars::cons].nComp(); for (int lev = finest_level-1; lev >= 0; lev--) { @@ -459,16 +459,16 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) const Box& bx = mfi.tilebox(); const Array4< Real> cons_arr = vars_new[lev][Vars::cons].array(mfi); const Array4 mapfac_arr = mapfac_m[lev]->const_array(mfi); - if (use_terrain) { - const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); + if (solverChoice.mesh_type == MeshType::ConstantDz) { ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - cons_arr(i,j,k,n) *= detJ_arr(i,j,k) / (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); + cons_arr(i,j,k,n) /= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); }); } else { + const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - cons_arr(i,j,k,n) /= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); + cons_arr(i,j,k,n) *= detJ_arr(i,j,k) / (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); }); } } // mfi @@ -482,16 +482,16 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) const Box& bx = mfi.tilebox(); const Array4< Real> cons_arr = vars_new[lev][Vars::cons].array(mfi); const Array4 mapfac_arr = mapfac_m[lev]->const_array(mfi); - if (use_terrain) { - const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); + if (solverChoice.mesh_type == MeshType::ConstantDz) { ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)) / detJ_arr(i,j,k); + cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); }); } else { + const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); + cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)) / detJ_arr(i,j,k); }); } } // mfi @@ -633,16 +633,16 @@ void ERF::InitData_post () { if (restart_chkfile.empty()) { - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { if (init_type == InitType::Ideal) { - Abort("We do not currently support init_type = ideal with terrain"); + Abort("We do not currently support init_type = ideal with non-constant dz"); } } // // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one // - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { for (int crse_lev = finest_level-1; crse_lev >= 0; crse_lev--) { average_down( *detJ_cc[crse_lev+1], *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev)); average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev)); @@ -704,8 +704,8 @@ ERF::InitData_post () { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 0, "Thin immersed body with refinement not currently supported."); - if (SolverChoice::terrain_type != TerrainType::None) { - amrex::Print() << "NOTE: Thin immersed body with terrain has not been tested." << std::endl; + if (SolverChoice::mesh_type != MeshType::ConstantDz) { + amrex::Print() << "NOTE: Thin immersed body with non-constant dz has not been tested." << std::endl; } } @@ -782,8 +782,8 @@ ERF::InitData_post () h_v_geos[lev], d_v_geos[lev], geom[lev], z_phys_cc[lev]); } else { - if (SolverChoice::terrain_type != TerrainType::None && !SolverChoice::terrain_is_flat) { - amrex::Print() << "Note: 1-D geostrophic wind profile input is only defined for grid stretching, not real terrain" << std::endl; + if (SolverChoice::mesh_type == MeshType::VariableDz) { + amrex::Print() << "Note: 1-D geostrophic wind profile input is not defined for real terrain" << std::endl; } init_geo_wind_profile(solverChoice.abl_geo_wind_table, h_u_geos[lev], d_u_geos[lev], @@ -875,7 +875,6 @@ ERF::InitData_post () // if (restart_chkfile == "") { - // Note -- this projection is only defined for no terrain if (solverChoice.project_initial_velocity) { Real dummy_dt = 1.0; for (int lev = 0; lev <= finest_level; ++lev) @@ -917,7 +916,7 @@ ERF::InitData_post () for (int lev = 0; lev <= finest_level; ++lev) { dz_min[lev] = geom[lev].CellSize(2); - if ( SolverChoice::terrain_type != TerrainType::None ) { + if ( SolverChoice::mesh_type != MeshType::ConstantDz ) { dz_min[lev] *= (*detJ_cc[lev]).min(0); } } diff --git a/Source/ERF_MakeNewArrays.cpp b/Source/ERF_MakeNewArrays.cpp index c0874e153..40a67bbb5 100644 --- a/Source/ERF_MakeNewArrays.cpp +++ b/Source/ERF_MakeNewArrays.cpp @@ -42,7 +42,8 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, // ******************************************************************************************** // Allocate terrain arrays // ******************************************************************************************** - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type == MeshType::StretchedDz || + SolverChoice::mesh_type == MeshType::VariableDz) { z_phys_cc[lev] = std::make_unique(ba,dm,1,1); if (solverChoice.terrain_type == TerrainType::Moving) @@ -452,7 +453,8 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap void ERF::init_zphys (int lev, Real time) { - if (SolverChoice::terrain_type != TerrainType::None) + if (SolverChoice::mesh_type == MeshType::StretchedDz || + SolverChoice::mesh_type == MeshType::VariableDz) { if (init_type != InitType::Real && init_type != InitType::Metgrid) { @@ -490,7 +492,7 @@ ERF::init_zphys (int lev, Real time) void ERF::remake_zphys (int lev, std::unique_ptr& temp_zphys_nd) { - if (lev > 0 && SolverChoice::terrain_type != TerrainType::None) + if (lev > 0 && SolverChoice::mesh_type == MeshType::VariableDz) { // // First interpolate from coarser level @@ -517,7 +519,8 @@ ERF::remake_zphys (int lev, std::unique_ptr& temp_zphys_nd) void ERF::update_terrain_arrays (int lev) { - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type == MeshType::StretchedDz || + SolverChoice::mesh_type == MeshType::VariableDz) { make_J(geom[lev],*z_phys_nd[lev],*detJ_cc[lev]); make_areas(geom[lev],*z_phys_nd[lev],*ax[lev],*ay[lev],*az[lev]); make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]); @@ -549,7 +552,7 @@ ERF::initialize_integrator (int lev, MultiFab& cons_mf, MultiFab& vel_mf) void ERF::make_physbcs (int lev) { - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type == MeshType::VariableDz) { AMREX_ALWAYS_ASSERT(z_phys_nd[lev] != nullptr); } diff --git a/Source/ERF_MakeNewLevel.cpp b/Source/ERF_MakeNewLevel.cpp index c524b2434..fdcfd4f4b 100644 --- a/Source/ERF_MakeNewLevel.cpp +++ b/Source/ERF_MakeNewLevel.cpp @@ -210,7 +210,7 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, // // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one // - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { for (int crse_lev = lev-1; crse_lev >= 0; crse_lev--) { average_down( *detJ_cc[crse_lev+1], *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev)); average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev)); @@ -339,7 +339,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one // - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { for (int crse_lev = lev-1; crse_lev >= 0; crse_lev--) { average_down( *detJ_cc[crse_lev+1], *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev)); average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev)); diff --git a/Source/ERF_Tagging.cpp b/Source/ERF_Tagging.cpp index e6d839911..db63f88eb 100644 --- a/Source/ERF_Tagging.cpp +++ b/Source/ERF_Tagging.cpp @@ -161,7 +161,7 @@ ERF::refinement_criteria_setup () jlo = static_cast((rbox_lo[1] - plo[1])/dx[1]); ihi = static_cast((rbox_hi[0] - plo[0])/dx[0]-1); jhi = static_cast((rbox_hi[1] - plo[1])/dx[1]-1); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { // Search for k indices corresponding to nominal grid // AGL heights const Box& domain = geom[lev_for_box].Domain(); diff --git a/Source/IO/ERF_Checkpoint.cpp b/Source/IO/ERF_Checkpoint.cpp index 72c3a9da1..d2294277f 100644 --- a/Source/IO/ERF_Checkpoint.cpp +++ b/Source/IO/ERF_Checkpoint.cpp @@ -145,7 +145,7 @@ ERF::WriteCheckpointFile () const MultiFab::Copy(base,base_state[lev],0,0,ncomp_base,ng_base); VisMF::Write(base, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "BaseState")); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { // Note that we also write the ghost cells of z_phys_nd IntVect ng = z_phys_nd[lev]->nGrowVect(); MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng); @@ -443,7 +443,7 @@ ERF::ReadCheckpointFile () } base_state[lev].FillBoundary(geom[lev].periodicity()); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { // Note that we also read the ghost cells of z_phys_nd IntVect ng = z_phys_nd[lev]->nGrowVect(); MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng); diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index 4b9ab1084..6378a9bbd 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -72,7 +72,8 @@ 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::terrain_type != TerrainType::None || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) { + if ( SolverChoice::mesh_type != MeshType::ConstantDz || + (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) { if ( (solverChoice.moisture_type == MoistureType::SAM || solverChoice.moisture_type == MoistureType::SAM_NoIce) || (derived_names[i] != "qi" && @@ -225,7 +226,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p // Vector of MultiFabs for nodal data Vector mf_nd(finest_level+1); - if (SolverChoice::terrain_type != TerrainType::None) { + if ( SolverChoice::mesh_type != MeshType::ConstantDz) { for (int lev = 0; lev <= finest_level; ++lev) { BoxArray nodal_grids(grids[lev]); nodal_grids.surroundingNodes(); mf_nd[lev].define(nodal_grids, dmap[lev], 3, 0); @@ -578,7 +579,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p const Array4& derdat = mf[lev].array(mfi); const Array4 & p_arr = pres.array(mfi); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -672,7 +673,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p const Array4& derdat = mf[lev].array(mfi); const Array4 & p_arr = pres.array(mfi); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -843,7 +844,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p mf_comp += 1; } // pres_hse_y - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { if (containerHasElement(plot_var_names, "z_phys")) { MultiFab::Copy(mf[lev],*z_phys_cc[lev],0,mf_comp,1,0); @@ -1350,7 +1351,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p #endif // Fill terrain distortion MF - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { for (int lev(0); lev <= finest_level; ++lev) { MultiFab::Copy(mf_nd[lev],*z_phys_nd[lev],0,2,1,0); Real dz = Geom()[lev].CellSizeArray()[2]; @@ -1399,7 +1400,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p if (plotfile_type == PlotFileType::Amrex) { Print() << "Writing native plotfile " << plotfilename << "\n"; - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, GetVecOfConstPtrs(mf), GetVecOfConstPtrs(mf_nd), @@ -1510,7 +1511,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p } Print() << "Writing plotfile " << plotfilename << "\n"; - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, GetVecOfConstPtrs(mf2), GetVecOfConstPtrs(mf_nd), @@ -1523,7 +1524,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p } } else { - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, GetVecOfConstPtrs(mf), GetVecOfConstPtrs(mf_nd), diff --git a/Source/IO/ERF_WriteScalarProfiles.cpp b/Source/IO/ERF_WriteScalarProfiles.cpp index d839640ab..a8ef59498 100644 --- a/Source/IO/ERF_WriteScalarProfiles.cpp +++ b/Source/IO/ERF_WriteScalarProfiles.cpp @@ -406,7 +406,7 @@ ERF::volWgtSumMF (int lev, auto const& dx = geom[lev].CellSizeArray(); Real cell_vol = dx[0]*dx[1]*dx[2]; volume.setVal(cell_vol); - if (solverChoice.terrain_type != TerrainType::None) { + if (solverChoice.mesh_type != MeshType::ConstantDz) { MultiFab::Multiply(volume, *detJ_cc[lev], 0, 0, 1, 0); } sum = MultiFab::Dot(tmp, 0, volume, 0, 1, 0, local); diff --git a/Source/Initialization/ERF_Init1D.cpp b/Source/Initialization/ERF_Init1D.cpp index 6af99e5c7..c1ad1e081 100644 --- a/Source/Initialization/ERF_Init1D.cpp +++ b/Source/Initialization/ERF_Init1D.cpp @@ -95,7 +95,7 @@ ERF::initHSE (int lev) std::unique_ptr new_z_phys_cc; std::unique_ptr new_z_phys_nd; - if (solverChoice.terrain_type != TerrainType::None) { + if (solverChoice.mesh_type != MeshType::ConstantDz) { new_z_phys_cc = std::make_unique(ba_new,dm_new,1,1); new_z_phys_cc->ParallelCopy(*z_phys_cc[lev],0,0,1,1,1); @@ -151,7 +151,7 @@ ERF::erf_enforce_hse (int lev, std::unique_ptr& z_cc) { Real l_gravity = solverChoice.gravity; - bool l_use_terrain = (solverChoice.terrain_type != TerrainType::None); + bool l_use_terrain = (solverChoice.mesh_type != MeshType::ConstantDz); const auto geomdata = geom[lev].data(); const Real dz = geomdata.CellSize(2); diff --git a/Source/LinearSolvers/ERF_ComputeDivergence.cpp b/Source/LinearSolvers/ERF_ComputeDivergence.cpp index 292b0a62d..fbbeb133e 100644 --- a/Source/LinearSolvers/ERF_ComputeDivergence.cpp +++ b/Source/LinearSolvers/ERF_ComputeDivergence.cpp @@ -11,8 +11,6 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Array& rho0w_arr = rho0_u_const[2]->const_array(mfi); const Array4& rhs_arr = rhs.array(mfi); - if (SolverChoice::terrain_is_flat) { // flat terrain + if (SolverChoice::mesh_type == MeshType::StretchedDz) { Real* stretched_dz_d_ptr = stretched_dz_d[lev].data(); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real dz = stretched_dz_d_ptr[k]; @@ -41,7 +43,7 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Array& mom_mf, Mult auto const dom_hi = ubound(geom[lev].Domain()); #endif - bool l_use_terrain = SolverChoice::terrain_type != TerrainType::None; - // Make sure the solver only sees the levels over which we are solving Vector ba_tmp; ba_tmp.push_back(mom_mf[Vars::cons].boxArray()); Vector dm_tmp; dm_tmp.push_back(mom_mf[Vars::cons].DistributionMap()); @@ -44,7 +42,8 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult // **************************************************************************** // #ifndef ERF_USE_EB - if (l_use_terrain && !SolverChoice::terrain_is_flat) { + if (solverChoice.mesh_type == MeshType::VariableDz) + { for ( MFIter mfi(rhs[0],TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Array4& rho0u_arr = mom_mf[IntVars::xmom].const_array(mfi); @@ -133,7 +132,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult // **************************************************************************** // No terrain or grid stretching // **************************************************************************** - if (!l_use_terrain) { + if (solverChoice.mesh_type == MeshType::ConstantDz) { #ifdef ERF_USE_FFT if (use_fft) { if (boxes_make_rectangle) { @@ -156,7 +155,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult // **************************************************************************** // Grid stretching (flat terrain) // **************************************************************************** - else if (l_use_terrain && SolverChoice::terrain_is_flat) { + else if (solverChoice.mesh_type == MeshType::StretchedDz) { #ifndef ERF_USE_FFT amrex::Abort("Rebuild with USE_FFT = TRUE so you can use the FFT solver"); #else @@ -174,7 +173,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult // **************************************************************************** // General terrain // **************************************************************************** - else if (l_use_terrain && !SolverChoice::terrain_is_flat) { + else if (solverChoice.mesh_type == MeshType::VariableDz) { #ifdef ERF_USE_FFT if (use_fft) { @@ -243,7 +242,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult // Now convert the rho0w MultiFab back to holding (rho0w) rather than Omega // **************************************************************************** // - if (l_use_terrain && !solverChoice.terrain_is_flat) + if (solverChoice.mesh_type == MeshType::VariableDz) { for (MFIter mfi(mom_mf[Vars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) { diff --git a/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp b/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp index 100fdf75b..107f6c0b7 100644 --- a/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp +++ b/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp @@ -20,7 +20,7 @@ bool ERF::projection_has_dirichlet (Array bcs) const void ERF::project_velocities_tb (int lev, Real l_dt, Vector& vmf, MultiFab& pmf) { BL_PROFILE("ERF::project_velocities_tb()"); - AMREX_ALWAYS_ASSERT(solverChoice.terrain_type == TerrainType::None); + AMREX_ALWAYS_ASSERT(solverChoice.mesh_type == MeshType::ConstantDz); // Make sure the solver only sees the levels over which we are solving Vector ba_tmp; ba_tmp.push_back(vmf[Vars::cons].boxArray()); diff --git a/Source/LinearSolvers/ERF_SolveWithFFT.cpp b/Source/LinearSolvers/ERF_SolveWithFFT.cpp index e499a7f07..22228ea40 100644 --- a/Source/LinearSolvers/ERF_SolveWithFFT.cpp +++ b/Source/LinearSolvers/ERF_SolveWithFFT.cpp @@ -12,8 +12,6 @@ void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array 0) { amrex::Print() << "Using the 3D FFT solver..." << std::endl; @@ -49,7 +46,8 @@ void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array 0) { amrex::Print() << "Using the hybrid FFT solver..." << std::endl; @@ -89,7 +87,7 @@ void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array const& fz_arr = fluxes[2].array(mfi); - if (l_use_terrain && SolverChoice::terrain_is_flat) { + if (solverChoice.mesh_type == MeshType::StretchedDz) { Real* stretched_dz_d_ptr = stretched_dz_d[lev].data(); ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { diff --git a/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp b/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp index e5ef6b0bc..f724406d7 100644 --- a/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp +++ b/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp @@ -14,7 +14,7 @@ using namespace amrex; * @param[in] S_stage_prim primitive variables (i.e. conserved variables divided by density) at the last stage * @param[in] pi_stage Exner function at the last stage * @param[in] geom Container for geometric information - * @param[in] terrain_type Are we using terrain-fitted coordinates + * @param[in] mesh_type Do we have constant dz? * @param[in] gravity Magnitude of gravity * @param[in] c_p Coefficient at constant pressure * @param[in] r0 Reference (hydrostatically stratified) density @@ -30,7 +30,7 @@ void make_fast_coeffs (int /*level*/, const MultiFab& pi_stage, // Exner function evaluated at least stage const amrex::Geometry geom, bool l_use_moisture, - TerrainType terrain_type, + MeshType mesh_type, Real gravity, Real c_p, std::unique_ptr& detJ_cc, const MultiFab* r0, const MultiFab* pi0, @@ -76,7 +76,7 @@ void make_fast_coeffs (int /*level*/, const Array4 & stage_cons = S_stage_data[IntVars::cons].const_array(mfi); const Array4 & prim = S_stage_prim.const_array(mfi); - const Array4& detJ = (terrain_type != TerrainType::None) ? detJ_cc->const_array(mfi) : Array4{}; + const Array4& detJ = (mesh_type != MeshType::ConstantDz) ? detJ_cc->const_array(mfi) : Array4{}; const Array4& r0_ca = r0->const_array(mfi); const Array4& pi0_ca = pi0->const_array(mfi); @@ -107,7 +107,7 @@ void make_fast_coeffs (int /*level*/, Real halfg = std::abs(0.5 * grav_gpu[2]); //Note we don't act on the bottom or top boundaries of the domain - if (terrain_type != TerrainType::None) + if (mesh_type != MeshType::ConstantDz) { ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k) { diff --git a/Source/TimeIntegration/ERF_SlowRhsPost.cpp b/Source/TimeIntegration/ERF_SlowRhsPost.cpp index 2cb7c9f60..aaee5f3fe 100644 --- a/Source/TimeIntegration/ERF_SlowRhsPost.cpp +++ b/Source/TimeIntegration/ERF_SlowRhsPost.cpp @@ -110,7 +110,7 @@ void erf_slow_rhs_post (int level, int finest_level, const MultiFab* t_mean_mf = nullptr; if (most) t_mean_mf = most->get_mac_avg(0,2); - const bool l_use_terrain = (solverChoice.terrain_type != TerrainType::None); + const bool l_use_terrain = (solverChoice.mesh_type != MeshType::ConstantDz); const bool l_moving_terrain = (solverChoice.terrain_type == TerrainType::Moving); const bool l_reflux = (solverChoice.coupling_type != CouplingType::OneWay); if (l_moving_terrain) AMREX_ALWAYS_ASSERT(l_use_terrain); diff --git a/Source/TimeIntegration/ERF_TI_fast_headers.H b/Source/TimeIntegration/ERF_TI_fast_headers.H index e21fd3c74..a0731419e 100644 --- a/Source/TimeIntegration/ERF_TI_fast_headers.H +++ b/Source/TimeIntegration/ERF_TI_fast_headers.H @@ -112,7 +112,7 @@ void make_fast_coeffs (int level, const amrex::MultiFab& pi_stage, const amrex::Geometry geom, const bool use_moisture, - const TerrainType terrain_type, + const MeshType mesh_type, const amrex::Real gravity, const amrex::Real c_p, std::unique_ptr& detJ_cc, diff --git a/Source/TimeIntegration/ERF_TI_substep_fun.H b/Source/TimeIntegration/ERF_TI_substep_fun.H index e336631eb..0717fd577 100644 --- a/Source/TimeIntegration/ERF_TI_substep_fun.H +++ b/Source/TimeIntegration/ERF_TI_substep_fun.H @@ -93,7 +93,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, // We have to call this each step since it depends on the substep time now // Note we pass in the *old* detJ here make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, - l_use_moisture, SolverChoice::terrain_type, solverChoice.gravity, solverChoice.c_p, + l_use_moisture, SolverChoice::mesh_type, solverChoice.gravity, solverChoice.c_p, detJ_cc[level], r0, pi0, dtau, beta_s, phys_bc_type); if (fast_step == 0) { @@ -121,12 +121,12 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, mapfac_m[level], mapfac_u[level], mapfac_v[level], fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping); } - } else if (solverChoice.terrain_type == TerrainType::Static) { + } else if (solverChoice.mesh_type != MeshType::ConstantDz) { if (fast_step == 0) { // If this is the first substep we make the coefficients since they are based only on stage data make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, - l_use_moisture, SolverChoice::terrain_type, solverChoice.gravity, solverChoice.c_p, + l_use_moisture, SolverChoice::mesh_type, solverChoice.gravity, solverChoice.c_p, detJ_cc[level], r0, pi0, dtau, beta_s, phys_bc_type); // If this is the first substep we pass in S_old as the previous step's solution @@ -152,7 +152,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, // If this is the first substep we make the coefficients since they are based only on stage data make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, - l_use_moisture, SolverChoice::terrain_type, solverChoice.gravity, solverChoice.c_p, + l_use_moisture, SolverChoice::mesh_type, solverChoice.gravity, solverChoice.c_p, detJ_cc[level], r0, pi0, dtau, beta_s, phys_bc_type); // If this is the first substep we pass in S_old as the previous step's solution diff --git a/Source/Utils/ERF_AverageDown.cpp b/Source/Utils/ERF_AverageDown.cpp index c080d1e1b..1c28e439b 100644 --- a/Source/Utils/ERF_AverageDown.cpp +++ b/Source/Utils/ERF_AverageDown.cpp @@ -55,7 +55,7 @@ ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT const Box& bx = mfi.tilebox(); const Array4< Real> cons_arr = vars_new[lev][Vars::cons].array(mfi); const Array4 mapfac_arr = mapfac_m[lev]->const_array(mfi); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { @@ -119,7 +119,7 @@ ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT const Box& bx = mfi.tilebox(); const Array4< Real> cons_arr = vars_new[lev][Vars::cons].array(mfi); const Array4 mapfac_arr = mapfac_m[lev]->const_array(mfi); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { From 2f908f31d6f7b150a89afdc43dddf85d6f3a776f Mon Sep 17 00:00:00 2001 From: Debojyoti Ghosh Date: Thu, 12 Dec 2024 13:06:47 -0800 Subject: [PATCH 4/9] added input parameter for choosing the stable Redistribute() option; moved fixing random seed to the beginning of the constructor (#2021) --- Source/ERF.cpp | 10 ++++++++++ Source/Initialization/ERF_InitCustom.cpp | 9 --------- Source/Particles/ERFPC.H | 2 ++ Source/Particles/ERFPCInitializations.cpp | 9 +++++++++ 4 files changed, 21 insertions(+), 9 deletions(-) diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 38c87fea4..5c7e898fe 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -84,6 +85,15 @@ Vector BCNames = {"xlo", "ylo", "zlo", "xhi", "yhi", "zhi"}; // - initializes BCRec boundary condition object ERF::ERF () { + int fix_random_seed = 0; + ParmParse pp("erf"); pp.query("fix_random_seed", fix_random_seed); + // Note that the value of 1024UL is not significant -- the point here is just to set the + // same seed for all MPI processes for the purpose of regression testing + if (fix_random_seed) { + Print() << "Fixing the random seed" << std::endl; + InitRandom(1024UL); + } + ERF_shared(); } diff --git a/Source/Initialization/ERF_InitCustom.cpp b/Source/Initialization/ERF_InitCustom.cpp index f7b7f1eb7..4d7d4ed44 100644 --- a/Source/Initialization/ERF_InitCustom.cpp +++ b/Source/Initialization/ERF_InitCustom.cpp @@ -42,15 +42,6 @@ ERF::init_custom (int lev) yvel_pert.setVal(0.); zvel_pert.setVal(0.); - int fix_random_seed = 0; - ParmParse pp("erf"); pp.query("fix_random_seed", fix_random_seed); - // Note that the value of 1024UL is not significant -- the point here is just to set the - // same seed for all MPI processes for the purpose of regression testing - if (fix_random_seed) { - Print() << "Fixing the random seed" << std::endl; - InitRandom(1024UL); - } - #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif diff --git a/Source/Particles/ERFPC.H b/Source/Particles/ERFPC.H index 986e04ac0..dd3ef7cd9 100644 --- a/Source/Particles/ERFPC.H +++ b/Source/Particles/ERFPC.H @@ -223,6 +223,8 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdxAoS::ncomps, std::string m_initialization_type; /*!< initial particle distribution type */ int m_ppc_init; /*!< initial number of particles per cell */ + bool m_stable_redistribute; /*!< use stable redistribute for deterministic simulations */ + /*! read inputs from file */ virtual void readInputs (); diff --git a/Source/Particles/ERFPCInitializations.cpp b/Source/Particles/ERFPCInitializations.cpp index feca8c476..3d4abcc91 100644 --- a/Source/Particles/ERFPCInitializations.cpp +++ b/Source/Particles/ERFPCInitializations.cpp @@ -49,6 +49,14 @@ void ERFPC::readInputs () m_advect_w_gravity = (m_name == ERFParticleNames::hydro ? true : false); pp.query("advect_with_gravity", m_advect_w_gravity); + m_stable_redistribute = false; + pp.query("stable_redistribute", m_stable_redistribute); + + if (m_stable_redistribute) { + Print() << "Note: using stable Redistribute() for particle container " << m_name << ".\n"; + } + setStableRedistribute(m_stable_redistribute); + return; } @@ -65,6 +73,7 @@ void ERFPC::InitializeParticles (const std::unique_ptr& a_height_ptr) << m_name << " particle species.\n"; Error("See error message!"); } + return; } From 9159695011cedc4b5411877974e7404249f903ce Mon Sep 17 00:00:00 2001 From: Harish Date: Thu, 12 Dec 2024 14:07:03 -0700 Subject: [PATCH 5/9] Immersed Forcing (#2017) * First commit * Delete Exec/ABL/Backtrace.0 * Some cleaning * Bug fix * Whitespace cleaning * Some more whitespace * Making terrain list GPU aware * Plot IB mask. * merge origin. --------- Co-authored-by: AMLattanzi --- CMake/BuildERFExe.cmake | 1 + Exec/ABL/erf_terrain_def | 1024 ++++++++++++++++++ Exec/ABL/inputs_terrain | 63 ++ Exec/ABL/write_terrain.py | 15 + Source/DataStructs/ERF_DataStruct.H | 2 + Source/ERF.H | 6 +- Source/ERF.cpp | 16 +- Source/ERF_MakeNewLevel.cpp | 3 + Source/IO/ERF_Plotfile.cpp | 18 + Source/SourceTerms/ERF_MakeMomSources.cpp | 57 + Source/SourceTerms/ERF_SrcHeaders.H | 1 + Source/SourceTerms/ERF_TerrainDrag.H | 33 + Source/SourceTerms/ERF_TerrainDrag.cpp | 97 ++ Source/SourceTerms/Make.package | 2 + Source/TimeIntegration/ERF_TI_slow_rhs_fun.H | 12 +- 15 files changed, 1345 insertions(+), 5 deletions(-) create mode 100644 Exec/ABL/erf_terrain_def create mode 100644 Exec/ABL/inputs_terrain create mode 100644 Exec/ABL/write_terrain.py create mode 100644 Source/SourceTerms/ERF_TerrainDrag.H create mode 100644 Source/SourceTerms/ERF_TerrainDrag.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 0a926b1bb..47e1e8cab 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -170,6 +170,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/SourceTerms/ERF_MoistSetRhs.cpp ${SRC_DIR}/SourceTerms/ERF_NumericalDiffusion.cpp ${SRC_DIR}/SourceTerms/ERF_ForestDrag.cpp + ${SRC_DIR}/SourceTerms/ERF_TerrainDrag.cpp ${SRC_DIR}/TimeIntegration/ERF_ComputeTimestep.cpp ${SRC_DIR}/TimeIntegration/ERF_Advance.cpp ${SRC_DIR}/TimeIntegration/ERF_TimeStep.cpp diff --git a/Exec/ABL/erf_terrain_def b/Exec/ABL/erf_terrain_def new file mode 100644 index 000000000..1a649774c --- /dev/null +++ b/Exec/ABL/erf_terrain_def @@ -0,0 +1,1024 @@ +768 768 108 +784 768 108 +800 768 108 +816 768 108 +832 768 108 +848 768 108 +864 768 108 +880 768 108 +896 768 108 +912 768 108 +928 768 108 +944 768 108 +960 768 108 +976 768 108 +992 768 108 +1008 768 108 +1024 768 108 +1040 768 108 +1056 768 108 +1072 768 108 +1088 768 108 +1104 768 108 +1120 768 108 +1136 768 108 +1152 768 108 +1168 768 108 +1184 768 108 +1200 768 108 +1216 768 108 +1232 768 108 +1248 768 108 +1264 768 108 +768 784 108 +784 784 108 +800 784 108 +816 784 108 +832 784 108 +848 784 108 +864 784 108 +880 784 108 +896 784 108 +912 784 108 +928 784 108 +944 784 108 +960 784 108 +976 784 108 +992 784 108 +1008 784 108 +1024 784 108 +1040 784 108 +1056 784 108 +1072 784 108 +1088 784 108 +1104 784 108 +1120 784 108 +1136 784 108 +1152 784 108 +1168 784 108 +1184 784 108 +1200 784 108 +1216 784 108 +1232 784 108 +1248 784 108 +1264 784 108 +768 800 108 +784 800 108 +800 800 108 +816 800 108 +832 800 108 +848 800 108 +864 800 108 +880 800 108 +896 800 108 +912 800 108 +928 800 108 +944 800 108 +960 800 108 +976 800 108 +992 800 108 +1008 800 108 +1024 800 108 +1040 800 108 +1056 800 108 +1072 800 108 +1088 800 108 +1104 800 108 +1120 800 108 +1136 800 108 +1152 800 108 +1168 800 108 +1184 800 108 +1200 800 108 +1216 800 108 +1232 800 108 +1248 800 108 +1264 800 108 +768 816 108 +784 816 108 +800 816 108 +816 816 108 +832 816 108 +848 816 108 +864 816 108 +880 816 108 +896 816 108 +912 816 108 +928 816 108 +944 816 108 +960 816 108 +976 816 108 +992 816 108 +1008 816 108 +1024 816 108 +1040 816 108 +1056 816 108 +1072 816 108 +1088 816 108 +1104 816 108 +1120 816 108 +1136 816 108 +1152 816 108 +1168 816 108 +1184 816 108 +1200 816 108 +1216 816 108 +1232 816 108 +1248 816 108 +1264 816 108 +768 832 108 +784 832 108 +800 832 108 +816 832 108 +832 832 108 +848 832 108 +864 832 108 +880 832 108 +896 832 108 +912 832 108 +928 832 108 +944 832 108 +960 832 108 +976 832 108 +992 832 108 +1008 832 108 +1024 832 108 +1040 832 108 +1056 832 108 +1072 832 108 +1088 832 108 +1104 832 108 +1120 832 108 +1136 832 108 +1152 832 108 +1168 832 108 +1184 832 108 +1200 832 108 +1216 832 108 +1232 832 108 +1248 832 108 +1264 832 108 +768 848 108 +784 848 108 +800 848 108 +816 848 108 +832 848 108 +848 848 108 +864 848 108 +880 848 108 +896 848 108 +912 848 108 +928 848 108 +944 848 108 +960 848 108 +976 848 108 +992 848 108 +1008 848 108 +1024 848 108 +1040 848 108 +1056 848 108 +1072 848 108 +1088 848 108 +1104 848 108 +1120 848 108 +1136 848 108 +1152 848 108 +1168 848 108 +1184 848 108 +1200 848 108 +1216 848 108 +1232 848 108 +1248 848 108 +1264 848 108 +768 864 108 +784 864 108 +800 864 108 +816 864 108 +832 864 108 +848 864 108 +864 864 108 +880 864 108 +896 864 108 +912 864 108 +928 864 108 +944 864 108 +960 864 108 +976 864 108 +992 864 108 +1008 864 108 +1024 864 108 +1040 864 108 +1056 864 108 +1072 864 108 +1088 864 108 +1104 864 108 +1120 864 108 +1136 864 108 +1152 864 108 +1168 864 108 +1184 864 108 +1200 864 108 +1216 864 108 +1232 864 108 +1248 864 108 +1264 864 108 +768 880 108 +784 880 108 +800 880 108 +816 880 108 +832 880 108 +848 880 108 +864 880 108 +880 880 108 +896 880 108 +912 880 108 +928 880 108 +944 880 108 +960 880 108 +976 880 108 +992 880 108 +1008 880 108 +1024 880 108 +1040 880 108 +1056 880 108 +1072 880 108 +1088 880 108 +1104 880 108 +1120 880 108 +1136 880 108 +1152 880 108 +1168 880 108 +1184 880 108 +1200 880 108 +1216 880 108 +1232 880 108 +1248 880 108 +1264 880 108 +768 896 108 +784 896 108 +800 896 108 +816 896 108 +832 896 108 +848 896 108 +864 896 108 +880 896 108 +896 896 108 +912 896 108 +928 896 108 +944 896 108 +960 896 108 +976 896 108 +992 896 108 +1008 896 108 +1024 896 108 +1040 896 108 +1056 896 108 +1072 896 108 +1088 896 108 +1104 896 108 +1120 896 108 +1136 896 108 +1152 896 108 +1168 896 108 +1184 896 108 +1200 896 108 +1216 896 108 +1232 896 108 +1248 896 108 +1264 896 108 +768 912 108 +784 912 108 +800 912 108 +816 912 108 +832 912 108 +848 912 108 +864 912 108 +880 912 108 +896 912 108 +912 912 108 +928 912 108 +944 912 108 +960 912 108 +976 912 108 +992 912 108 +1008 912 108 +1024 912 108 +1040 912 108 +1056 912 108 +1072 912 108 +1088 912 108 +1104 912 108 +1120 912 108 +1136 912 108 +1152 912 108 +1168 912 108 +1184 912 108 +1200 912 108 +1216 912 108 +1232 912 108 +1248 912 108 +1264 912 108 +768 928 108 +784 928 108 +800 928 108 +816 928 108 +832 928 108 +848 928 108 +864 928 108 +880 928 108 +896 928 108 +912 928 108 +928 928 108 +944 928 108 +960 928 108 +976 928 108 +992 928 108 +1008 928 108 +1024 928 108 +1040 928 108 +1056 928 108 +1072 928 108 +1088 928 108 +1104 928 108 +1120 928 108 +1136 928 108 +1152 928 108 +1168 928 108 +1184 928 108 +1200 928 108 +1216 928 108 +1232 928 108 +1248 928 108 +1264 928 108 +768 944 108 +784 944 108 +800 944 108 +816 944 108 +832 944 108 +848 944 108 +864 944 108 +880 944 108 +896 944 108 +912 944 108 +928 944 108 +944 944 108 +960 944 108 +976 944 108 +992 944 108 +1008 944 108 +1024 944 108 +1040 944 108 +1056 944 108 +1072 944 108 +1088 944 108 +1104 944 108 +1120 944 108 +1136 944 108 +1152 944 108 +1168 944 108 +1184 944 108 +1200 944 108 +1216 944 108 +1232 944 108 +1248 944 108 +1264 944 108 +768 960 108 +784 960 108 +800 960 108 +816 960 108 +832 960 108 +848 960 108 +864 960 108 +880 960 108 +896 960 108 +912 960 108 +928 960 108 +944 960 108 +960 960 108 +976 960 108 +992 960 108 +1008 960 108 +1024 960 108 +1040 960 108 +1056 960 108 +1072 960 108 +1088 960 108 +1104 960 108 +1120 960 108 +1136 960 108 +1152 960 108 +1168 960 108 +1184 960 108 +1200 960 108 +1216 960 108 +1232 960 108 +1248 960 108 +1264 960 108 +768 976 108 +784 976 108 +800 976 108 +816 976 108 +832 976 108 +848 976 108 +864 976 108 +880 976 108 +896 976 108 +912 976 108 +928 976 108 +944 976 108 +960 976 108 +976 976 108 +992 976 108 +1008 976 108 +1024 976 108 +1040 976 108 +1056 976 108 +1072 976 108 +1088 976 108 +1104 976 108 +1120 976 108 +1136 976 108 +1152 976 108 +1168 976 108 +1184 976 108 +1200 976 108 +1216 976 108 +1232 976 108 +1248 976 108 +1264 976 108 +768 992 108 +784 992 108 +800 992 108 +816 992 108 +832 992 108 +848 992 108 +864 992 108 +880 992 108 +896 992 108 +912 992 108 +928 992 108 +944 992 108 +960 992 108 +976 992 108 +992 992 108 +1008 992 108 +1024 992 108 +1040 992 108 +1056 992 108 +1072 992 108 +1088 992 108 +1104 992 108 +1120 992 108 +1136 992 108 +1152 992 108 +1168 992 108 +1184 992 108 +1200 992 108 +1216 992 108 +1232 992 108 +1248 992 108 +1264 992 108 +768 1008 108 +784 1008 108 +800 1008 108 +816 1008 108 +832 1008 108 +848 1008 108 +864 1008 108 +880 1008 108 +896 1008 108 +912 1008 108 +928 1008 108 +944 1008 108 +960 1008 108 +976 1008 108 +992 1008 108 +1008 1008 108 +1024 1008 108 +1040 1008 108 +1056 1008 108 +1072 1008 108 +1088 1008 108 +1104 1008 108 +1120 1008 108 +1136 1008 108 +1152 1008 108 +1168 1008 108 +1184 1008 108 +1200 1008 108 +1216 1008 108 +1232 1008 108 +1248 1008 108 +1264 1008 108 +768 1024 108 +784 1024 108 +800 1024 108 +816 1024 108 +832 1024 108 +848 1024 108 +864 1024 108 +880 1024 108 +896 1024 108 +912 1024 108 +928 1024 108 +944 1024 108 +960 1024 108 +976 1024 108 +992 1024 108 +1008 1024 108 +1024 1024 108 +1040 1024 108 +1056 1024 108 +1072 1024 108 +1088 1024 108 +1104 1024 108 +1120 1024 108 +1136 1024 108 +1152 1024 108 +1168 1024 108 +1184 1024 108 +1200 1024 108 +1216 1024 108 +1232 1024 108 +1248 1024 108 +1264 1024 108 +768 1040 108 +784 1040 108 +800 1040 108 +816 1040 108 +832 1040 108 +848 1040 108 +864 1040 108 +880 1040 108 +896 1040 108 +912 1040 108 +928 1040 108 +944 1040 108 +960 1040 108 +976 1040 108 +992 1040 108 +1008 1040 108 +1024 1040 108 +1040 1040 108 +1056 1040 108 +1072 1040 108 +1088 1040 108 +1104 1040 108 +1120 1040 108 +1136 1040 108 +1152 1040 108 +1168 1040 108 +1184 1040 108 +1200 1040 108 +1216 1040 108 +1232 1040 108 +1248 1040 108 +1264 1040 108 +768 1056 108 +784 1056 108 +800 1056 108 +816 1056 108 +832 1056 108 +848 1056 108 +864 1056 108 +880 1056 108 +896 1056 108 +912 1056 108 +928 1056 108 +944 1056 108 +960 1056 108 +976 1056 108 +992 1056 108 +1008 1056 108 +1024 1056 108 +1040 1056 108 +1056 1056 108 +1072 1056 108 +1088 1056 108 +1104 1056 108 +1120 1056 108 +1136 1056 108 +1152 1056 108 +1168 1056 108 +1184 1056 108 +1200 1056 108 +1216 1056 108 +1232 1056 108 +1248 1056 108 +1264 1056 108 +768 1072 108 +784 1072 108 +800 1072 108 +816 1072 108 +832 1072 108 +848 1072 108 +864 1072 108 +880 1072 108 +896 1072 108 +912 1072 108 +928 1072 108 +944 1072 108 +960 1072 108 +976 1072 108 +992 1072 108 +1008 1072 108 +1024 1072 108 +1040 1072 108 +1056 1072 108 +1072 1072 108 +1088 1072 108 +1104 1072 108 +1120 1072 108 +1136 1072 108 +1152 1072 108 +1168 1072 108 +1184 1072 108 +1200 1072 108 +1216 1072 108 +1232 1072 108 +1248 1072 108 +1264 1072 108 +768 1088 108 +784 1088 108 +800 1088 108 +816 1088 108 +832 1088 108 +848 1088 108 +864 1088 108 +880 1088 108 +896 1088 108 +912 1088 108 +928 1088 108 +944 1088 108 +960 1088 108 +976 1088 108 +992 1088 108 +1008 1088 108 +1024 1088 108 +1040 1088 108 +1056 1088 108 +1072 1088 108 +1088 1088 108 +1104 1088 108 +1120 1088 108 +1136 1088 108 +1152 1088 108 +1168 1088 108 +1184 1088 108 +1200 1088 108 +1216 1088 108 +1232 1088 108 +1248 1088 108 +1264 1088 108 +768 1104 108 +784 1104 108 +800 1104 108 +816 1104 108 +832 1104 108 +848 1104 108 +864 1104 108 +880 1104 108 +896 1104 108 +912 1104 108 +928 1104 108 +944 1104 108 +960 1104 108 +976 1104 108 +992 1104 108 +1008 1104 108 +1024 1104 108 +1040 1104 108 +1056 1104 108 +1072 1104 108 +1088 1104 108 +1104 1104 108 +1120 1104 108 +1136 1104 108 +1152 1104 108 +1168 1104 108 +1184 1104 108 +1200 1104 108 +1216 1104 108 +1232 1104 108 +1248 1104 108 +1264 1104 108 +768 1120 108 +784 1120 108 +800 1120 108 +816 1120 108 +832 1120 108 +848 1120 108 +864 1120 108 +880 1120 108 +896 1120 108 +912 1120 108 +928 1120 108 +944 1120 108 +960 1120 108 +976 1120 108 +992 1120 108 +1008 1120 108 +1024 1120 108 +1040 1120 108 +1056 1120 108 +1072 1120 108 +1088 1120 108 +1104 1120 108 +1120 1120 108 +1136 1120 108 +1152 1120 108 +1168 1120 108 +1184 1120 108 +1200 1120 108 +1216 1120 108 +1232 1120 108 +1248 1120 108 +1264 1120 108 +768 1136 108 +784 1136 108 +800 1136 108 +816 1136 108 +832 1136 108 +848 1136 108 +864 1136 108 +880 1136 108 +896 1136 108 +912 1136 108 +928 1136 108 +944 1136 108 +960 1136 108 +976 1136 108 +992 1136 108 +1008 1136 108 +1024 1136 108 +1040 1136 108 +1056 1136 108 +1072 1136 108 +1088 1136 108 +1104 1136 108 +1120 1136 108 +1136 1136 108 +1152 1136 108 +1168 1136 108 +1184 1136 108 +1200 1136 108 +1216 1136 108 +1232 1136 108 +1248 1136 108 +1264 1136 108 +768 1152 108 +784 1152 108 +800 1152 108 +816 1152 108 +832 1152 108 +848 1152 108 +864 1152 108 +880 1152 108 +896 1152 108 +912 1152 108 +928 1152 108 +944 1152 108 +960 1152 108 +976 1152 108 +992 1152 108 +1008 1152 108 +1024 1152 108 +1040 1152 108 +1056 1152 108 +1072 1152 108 +1088 1152 108 +1104 1152 108 +1120 1152 108 +1136 1152 108 +1152 1152 108 +1168 1152 108 +1184 1152 108 +1200 1152 108 +1216 1152 108 +1232 1152 108 +1248 1152 108 +1264 1152 108 +768 1168 108 +784 1168 108 +800 1168 108 +816 1168 108 +832 1168 108 +848 1168 108 +864 1168 108 +880 1168 108 +896 1168 108 +912 1168 108 +928 1168 108 +944 1168 108 +960 1168 108 +976 1168 108 +992 1168 108 +1008 1168 108 +1024 1168 108 +1040 1168 108 +1056 1168 108 +1072 1168 108 +1088 1168 108 +1104 1168 108 +1120 1168 108 +1136 1168 108 +1152 1168 108 +1168 1168 108 +1184 1168 108 +1200 1168 108 +1216 1168 108 +1232 1168 108 +1248 1168 108 +1264 1168 108 +768 1184 108 +784 1184 108 +800 1184 108 +816 1184 108 +832 1184 108 +848 1184 108 +864 1184 108 +880 1184 108 +896 1184 108 +912 1184 108 +928 1184 108 +944 1184 108 +960 1184 108 +976 1184 108 +992 1184 108 +1008 1184 108 +1024 1184 108 +1040 1184 108 +1056 1184 108 +1072 1184 108 +1088 1184 108 +1104 1184 108 +1120 1184 108 +1136 1184 108 +1152 1184 108 +1168 1184 108 +1184 1184 108 +1200 1184 108 +1216 1184 108 +1232 1184 108 +1248 1184 108 +1264 1184 108 +768 1200 108 +784 1200 108 +800 1200 108 +816 1200 108 +832 1200 108 +848 1200 108 +864 1200 108 +880 1200 108 +896 1200 108 +912 1200 108 +928 1200 108 +944 1200 108 +960 1200 108 +976 1200 108 +992 1200 108 +1008 1200 108 +1024 1200 108 +1040 1200 108 +1056 1200 108 +1072 1200 108 +1088 1200 108 +1104 1200 108 +1120 1200 108 +1136 1200 108 +1152 1200 108 +1168 1200 108 +1184 1200 108 +1200 1200 108 +1216 1200 108 +1232 1200 108 +1248 1200 108 +1264 1200 108 +768 1216 108 +784 1216 108 +800 1216 108 +816 1216 108 +832 1216 108 +848 1216 108 +864 1216 108 +880 1216 108 +896 1216 108 +912 1216 108 +928 1216 108 +944 1216 108 +960 1216 108 +976 1216 108 +992 1216 108 +1008 1216 108 +1024 1216 108 +1040 1216 108 +1056 1216 108 +1072 1216 108 +1088 1216 108 +1104 1216 108 +1120 1216 108 +1136 1216 108 +1152 1216 108 +1168 1216 108 +1184 1216 108 +1200 1216 108 +1216 1216 108 +1232 1216 108 +1248 1216 108 +1264 1216 108 +768 1232 108 +784 1232 108 +800 1232 108 +816 1232 108 +832 1232 108 +848 1232 108 +864 1232 108 +880 1232 108 +896 1232 108 +912 1232 108 +928 1232 108 +944 1232 108 +960 1232 108 +976 1232 108 +992 1232 108 +1008 1232 108 +1024 1232 108 +1040 1232 108 +1056 1232 108 +1072 1232 108 +1088 1232 108 +1104 1232 108 +1120 1232 108 +1136 1232 108 +1152 1232 108 +1168 1232 108 +1184 1232 108 +1200 1232 108 +1216 1232 108 +1232 1232 108 +1248 1232 108 +1264 1232 108 +768 1248 108 +784 1248 108 +800 1248 108 +816 1248 108 +832 1248 108 +848 1248 108 +864 1248 108 +880 1248 108 +896 1248 108 +912 1248 108 +928 1248 108 +944 1248 108 +960 1248 108 +976 1248 108 +992 1248 108 +1008 1248 108 +1024 1248 108 +1040 1248 108 +1056 1248 108 +1072 1248 108 +1088 1248 108 +1104 1248 108 +1120 1248 108 +1136 1248 108 +1152 1248 108 +1168 1248 108 +1184 1248 108 +1200 1248 108 +1216 1248 108 +1232 1248 108 +1248 1248 108 +1264 1248 108 +768 1264 108 +784 1264 108 +800 1264 108 +816 1264 108 +832 1264 108 +848 1264 108 +864 1264 108 +880 1264 108 +896 1264 108 +912 1264 108 +928 1264 108 +944 1264 108 +960 1264 108 +976 1264 108 +992 1264 108 +1008 1264 108 +1024 1264 108 +1040 1264 108 +1056 1264 108 +1072 1264 108 +1088 1264 108 +1104 1264 108 +1120 1264 108 +1136 1264 108 +1152 1264 108 +1168 1264 108 +1184 1264 108 +1200 1264 108 +1216 1264 108 +1232 1264 108 +1248 1264 108 +1264 1264 108 diff --git a/Exec/ABL/inputs_terrain b/Exec/ABL/inputs_terrain new file mode 100644 index 000000000..e91f8c3f1 --- /dev/null +++ b/Exec/ABL/inputs_terrain @@ -0,0 +1,63 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 1 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 2048 2048 1024 +amr.n_cell = 64 64 32 + +geometry.is_periodic = 1 1 0 + +zlo.type = "SlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.cfl = 0.5 + +# 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 = 512 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 512 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta terrain_IB_mask + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "None" +erf.les_type = "Smagorinsky" +erf.Cs = 0.1 + +erf.init_type = "uniform" + +erf.terrain_file = erf_terrain_def + +# PROBLEM PARAMETERS +prob.rho_0 = 1.0 +prob.A_0 = 1.0 + +prob.U_0 = 10.0 +prob.V_0 = 0.0 +prob.W_0 = 0.0 +prob.T_0 = 300.0 + +# Higher values of perturbations lead to instability +# Instability seems to be coming from BC +prob.U_0_Pert_Mag = 0.0 +prob.V_0_Pert_Mag = 0.0 # +prob.W_0_Pert_Mag = 0.0 diff --git a/Exec/ABL/write_terrain.py b/Exec/ABL/write_terrain.py new file mode 100644 index 000000000..e85e89ff7 --- /dev/null +++ b/Exec/ABL/write_terrain.py @@ -0,0 +1,15 @@ +import numpy as np + +target=open("erf_terrain_def","w") + +xc=1024 +yc=1024 +x=np.arange(xc-256,xc+256,16) +y=np.arange(yc-256,yc+256,16) + +X,Y=np.meshgrid(x,y) + +for i in range(0,X.shape[0]): + for j in range(0,X.shape[1]): + target.write("%g %g %g\n"%(X[i,j],Y[i,j],108)) +target.close() diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 865db5f64..bae93281f 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -690,5 +690,7 @@ struct SolverChoice { // Flag for valid canopy model bool do_forest {false}; + // Immersed Forcing + bool do_terrain {false}; }; #endif diff --git a/Source/ERF.H b/Source/ERF.H index 2523005a5..0d38af78d 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -46,6 +46,7 @@ #include #include #include +#include #ifdef ERF_USE_PARTICLES #include "ERF_ParticleData.H" @@ -980,7 +981,9 @@ private: "Lturb", // moisture vars "qt", "qv", "qc", "qi", "qp", "qrain", "qsnow", "qgraup", "qsat", - "rain_accum", "snow_accum", "graup_accum" + "rain_accum", "snow_accum", "graup_accum", + // Terrain IB mask + "terrain_IB_mask" #ifdef ERF_USE_EB // EB variables ,"volfrac", @@ -1155,6 +1158,7 @@ private: std::unique_ptr m_r2d = nullptr; std::unique_ptr m_most = nullptr; amrex::Vector> m_forest; + amrex::Vector> m_terrain; // // Holds info for dynamically generated tagging criteria diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 5c7e898fe..fbf06fd48 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -141,7 +141,12 @@ ERF::ERF_shared () // NOTE: size canopy model before readparams (if file exists, we construct) m_forest.resize(nlevs_max); - for (int lev = 0; lev < max_level; ++lev) { m_forest[lev] = nullptr; } + for (int lev = 0; lev < max_level; ++lev) { m_forest[lev] = nullptr;} + + // Immersed Forcing + m_terrain.resize(nlevs_max); + for (int lev = 0; lev < max_level; ++lev) { m_terrain[lev] = nullptr;} + ReadParameters(); initializeMicrophysics(nlevs_max); @@ -1630,6 +1635,15 @@ ERF::ReadParameters () } } + //Query the terrain file name + std::string terrainfile; + solverChoice.do_terrain = pp.query("terrain_file", terrainfile); + if (solverChoice.do_terrain) { + for (int lev = 0; lev <= max_level; ++lev) { + m_terrain[lev] = std::make_unique(terrainfile); + } + } + // AmrMesh iterate on grids? bool iterate(true); pp_amr.query("iterate_grids",iterate); diff --git a/Source/ERF_MakeNewLevel.cpp b/Source/ERF_MakeNewLevel.cpp index fdcfd4f4b..1f2f5993b 100644 --- a/Source/ERF_MakeNewLevel.cpp +++ b/Source/ERF_MakeNewLevel.cpp @@ -143,6 +143,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, // ******************************************************************************************** if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_terrain) { m_terrain[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } //******************************************************************************************** // Microphysics // ******************************************************************************************* @@ -222,6 +223,7 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, // ******************************************************************************************** if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_terrain) { m_terrain[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } //******************************************************************************************** // Microphysics // ******************************************************************************************* @@ -351,6 +353,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // ******************************************************************************************** if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_terrain) { m_terrain[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } // ***************************************************************************************************** // Create the physbcs objects (after initializing the terrain but before calling FillCoarsePatch // ***************************************************************************************************** diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index 6378a9bbd..8d558b6e3 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -489,6 +489,24 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p mf_comp ++; } + if (containerHasElement(plot_var_names, "terrain_IB_mask")) + { + MultiFab* terrain_blank = m_terrain[lev]->get_terrain_blank_field(); +#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& src = terrain_blank->const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + derdat(i, j, k, mf_comp) = src(i,j,k); + }); + } + mf_comp ++; + } + #ifdef ERF_USE_WINDFARM if (containerHasElement(plot_var_names, "num_turb") and (solverChoice.windfarm_type == WindFarmType::Fitch or solverChoice.windfarm_type == WindFarmType::EWP or diff --git a/Source/SourceTerms/ERF_MakeMomSources.cpp b/Source/SourceTerms/ERF_MakeMomSources.cpp index 8133f0482..7dd8d4d3e 100644 --- a/Source/SourceTerms/ERF_MakeMomSources.cpp +++ b/Source/SourceTerms/ERF_MakeMomSources.cpp @@ -53,6 +53,7 @@ void make_mom_sources (int level, MultiFab & zmom_src, const MultiFab& base_state, MultiFab* forest_drag, + MultiFab* terrain_blank, const Geometry geom, const SolverChoice& solverChoice, std::unique_ptr& /*mapfac_m*/, @@ -86,10 +87,21 @@ void make_mom_sources (int level, // 6. nudging towards input sounding data // 7. numerical diffusion for (xmom,ymom,zmom) // 8. sponge + // 9. Forest canopy + // 10. Immersed Forcing // ***************************************************************************** const bool l_use_ndiff = solverChoice.use_NumDiff; const bool use_terrain = solverChoice.terrain_type != TerrainType::None; const bool l_do_forest = solverChoice.do_forest; + const bool l_do_terrain = solverChoice.do_terrain; + + // Check if terrain and immersed terrain clash + if(use_terrain && l_do_terrain){ + amrex::Error(" Cannot use immersed forcing with terrain"); + } + if(l_do_forest && l_do_terrain){ + amrex::Error(" Currently forest canopy cannot be used with immersed forcing"); + } // ***************************************************************************** // Data for Coriolis forcing @@ -228,6 +240,8 @@ void make_mom_sources (int level, const Array4& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) : Array4{}; + const Array4& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) : + Array4{}; const Array4& z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) : Array4{}; @@ -526,5 +540,48 @@ void make_mom_sources (int level, zmom_src_arr(i, j, k) -= f_drag * uz * windspeed; }); } + // ***************************************************************************** + // 10. Add Immersed source terms + // ***************************************************************************** + if (l_do_terrain) { + const Real drag_coefficient=10.0/dz; + const Real tiny = std::numeric_limits::epsilon(); + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const Real ux = u(i, j, k); + const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k ) + + v(i, j+1, k ) + v(i-1, j+1, k ) ); + const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k ) + + w(i, j , k+1) + w(i-1, j , k+1) ); + const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz); + const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k)); + const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0); + xmom_src_arr(i, j, k) -= t_blank * CdM * ux * windspeed; + }); + ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k ) + + u(i+1, j , k ) + u(i+1, j-1, k ) ); + const Real uy = v(i, j, k); + const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k ) + + w(i , j , k+1) + w(i , j-1, k+1) ); + const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz); + const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k)); + const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0); + ymom_src_arr(i, j, k) -= t_blank * CdM * uy * windspeed; + }); + ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k ) + + u(i , j , k-1) + u(i+1, j , k-1) ); + const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k ) + + v(i , j , k-1) + v(i , j+1, k-1) ); + const amrex::Real uz = w(i, j, k); + const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz); + const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1)); + const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0); + zmom_src_arr(i, j, k) -= t_blank * CdM * uz * windspeed; + }); + } } // mfi } diff --git a/Source/SourceTerms/ERF_SrcHeaders.H b/Source/SourceTerms/ERF_SrcHeaders.H index 5ea93913d..4992e80cd 100644 --- a/Source/SourceTerms/ERF_SrcHeaders.H +++ b/Source/SourceTerms/ERF_SrcHeaders.H @@ -59,6 +59,7 @@ void make_mom_sources (int level, int nrk, amrex::MultiFab& zmom_source, const amrex::MultiFab& base_state, amrex::MultiFab* forest_drag, + amrex::MultiFab* terrain_blank, const amrex::Geometry geom, const SolverChoice& solverChoice, std::unique_ptr& mapfac_m, diff --git a/Source/SourceTerms/ERF_TerrainDrag.H b/Source/SourceTerms/ERF_TerrainDrag.H new file mode 100644 index 000000000..289db4a55 --- /dev/null +++ b/Source/SourceTerms/ERF_TerrainDrag.H @@ -0,0 +1,33 @@ +#ifndef ERF_TERRAINDRAG_H_ +#define ERF_TERRAINDRAG_H_ + +#include +#include + +/* + Adapted from: https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020MS002141 + */ +class TerrainDrag +{ +public: + + explicit TerrainDrag (std::string terrainfile); + + ~TerrainDrag () = default; + + void + define_terrain_blank_field (const amrex::BoxArray& ba, + const amrex::DistributionMapping& dm, + amrex::Geometry& geom, + amrex::MultiFab* z_phys_nd); + + amrex::MultiFab* + get_terrain_blank_field () { return m_terrain_blank.get(); } + +private: + amrex::Vector m_x_terrain; + amrex::Vector m_y_terrain; + amrex::Vector m_height_terrain; + std::unique_ptr m_terrain_blank; +}; +#endif diff --git a/Source/SourceTerms/ERF_TerrainDrag.cpp b/Source/SourceTerms/ERF_TerrainDrag.cpp new file mode 100644 index 000000000..36bc04b07 --- /dev/null +++ b/Source/SourceTerms/ERF_TerrainDrag.cpp @@ -0,0 +1,97 @@ +#include + +using namespace amrex; + +/* + Constructor to get the terrain parameters: + TreeType xc, yc, height +*/ +TerrainDrag::TerrainDrag(std::string terrainfile) +{ + std::ifstream file(terrainfile, std::ios::in); + if (!file.good()) { + Abort("Cannot find terrain file: " + terrainfile); + } + // xc yc height + Real value1, value2, value3; + while (file >> value1 >> value2 >> value3) { + m_x_terrain.push_back(value1); + m_y_terrain.push_back(value2); + m_height_terrain.push_back(value3); + } + file.close(); +} + +void +TerrainDrag::define_terrain_blank_field( + const BoxArray& ba, + const DistributionMapping& dm, + Geometry& geom, + MultiFab* z_phys_nd) +{ + // Geometry params + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + + // Allocate the terrain blank MF + // NOTE: 1 ghost cell for averaging to faces + m_terrain_blank.reset(); + m_terrain_blank = std::make_unique(ba, dm, 1, 1); + m_terrain_blank->setVal(0.); + const auto xterrain_size = m_x_terrain.size(); + amrex::Gpu::DeviceVector d_xterrain(xterrain_size); + amrex::Gpu::DeviceVector d_yterrain(xterrain_size); + amrex::Gpu::DeviceVector d_zterrain(xterrain_size); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_x_terrain.begin(), m_x_terrain.end(), + d_xterrain.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_y_terrain.begin(), m_y_terrain.end(), + d_yterrain.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_height_terrain.begin(), m_height_terrain.end(), + d_zterrain.begin()); + const auto* xterrain_ptr = d_xterrain.data(); + const auto* yterrain_ptr = d_yterrain.data(); + const auto* zterrain_ptr = d_zterrain.data(); + // Set the terrain blank data + for (MFIter mfi(*m_terrain_blank); mfi.isValid(); ++mfi) { + Box gtbx = mfi.growntilebox(); + const Array4& levelBlank = m_terrain_blank->array(mfi); + const Array4& z_nd = + (z_phys_nd) ? z_phys_nd->const_array(mfi) : Array4{}; + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Loop over terrain points + for (unsigned ii = 0; ii < xterrain_size; ++ii) { + Real ht = zterrain_ptr[ii]; + Real xt = xterrain_ptr[ii]; + Real yt = yterrain_ptr[ii]; + // Physical positions of cell-centers + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + Real z = prob_lo[2] + (k + 0.5) * dx[2]; + if (z_nd) { + z = 0.125 * (z_nd(i, j, k) + z_nd(i + 1, j, k) + z_nd(i, j + 1, k) + + z_nd(i + 1, j + 1, k) + z_nd(i, j, k + 1) + + z_nd(i + 1, j, k + 1) + z_nd(i, j + 1, k + 1) + + z_nd(i + 1, j + 1, k + 1)); + } + z = std::max(z, 0.0); + const Real cell_radius = std::sqrt(dx[0] * dx[0] + dx[1] * dx[1]); + // Proximity to the terrain + const Real radius = + std::sqrt((x - xt) * (x - xt) + (y - yt) * (y - yt)); + const Real terrain_point = + (radius <= cell_radius && z <= ht) ? 1.0 : 0.0; + levelBlank(i, j, k) = terrain_point; + if(terrain_point==1){ + break; + } + } + }); + } // mfi + +// Fillboundary for periodic ghost cell copy +m_terrain_blank->FillBoundary(geom.periodicity()); + +} // init_terrain_blank_field diff --git a/Source/SourceTerms/Make.package b/Source/SourceTerms/Make.package index d70c5c9e0..92b757433 100644 --- a/Source/SourceTerms/Make.package +++ b/Source/SourceTerms/Make.package @@ -6,6 +6,7 @@ CEXE_sources += ERF_ApplySpongeZoneBCs.cpp CEXE_sources += ERF_ApplySpongeZoneBCs_ReadFromFile.cpp CEXE_sources += ERF_NumericalDiffusion.cpp CEXE_sources += ERF_ForestDrag.cpp +CEXE_sources += ERF_TerrainDrag.cpp ifeq ($(USE_NETCDF),TRUE) CEXE_sources += ERF_MoistSetRhs.cpp @@ -15,3 +16,4 @@ CEXE_headers += ERF_NumericalDiffusion.H CEXE_headers += ERF_SrcHeaders.H CEXE_headers += ERF_BuoyancyUtils.H CEXE_headers += ERF_ForestDrag.H +CEXE_headers += ERF_TerrainDrag.H \ No newline at end of file diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H index 314e1c766..77c8998f4 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H @@ -63,6 +63,9 @@ // Canopy data for mom sources MultiFab* forest_drag = nullptr; if (solverChoice.do_forest) { forest_drag = m_forest[level]->get_drag_field(); } + // Immersed Forcing + MultiFab* terrain_blank = nullptr; + if(solverChoice.do_terrain) { terrain_blank = m_terrain[level]->get_terrain_blank_field(); } // Moving terrain if ( solverChoice.terrain_type == TerrainType::Moving ) @@ -131,7 +134,7 @@ z_phys_nd[level], z_phys_cc[level], xvel_new, yvel_new, zvel_new, xmom_src, ymom_src, zmom_src, - base_state_new[level], forest_drag, fine_geom, solverChoice, + base_state_new[level], forest_drag, terrain_blank, fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, @@ -240,7 +243,7 @@ z_phys_nd[level], z_phys_cc[level], xvel_new, yvel_new, zvel_new, xmom_src, ymom_src, zmom_src, - base_state[level], forest_drag, fine_geom, solverChoice, + base_state[level], forest_drag, terrain_blank, fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, @@ -433,6 +436,9 @@ // Canopy data for mom sources MultiFab* forest_drag = nullptr; if (solverChoice.do_forest) { forest_drag = m_forest[level]->get_drag_field(); } + // Immersed Forcing + MultiFab* terrain_blank = nullptr; + if(solverChoice.do_terrain) { terrain_blank = m_terrain[level]->get_terrain_blank_field(); } make_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, cc_src, z_phys_cc[level], #if defined(ERF_USE_RRTMGP) @@ -449,7 +455,7 @@ z_phys_nd[level], z_phys_cc[level], xvel_new, yvel_new, zvel_new, xmom_src, ymom_src, zmom_src, - base_state[level], forest_drag, fine_geom, solverChoice, + base_state[level], forest_drag, terrain_blank, fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, From 617c1ca5a664f457a604741c733344f059f76bb6 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Thu, 12 Dec 2024 13:29:38 -0800 Subject: [PATCH 6/9] forest-->forest_drag and terrain-->terrain_drag (#2022) --- Source/DataStructs/ERF_DataStruct.H | 9 ++++---- Source/ERF.H | 4 ++-- Source/ERF.cpp | 22 ++++++++++---------- Source/ERF_MakeNewLevel.cpp | 12 +++++------ Source/IO/ERF_Plotfile.cpp | 2 +- Source/SourceTerms/ERF_MakeMomSources.cpp | 22 ++++++++++---------- Source/TimeIntegration/ERF_TI_slow_rhs_fun.H | 8 +++---- 7 files changed, 40 insertions(+), 39 deletions(-) diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index bae93281f..0be9f91e5 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -688,9 +688,10 @@ struct SolverChoice { amrex::Real windfarm_x_shift = -1.0; amrex::Real windfarm_y_shift = -1.0; - // Flag for valid canopy model - bool do_forest {false}; - // Immersed Forcing - bool do_terrain {false}; + // Use forest canopy model? + bool do_forest_drag {false}; + + // Use immersed forcing representation of terrain? + bool do_terrain_drag {false}; }; #endif diff --git a/Source/ERF.H b/Source/ERF.H index 0d38af78d..06370d69b 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -1157,8 +1157,8 @@ private: std::unique_ptr m_w2d = nullptr; std::unique_ptr m_r2d = nullptr; std::unique_ptr m_most = nullptr; - amrex::Vector> m_forest; - amrex::Vector> m_terrain; + amrex::Vector> m_forest_drag; + amrex::Vector> m_terrain_drag; // // Holds info for dynamically generated tagging criteria diff --git a/Source/ERF.cpp b/Source/ERF.cpp index fbf06fd48..1c29ccaad 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -140,12 +140,12 @@ ERF::ERF_shared () lsm_flux.resize(nlevs_max); // NOTE: size canopy model before readparams (if file exists, we construct) - m_forest.resize(nlevs_max); - for (int lev = 0; lev < max_level; ++lev) { m_forest[lev] = nullptr;} + m_forest_drag.resize(nlevs_max); + for (int lev = 0; lev < max_level; ++lev) { m_forest_drag[lev] = nullptr;} - // Immersed Forcing - m_terrain.resize(nlevs_max); - for (int lev = 0; lev < max_level; ++lev) { m_terrain[lev] = nullptr;} + // Immersed Forcing Representation of Terrain + m_terrain_drag.resize(nlevs_max); + for (int lev = 0; lev < max_level; ++lev) { m_terrain_drag[lev] = nullptr;} ReadParameters(); @@ -1628,19 +1628,19 @@ ERF::ReadParameters () // Query the canopy model file name std::string forestfile; - solverChoice.do_forest = pp.query("forest_file", forestfile); - if (solverChoice.do_forest) { + solverChoice.do_forest_drag = pp.query("forest_file", forestfile); + if (solverChoice.do_forest_drag) { for (int lev = 0; lev <= max_level; ++lev) { - m_forest[lev] = std::make_unique(forestfile); + m_forest_drag[lev] = std::make_unique(forestfile); } } //Query the terrain file name std::string terrainfile; - solverChoice.do_terrain = pp.query("terrain_file", terrainfile); - if (solverChoice.do_terrain) { + solverChoice.do_terrain_drag = pp.query("terrain_file", terrainfile); + if (solverChoice.do_terrain_drag) { for (int lev = 0; lev <= max_level; ++lev) { - m_terrain[lev] = std::make_unique(terrainfile); + m_terrain_drag[lev] = std::make_unique(terrainfile); } } diff --git a/Source/ERF_MakeNewLevel.cpp b/Source/ERF_MakeNewLevel.cpp index 1f2f5993b..51ceac884 100644 --- a/Source/ERF_MakeNewLevel.cpp +++ b/Source/ERF_MakeNewLevel.cpp @@ -141,9 +141,9 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, // ******************************************************************************************** // Build the data structures for canopy model (depends upon z_phys) // ******************************************************************************************** - if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_forest_drag) { m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } - if (solverChoice.do_terrain) { m_terrain[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_terrain_drag) { m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } //******************************************************************************************** // Microphysics // ******************************************************************************************* @@ -221,9 +221,9 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, // ******************************************************************************************** // Build the data structures for canopy model (depends upon z_phys) // ******************************************************************************************** - if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_forest_drag) { m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } - if (solverChoice.do_terrain) { m_terrain[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_terrain_drag) { m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } //******************************************************************************************** // Microphysics // ******************************************************************************************* @@ -351,9 +351,9 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // ******************************************************************************************** // Build the data structures for canopy model (depends upon z_phys) // ******************************************************************************************** - if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_forest_drag) { m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } - if (solverChoice.do_terrain) { m_terrain[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_terrain_drag) { m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } // ***************************************************************************************************** // Create the physbcs objects (after initializing the terrain but before calling FillCoarsePatch // ***************************************************************************************************** diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index 8d558b6e3..fb1d331e6 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -491,7 +491,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p if (containerHasElement(plot_var_names, "terrain_IB_mask")) { - MultiFab* terrain_blank = m_terrain[lev]->get_terrain_blank_field(); + MultiFab* terrain_blank = m_terrain_drag[lev]->get_terrain_blank_field(); #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif diff --git a/Source/SourceTerms/ERF_MakeMomSources.cpp b/Source/SourceTerms/ERF_MakeMomSources.cpp index 7dd8d4d3e..a6c05c654 100644 --- a/Source/SourceTerms/ERF_MakeMomSources.cpp +++ b/Source/SourceTerms/ERF_MakeMomSources.cpp @@ -90,16 +90,16 @@ void make_mom_sources (int level, // 9. Forest canopy // 10. Immersed Forcing // ***************************************************************************** - const bool l_use_ndiff = solverChoice.use_NumDiff; - const bool use_terrain = solverChoice.terrain_type != TerrainType::None; - const bool l_do_forest = solverChoice.do_forest; - const bool l_do_terrain = solverChoice.do_terrain; + const bool l_use_ndiff = solverChoice.use_NumDiff; + const bool l_use_zphys = (solverChoice.mesh_type != MeshType::ConstantDz); + const bool l_do_forest_drag = solverChoice.do_forest_drag; + const bool l_do_terrain_drag = solverChoice.do_terrain_drag; // Check if terrain and immersed terrain clash - if(use_terrain && l_do_terrain){ - amrex::Error(" Cannot use immersed forcing with terrain"); + if(l_use_zphys && l_do_terrain_drag){ + amrex::Error(" Cannot use immersed forcing for terrain with terrain-fitted coordinates"); } - if(l_do_forest && l_do_terrain){ + if(l_do_forest_drag && l_do_terrain_drag){ amrex::Error(" Currently forest canopy cannot be used with immersed forcing"); } @@ -243,9 +243,9 @@ void make_mom_sources (int level, const Array4& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) : Array4{}; - const Array4& z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) : + const Array4& z_nd_arr = (l_use_zphys) ? z_phys_nd->const_array(mfi) : Array4{}; - const Array4& z_cc_arr = (use_terrain) ? z_phys_cc->const_array(mfi) : + const Array4& z_cc_arr = (l_use_zphys) ? z_phys_cc->const_array(mfi) : Array4{}; // ***************************************************************************** @@ -505,7 +505,7 @@ void make_mom_sources (int level, // ***************************************************************************** // 9. Add CANOPY source terms // ***************************************************************************** - if (l_do_forest) { + if (l_do_forest_drag) { ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const Real ux = u(i, j, k); @@ -543,7 +543,7 @@ void make_mom_sources (int level, // ***************************************************************************** // 10. Add Immersed source terms // ***************************************************************************** - if (l_do_terrain) { + if (l_do_terrain_drag) { const Real drag_coefficient=10.0/dz; const Real tiny = std::numeric_limits::epsilon(); ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H index 77c8998f4..5a6b09a6f 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H @@ -62,10 +62,10 @@ // Canopy data for mom sources MultiFab* forest_drag = nullptr; - if (solverChoice.do_forest) { forest_drag = m_forest[level]->get_drag_field(); } + if (solverChoice.do_forest_drag) { forest_drag = m_forest_drag[level]->get_drag_field(); } // Immersed Forcing MultiFab* terrain_blank = nullptr; - if(solverChoice.do_terrain) { terrain_blank = m_terrain[level]->get_terrain_blank_field(); } + if(solverChoice.do_terrain_drag) { terrain_blank = m_terrain_drag[level]->get_terrain_blank_field(); } // Moving terrain if ( solverChoice.terrain_type == TerrainType::Moving ) @@ -435,10 +435,10 @@ // Canopy data for mom sources MultiFab* forest_drag = nullptr; - if (solverChoice.do_forest) { forest_drag = m_forest[level]->get_drag_field(); } + if (solverChoice.do_forest_drag) { forest_drag = m_forest_drag[level]->get_drag_field(); } // Immersed Forcing MultiFab* terrain_blank = nullptr; - if(solverChoice.do_terrain) { terrain_blank = m_terrain[level]->get_terrain_blank_field(); } + if(solverChoice.do_terrain_drag) { terrain_blank = m_terrain_drag[level]->get_terrain_blank_field(); } make_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, cc_src, z_phys_cc[level], #if defined(ERF_USE_RRTMGP) From e9cec7237bbbf2ddad76621e8cdd5fde01b25860 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Thu, 12 Dec 2024 14:37:00 -0800 Subject: [PATCH 7/9] Saturation Adjustment Micro Model (#2020) * SatAdj micro model tested on moist bubble rise. * Add docs. --------- Co-authored-by: Ann Almgren --- CMake/BuildERFExe.cmake | 4 + Docs/sphinx_doc/Inputs.rst | 2 +- Docs/sphinx_doc/theory/Microphysics.rst | 5 + Exec/Make.ERF | 5 + Exec/MoistRegTests/Bubble/ERF_Prob.cpp | 3 - Source/DataStructs/ERF_DataStruct.H | 14 +- Source/IO/ERF_Plotfile.cpp | 2 +- .../Microphysics/ERF_EulerianMicrophysics.H | 13 +- Source/Microphysics/ERF_Microphysics.H | 6 +- Source/Microphysics/Kessler/ERF_Kessler.H | 6 +- Source/Microphysics/SAM/ERF_SAM.H | 8 +- Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp | 76 +++++++ Source/Microphysics/SatAdj/ERF_SatAdj.H | 214 ++++++++++++++++++ Source/Microphysics/SatAdj/ERF_SatAdj.cpp | 83 +++++++ .../Microphysics/SatAdj/ERF_UpdateSatAdj.cpp | 37 +++ Source/Microphysics/SatAdj/Make.package | 4 + 16 files changed, 459 insertions(+), 23 deletions(-) create mode 100644 Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp create mode 100644 Source/Microphysics/SatAdj/ERF_SatAdj.H create mode 100644 Source/Microphysics/SatAdj/ERF_SatAdj.cpp create mode 100644 Source/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp create mode 100644 Source/Microphysics/SatAdj/Make.package diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 47e1e8cab..dc4888ff7 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -207,6 +207,9 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Microphysics/Kessler/ERF_InitKessler.cpp ${SRC_DIR}/Microphysics/Kessler/ERF_Kessler.cpp ${SRC_DIR}/Microphysics/Kessler/ERF_UpdateKessler.cpp + ${SRC_DIR}/Microphysics/SatAdj/ERF_InitSatAdj.cpp + ${SRC_DIR}/Microphysics/SatAdj/ERF_SatAdj.cpp + ${SRC_DIR}/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp ${SRC_DIR}/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp ${SRC_DIR}/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp ${SRC_DIR}/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp @@ -255,6 +258,7 @@ endif() target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) + target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 46ec4148a..f654250d3 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -1443,7 +1443,7 @@ List of Parameters | | | Values | | +=============================+==========================+====================+============+ | **erf.moisture_model** | Name of moisture model | "SAM", "Kessler", | "Null" | -| | | "FastEddy" | | +| | | "SatAdj" | | +-----------------------------+--------------------------+--------------------+------------+ | **erf.do_cloud** | use basic moisture model | true / false | true | +-----------------------------+--------------------------+--------------------+------------+ diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index b62d6704a..afe563eff 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -170,3 +170,8 @@ The evaporation rate of rain is .. math:: Q_{revp} = 2\pi(S-1)n_{0R}[0.78\lambda_{R}^{-2}+0.31S_{c}^{1/3}\Gamma[(b+5)/2]a^{1/2}\mu^{-1/2}(\frac{\rho_{0}}{\rho})^{1/4}\lambda_{R}^{(b+5)/2}](\frac{1}{\rho})(\frac{L_{v}^{2}}{K_{0}R_{w}T^{2}}+\frac{1}{\rho r_{s}\psi})^{-1} + + Saturation Adjustment Microphysics Model +---------------------------------- +The saturation adjustment microphysics model is the simplest possible moisture model and only transports the +water vapor mixing ratio, :math:`q_v`, and the cloud water mixing ration, :math:`q_c`. Evaporation, :math:`q_v \longrightarrow q_c`, and condensation, :math:`q_c \longrightarrow q_v`, are the only relevant mechanisms. The final saturation state, :math:`q_v = q_{vs}(T)` is obtained from Newton-Raphson iterations on the thermal temperature :math:`T`. diff --git a/Exec/Make.ERF b/Exec/Make.ERF index d90ccc172..d25135f68 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -149,6 +149,11 @@ include $(ERF_MOISTURE_KESSLER_DIR)/Make.package VPATH_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) INCLUDE_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) +ERF_MOISTURE_SATADJ_DIR = $(ERF_SOURCE_DIR)/Microphysics/SatAdj +include $(ERF_MOISTURE_SATADJ_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_MOISTURE_SATADJ_DIR) +INCLUDE_LOCATIONS += $(ERF_MOISTURE_SATADJ_DIR) + # If using windfarm parametrization, then compile all models and choose # at runtime from the inputs ifeq ($(USE_WINDFARM), TRUE) diff --git a/Exec/MoistRegTests/Bubble/ERF_Prob.cpp b/Exec/MoistRegTests/Bubble/ERF_Prob.cpp index a28cdfdaa..f079ba5d9 100644 --- a/Exec/MoistRegTests/Bubble/ERF_Prob.cpp +++ b/Exec/MoistRegTests/Bubble/ERF_Prob.cpp @@ -303,7 +303,6 @@ Problem::init_custom_pert( if (use_moisture) { state_pert(i, j, k, RhoQ1_comp) = 0.0; state_pert(i, j, k, RhoQ2_comp) = 0.0; - state_pert(i, j, k, RhoQ3_comp) = 0.0; } }); } // do_moist_bubble @@ -389,7 +388,6 @@ Problem::init_custom_pert( // mean states state_pert(i, j, k, RhoQ1_comp) = rho*q_v_hot; state_pert(i, j, k, RhoQ2_comp) = rho*(parms_d.qt_init - q_v_hot); - state_pert(i, j, k, RhoQ3_comp) = 0.0; // Cold microphysics are present int nstate = state_pert.ncomp; @@ -429,7 +427,6 @@ Problem::init_custom_pert( if (use_moisture) { state_pert(i, j, k, RhoQ1_comp) = 0.0; state_pert(i, j, k, RhoQ2_comp) = 0.0; - state_pert(i, j, k, RhoQ3_comp) = 0.0; } }); } // do_moist_bubble diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 0be9f91e5..c7e9e6b32 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -42,7 +42,7 @@ AMREX_ENUM(MoistureModelType, ); AMREX_ENUM(MoistureType, - Kessler, SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler_NoRain, None + SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler, Kessler_NoRain, SatAdj, None ); AMREX_ENUM(WindFarmType, @@ -130,15 +130,19 @@ struct SolverChoice { } else if (moisture_type == MoistureType::Kessler_NoRain) { RhoQv_comp = RhoQ1_comp; RhoQc_comp = RhoQ2_comp; + } else if (moisture_type == MoistureType::SatAdj) { + RhoQv_comp = RhoQ1_comp; + RhoQc_comp = RhoQ2_comp; } // TODO: should we set default for dry?? // 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_NoIce || - moisture_type == MoistureType::SAM_NoPrecip_NoIce) + if (moisture_type == MoistureType::Kessler_NoRain || + moisture_type == MoistureType::SAM || + moisture_type == MoistureType::SAM_NoIce || + moisture_type == MoistureType::SAM_NoPrecip_NoIce || + moisture_type == MoistureType::SatAdj) { buoyancy_type = 1; // uses Rhoprime } else { diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index fb1d331e6..26ed76e16 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -1072,7 +1072,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p // Precipitating components //-------------------------------------------------------------------------- - if(containerHasElement(plot_var_names, "qp")) + if(containerHasElement(plot_var_names, "qp") && (n_qstate >= 3)) { int n_start = (n_qstate > 3) ? RhoQ4_comp : RhoQ3_comp; int n_end = ncomp_cons - 1; diff --git a/Source/Microphysics/ERF_EulerianMicrophysics.H b/Source/Microphysics/ERF_EulerianMicrophysics.H index cb91854b4..edf9c34cd 100644 --- a/Source/Microphysics/ERF_EulerianMicrophysics.H +++ b/Source/Microphysics/ERF_EulerianMicrophysics.H @@ -8,6 +8,7 @@ #include "ERF_NullMoist.H" #include "ERF_SAM.H" #include "ERF_Kessler.H" +#include "ERF_SatAdj.H" #include "ERF_Microphysics.H" /*! \brief Eulerian microphysics interface */ @@ -34,6 +35,8 @@ public: } else if (a_model_type == MoistureType::Kessler || a_model_type == MoistureType::Kessler_NoRain) { SetModel(); + } else if (a_model_type == MoistureType::SatAdj) { + SetModel(); } else if (a_model_type == MoistureType::None) { SetModel(); } else { @@ -108,12 +111,12 @@ public: /*! \brief get the indices and names of moisture model variables for restart at a given level */ - void Get_Qmoist_Restart_Vars ( const int a_lev, /*!< level */ - const SolverChoice& a_sc, /*!< Solver choice object */ - std::vector& a_idx, /*!< indices */ - std::vector& a_names /*!< names */ ) const override + void Get_Qmoist_Restart_Vars (const int a_lev, /*!< level */ + const SolverChoice& a_sc, /*!< Solver choice object */ + std::vector& a_idx, /*!< indices */ + std::vector& a_names /*!< names */ ) const override { - m_moist_model[a_lev]->Qmoist_Restart_Vars( a_sc, a_idx, a_names ); + m_moist_model[a_lev]->Qmoist_Restart_Vars(a_sc, a_idx, a_names); } protected: diff --git a/Source/Microphysics/ERF_Microphysics.H b/Source/Microphysics/ERF_Microphysics.H index 29eab027d..084572e0f 100644 --- a/Source/Microphysics/ERF_Microphysics.H +++ b/Source/Microphysics/ERF_Microphysics.H @@ -58,7 +58,10 @@ public: /*! \brief get the indices and names of moisture model variables for restart at a given level */ - virtual void Get_Qmoist_Restart_Vars ( int, const SolverChoice&, std::vector&, std::vector& ) const = 0; + virtual void Get_Qmoist_Restart_Vars (int, + const SolverChoice&, + std::vector&, + std::vector& ) const = 0; /*! \brief query if a specified moisture model is Eulerian or Lagrangian */ static MoistureModelType modelType (const MoistureType a_moisture_type) @@ -68,6 +71,7 @@ public: || (a_moisture_type == MoistureType::SAM_NoPrecip_NoIce) || (a_moisture_type == MoistureType::Kessler) || (a_moisture_type == MoistureType::Kessler_NoRain) + || (a_moisture_type == MoistureType::SatAdj) || (a_moisture_type == MoistureType::None) ) { return MoistureModelType::Eulerian; } else { diff --git a/Source/Microphysics/Kessler/ERF_Kessler.H b/Source/Microphysics/Kessler/ERF_Kessler.H index f361f013a..f890ef35b 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.H +++ b/Source/Microphysics/Kessler/ERF_Kessler.H @@ -125,9 +125,9 @@ public: Qstate_Size () override { return Kessler::m_qstate_size; } void - Qmoist_Restart_Vars ( const SolverChoice& a_sc, - std::vector& a_idx, - std::vector& a_names) const override + Qmoist_Restart_Vars (const SolverChoice& a_sc, + std::vector& a_idx, + std::vector& a_names) const override { a_idx.clear(); a_names.clear(); diff --git a/Source/Microphysics/SAM/ERF_SAM.H b/Source/Microphysics/SAM/ERF_SAM.H index e514b509a..7652b3093 100644 --- a/Source/Microphysics/SAM/ERF_SAM.H +++ b/Source/Microphysics/SAM/ERF_SAM.H @@ -264,9 +264,9 @@ public: } void - Qmoist_Restart_Vars ( const SolverChoice& a_sc, - std::vector& a_idx, - std::vector& a_names) const override + Qmoist_Restart_Vars (const SolverChoice& a_sc, + std::vector& a_idx, + std::vector& a_names) const override { a_idx.clear(); a_names.clear(); @@ -283,7 +283,7 @@ private: // Number of qmoist variables (qt, qv, qcl, qci, qp, qpr, qps, qpg) int m_qmoist_size = 11; - // Number of qmoist variables + // Number of qstate variables int m_qstate_size = 6; // CFL MAX for vertical advection diff --git a/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp b/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp new file mode 100644 index 000000000..5ce0b67c9 --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp @@ -0,0 +1,76 @@ +#include "ERF_SatAdj.H" + +using namespace amrex; + + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + * @param[in] qc_in Cloud variables input + * @param[in,out] qv_in Vapor variables input + * @param[in] qi_in Ice variables input + * @param[in] grids The boxes on which we will evolve the solution + * @param[in] geom Geometry associated with these MultiFabs and grids + * @param[in] dt_advance Timestep for the advance + */ +void SatAdj::Init (const MultiFab& cons_in, + const BoxArray& /*grids*/, + const Geometry& geom, + const Real& dt_advance, + std::unique_ptr& /*z_phys_nd*/, + std::unique_ptr& /*detJ_cc*/) +{ + dt = dt_advance; + m_geom = geom; + + MicVarMap.resize(m_qmoist_size); + MicVarMap = {MicVar_SatAdj::qv, MicVar_SatAdj::qc}; + + // initialize microphysics variables + for (auto ivar = 0; ivar < MicVar_SatAdj::NumVars; ++ivar) { + mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), + 1, cons_in.nGrowVect()); + mic_fab_vars[ivar]->setVal(0.); + } +} + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + */ +void SatAdj::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) { + const auto& tbx = mfi.tilebox(); + + auto states_array = cons_in.array(mfi); + + auto qv_array = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi); + + auto rho_array = mic_fab_vars[MicVar_SatAdj::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar_SatAdj::tabs]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_SatAdj::pres]->array(mfi); + + // Get pressure, theta, temperature, density, and qt, qp + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + rho_array(i,j,k) = states_array(i,j,k,Rho_comp); + theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); + qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); + qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); + + 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)); + + // Pressure in [mbar] for qsat evaluation + 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/SatAdj/ERF_SatAdj.H b/Source/Microphysics/SatAdj/ERF_SatAdj.H new file mode 100644 index 000000000..8b42743b5 --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.H @@ -0,0 +1,214 @@ +#ifndef ERF_SATADJ_H +#define ERF_SATADJ_H + +/* + * Implementation saturation adjustment microphysics model + * The model transports qv and qc and does Newton iterations + * to complete condensation/evaporation. + */ + +#include +#include +#include + +#include +#include +#include +#include + +#include "ERF_EOS.H" +#include "ERF_Constants.H" +#include "ERF_MicrophysicsUtils.H" +#include "ERF_IndexDefines.H" +#include "ERF_DataStruct.H" +#include "ERF_NullMoist.H" +#include "ERF_TileNoZ.H" + +namespace MicVar_SatAdj { + enum { + // independent variables + rho=0, // density + theta, // liquid/ice water potential temperature + tabs, // temperature + pres, // pressure + // non-precipitating vars + qv, // cloud vapor + qc, // cloud water + NumVars + }; +} + +class SatAdj : public NullMoist { + + using FabPtr = std::shared_ptr; + +public: + // constructor + SatAdj () {} + + // destructor + virtual ~SatAdj () = default; + + // cloud physics + void AdvanceSatAdj (const SolverChoice& /*solverChoice*/); + + // Set up for first time + void + Define (SolverChoice& sc) override + { + m_fac_cond = lcond / sc.c_p; + m_rdOcp = sc.rdOcp; + } + + // init + void + Init (const amrex::MultiFab& cons_in, + const amrex::BoxArray& /*grids*/, + const amrex::Geometry& geom, + const amrex::Real& dt_advance, + std::unique_ptr& /*z_phys_nd*/, + std::unique_ptr& /*detJ_cc*/) override; + + // Copy state into micro vars + void + Copy_State_to_Micro (const amrex::MultiFab& cons_in) override; + + // Copy state into micro vars + void + Copy_Micro_to_State (amrex::MultiFab& cons_in) override; + + // update micro vars + void + Update_Micro_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_State_to_Micro(cons_in); + } + + // update state vars + void + Update_State_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_Micro_to_State(cons_in); + } + + // wrapper to advance micro vars + void + Advance (const amrex::Real& dt_advance, + const SolverChoice& solverChoice) override + { + dt = dt_advance; + + this->AdvanceSatAdj(solverChoice); + } + + amrex::MultiFab* + Qmoist_Ptr (const int& varIdx) override + { + AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size); + return mic_fab_vars[MicVarMap[varIdx]].get(); + } + + int + Qmoist_Size () override { return SatAdj::m_qmoist_size; } + + int + Qstate_Size () override { return SatAdj::m_qstate_size; } + + void + Qmoist_Restart_Vars (const SolverChoice& /*a_sc*/, + std::vector& a_idx, + std::vector& a_names) const override + { + a_idx.clear(); + a_names.clear(); + } + + AMREX_GPU_HOST_DEVICE + AMREX_FORCE_INLINE + static amrex::Real + NewtonIterSat (int& i, int& j, int& k, + const amrex::Real& fac_cond, + const amrex::Array4& tabs_array, + const amrex::Array4& pres_array, + const amrex::Array4& qv_array, + const amrex::Array4& qc_array) + { + // Solution tolerance + amrex::Real tol = 1.0e-8; + + // Saturation moisture fractions + amrex::Real qsat, dqsat; + + // Newton iteration vars + int niter; + amrex::Real fff, dfff; + amrex::Real dtabs, delta_qv; + + // 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 { + // Saturation moisture fractions + erf_qsatw(tabs, pres, qsat); + erf_dtqsatw(tabs, pres, dqsat); + + // Function for root finding: + // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat) + fff = -tabs + tabs_array(i,j,k) + fac_cond*(qv_array(i,j,k) - qsat); + + // Derivative of function (T_new iterated on) + dfff = -1.0 - fac_cond*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; + + // Partition the change in non-precipitating q + qv_array(i,j,k) = qsat; + qc_array(i,j,k) += delta_qv; + + // Return to temperature + return tabs; + } + +private: + // Number of qmoist variables (qv, qc) + int m_qmoist_size = 2; + + // Number of qstate variables + int m_qstate_size = 2; + + // MicVar map (Qmoist indices -> MicVar enum) + amrex::Vector MicVarMap; + + // geometry + amrex::Geometry m_geom; + + // timestep + amrex::Real dt; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_rdOcp ; + + // independent variables + amrex::Array mic_fab_vars; +}; +#endif diff --git a/Source/Microphysics/SatAdj/ERF_SatAdj.cpp b/Source/Microphysics/SatAdj/ERF_SatAdj.cpp new file mode 100644 index 000000000..0956e4156 --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.cpp @@ -0,0 +1,83 @@ +#include "ERF_SatAdj.H" + +using namespace amrex; + +/** + * Compute Precipitation-related Microphysics quantities. + */ +void SatAdj::AdvanceSatAdj (const SolverChoice& /*solverChoice*/) +{ + auto tabs = mic_fab_vars[MicVar_SatAdj::tabs]; + + // Expose for GPU + Real d_fac_cond = m_fac_cond; + Real rdOcp = m_rdOcp; + + // get the temperature, dentisy, theta, qt and qc from input + for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const auto& tbx = mfi.tilebox(); + + auto qv_array = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar_SatAdj::tabs]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_SatAdj::pres]->array(mfi); + + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); + + //------- Evaporation/condensation + Real qsat; + erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsat); + + // There is enough moisutre to drive to equilibrium + if ((qv_array(i,j,k)+qc_array(i,j,k)) > qsat) { + + // Update temperature + tabs_array(i,j,k) = NewtonIterSat(i, j, k , + d_fac_cond, tabs_array, pres_array, + qv_array , qc_array ); + + // Update theta (constant pressure) + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); + + // + // 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 + Real delta_qc = qc_array(i,j,k); + + // Partition the change in non-precipitating q + qv_array(i,j,k) += qc_array(i,j,k); + qc_array(i,j,k) = 0.0; + + // Update temperature (endothermic since we evap/sublime) + tabs_array(i,j,k) -= d_fac_cond * delta_qc; + + // Update theta + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); + + // Verify assumption that qv > qsat does not occur + erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsat); + if (qv_array(i,j,k) > qsat) { + + // Update temperature + tabs_array(i,j,k) = NewtonIterSat(i, j, k , + d_fac_cond , tabs_array, pres_array, + qv_array , qc_array ); + + // 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/SatAdj/ERF_UpdateSatAdj.cpp b/Source/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp new file mode 100644 index 000000000..b45d619ef --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp @@ -0,0 +1,37 @@ +#include "ERF_SatAdj.H" + +using namespace amrex; + +/** + * Updates conserved and microphysics variables in the provided MultiFabs from + * the internal MultiFabs that store Microphysics module data. + * + * @param[out] cons Conserved variables + * @param[out] qmoist: qv, qc + */ +void SatAdj::Copy_Micro_to_State (MultiFab& cons) +{ + // Get the temperature, density, theta, qt and qp from input + for (MFIter mfi(cons,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const auto& tbx = mfi.tilebox(); + + auto states_arr = cons.array(mfi); + + auto rho_arr = mic_fab_vars[MicVar_SatAdj::rho]->array(mfi); + auto theta_arr = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi); + auto qv_arr = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi); + auto qc_arr = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi); + + // get potential total density, temperature, qt, qp + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); + states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k); + states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k); + }); + } + + // Fill interior ghost cells and periodic boundaries + cons.FillBoundary(m_geom.periodicity()); +} + diff --git a/Source/Microphysics/SatAdj/Make.package b/Source/Microphysics/SatAdj/Make.package new file mode 100644 index 000000000..5b9ff43e4 --- /dev/null +++ b/Source/Microphysics/SatAdj/Make.package @@ -0,0 +1,4 @@ +CEXE_sources += ERF_InitSatAdj.cpp +CEXE_sources += ERF_UpdateSatAdj.cpp +CEXE_sources += ERF_SatAdj.cpp +CEXE_headers += ERF_SatAdj.H From 7edaee3529e96203491aeabf81c10e7d3dfe6ea8 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Thu, 12 Dec 2024 16:23:40 -0800 Subject: [PATCH 8/9] fix for GPU (#2023) --- Source/LinearSolvers/ERF_TerrainPoisson.cpp | 78 +++++++++++++-------- 1 file changed, 50 insertions(+), 28 deletions(-) diff --git a/Source/LinearSolvers/ERF_TerrainPoisson.cpp b/Source/LinearSolvers/ERF_TerrainPoisson.cpp index 7eea971f7..f7bf18849 100644 --- a/Source/LinearSolvers/ERF_TerrainPoisson.cpp +++ b/Source/LinearSolvers/ERF_TerrainPoisson.cpp @@ -60,48 +60,70 @@ void TerrainPoisson::apply_bcs (MultiFab& phi) { Box bx = mfi.tilebox(); const Array4& phi_arr = phi.array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - if (i == domlo.x) { - if (bc_fft[0].first == FFT::Boundary::even) { + if (bx.smallEnd(0) <= domlo.x) { + if (bc_fft[0].first == FFT::Boundary::even) { + ParallelFor(makeSlab(bx,0,domlo.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { phi_arr(i-1,j,k) = phi_arr(i,j,k); - } else if (bc_fft[0].first == FFT::Boundary::odd) { - phi_arr(i-1,j,k) = -phi_arr(i,j,k); - } - } else if (i == domhi.x) { - if (bc_fft[0].second == FFT::Boundary::even) { + }); + } else if (bc_fft[0].first == FFT::Boundary::odd) { + ParallelFor(makeSlab(bx,0,domlo.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + phi_arr(i-1,j,k) = -phi_arr(i,j,k); + }); + } + } // lo x + if (bx.bigEnd(0) >= domhi.x) { + if (bc_fft[0].second == FFT::Boundary::even) { + ParallelFor(makeSlab(bx,0,domhi.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { phi_arr(i+1,j,k) = phi_arr(i,j,k); - } else if (bc_fft[0].second == FFT::Boundary::odd) { - phi_arr(i+1,j,k) = -phi_arr(i,j,k); - } + }); + } else if (bc_fft[0].second == FFT::Boundary::odd) { + ParallelFor(makeSlab(bx,0,domhi.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + phi_arr(i+1,j,k) = -phi_arr(i,j,k); + }); } - }); + } // lo x } // mfi - } + } // not periodic in x + if (!m_geom.isPeriodic(1)) { for (MFIter mfi(phi,true); mfi.isValid(); ++mfi) { Box bx = mfi.tilebox(); Box bx2(bx); bx2.grow(0,1); const Array4& phi_arr = phi.array(mfi); - ParallelFor(bx2, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - if (j == domlo.y) { - if (bc_fft[1].first == FFT::Boundary::even) { + if (bx.smallEnd(1) <= domlo.y) { + if (bc_fft[1].first == FFT::Boundary::even) { + ParallelFor(makeSlab(bx2,1,domlo.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { phi_arr(i,j-1,k) = phi_arr(i,j,k); - } else if (bc_fft[1].first == FFT::Boundary::odd) { - phi_arr(i,j-1,k) = -phi_arr(i,j,k); - } - } else if (j == domhi.y) { - if (bc_fft[1].second == FFT::Boundary::even) { + }); + } else if (bc_fft[1].first == FFT::Boundary::odd) { + ParallelFor(makeSlab(bx2,1,domlo.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + phi_arr(i,j-1,k) = -phi_arr(i,j,k); + }); + } + } // lo x + if (bx.bigEnd(1) >= domhi.y) { + if (bc_fft[1].second == FFT::Boundary::even) { + ParallelFor(makeSlab(bx2,1,domhi.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { phi_arr(i,j+1,k) = phi_arr(i,j,k); - } else if (bc_fft[1].second == FFT::Boundary::odd) { - phi_arr(i,j+1,k) = -phi_arr(i,j,k); - } + }); + } else if (bc_fft[1].second == FFT::Boundary::odd) { + ParallelFor(makeSlab(bx2,1,domhi.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + phi_arr(i,j+1,k) = -phi_arr(i,j,k); + }); } - }); + } // lo x + } // mfi - } + } // not periodic in y for (MFIter mfi(phi,true); mfi.isValid(); ++mfi) { From f4beeafc10d63ac1bd2bf284bfca4982040fd81f Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Fri, 13 Dec 2024 09:36:28 -0700 Subject: [PATCH 9/9] Add erf.substepping_cfl param (#2024) --- Docs/sphinx_doc/Inputs.rst | 10 ++++++---- Source/ERF.H | 1 + Source/ERF.cpp | 2 ++ Source/TimeIntegration/ERF_ComputeTimestep.cpp | 5 +++-- 4 files changed, 12 insertions(+), 6 deletions(-) diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index f654250d3..abf362ad9 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -432,10 +432,12 @@ List of Parameters | **erf.no_substepping** | Should we turn off | int (0 or 1) | 0 | | | substepping in time? | | | +----------------------------+----------------------+----------------+-------------------+ -| **erf.cfl** | CFL number for | Real > 0 and | 0.8 | -| | hydro | <= 1 | | -| | | | | -| | | | | +| **erf.cfl** | CFL number used to | Real > 0 and | 0.8 | +| | compute level 0 dt | <= 1 | | ++----------------------------+----------------------+----------------+-------------------+ +| **erf.substepping_cfl** | CFL number used to | Real > 0 and | 1.0 | +| | compute the number | <= 1 | | +| | of substeps | | | +----------------------------+----------------------+----------------+-------------------+ | **erf.fixed_dt** | set level 0 dt | Real > 0 | unused if not | | | as this value | | set | diff --git a/Source/ERF.H b/Source/ERF.H index 06370d69b..27dd95b14 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -922,6 +922,7 @@ private: // Time step controls static amrex::Real cfl; + static amrex::Real sub_cfl; static amrex::Real init_shrink; static amrex::Real change_max; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 1c29ccaad..3759ede46 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -26,6 +26,7 @@ SolverChoice ERF::solverChoice; // Time step control Real ERF::cfl = 0.8; +Real ERF::sub_cfl = 1.0; Real ERF::init_shrink = 1.0; Real ERF::change_max = 1.1; int ERF::fixed_mri_dt_ratio = 0; @@ -1455,6 +1456,7 @@ ERF::ReadParameters () // Time step controls pp.query("cfl", cfl); + pp.query("substepping_cfl", sub_cfl); pp.query("init_shrink", init_shrink); pp.query("change_max", change_max); diff --git a/Source/TimeIntegration/ERF_ComputeTimestep.cpp b/Source/TimeIntegration/ERF_ComputeTimestep.cpp index 6e91a1661..5c79fb85a 100644 --- a/Source/TimeIntegration/ERF_ComputeTimestep.cpp +++ b/Source/TimeIntegration/ERF_ComputeTimestep.cpp @@ -202,8 +202,9 @@ ERF::estTimeStep (int level, long& dt_fast_ratio) const if (fixed_dt[level] > 0. && fixed_fast_dt[level] > 0.) { dt_fast_ratio = static_cast( fixed_dt[level] / fixed_fast_dt[level] ); } else if (fixed_dt[level] > 0.) { - // Max CFL_c = 1.0 for substeps, but we enforce a min of 4 substeps - auto dt_sub_max = (estdt_comp/cfl); + // Max CFL_c = 1.0 for substeps by default, but we enforce a min of 4 substeps + auto dt_sub_max = (estdt_comp/cfl * sub_cfl); + Print() << "fixed_dt="<