From 6997ffc5c62cf83b107d3c6cb4beb7f99b8ebf19 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Thu, 12 Dec 2024 13:06:33 -0800 Subject: [PATCH] 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 {