Skip to content

Commit

Permalink
start to replace terrain_type by mesh_type (#2018)
Browse files Browse the repository at this point in the history 
* start to replace terrain_type by mesh_type

* if we read in terrain_type = Static, set mesh_type to VariableDz

* update for FFT
  • Loading branch information
@asalmgren
asalmgren authored Dec 12, 2024
1 parent 99a6690 commit 6997ffc
Show file tree
Hide file tree
Showing 19 changed files with 112 additions and 81 deletions.
35 changes: 33 additions & 2 deletions Source/DataStructs/ERF_DataStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,10 @@ AMREX_ENUM(SubsteppingType,
None, Explicit, Implicit
);

AMREX_ENUM(MeshType,
ConstantDz, StretchedDz, VariableDz
);

AMREX_ENUM(TerrainType,
None, Static, Moving
);
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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?
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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 ()
{
Expand Down
5 changes: 4 additions & 1 deletion Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -826,8 +826,10 @@ private:
amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_q1fx1_lev, SFS_q1fx2_lev, SFS_q1fx3_lev;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_q2fx3_lev;

// Terrain / grid stretching
// Grid stretching
amrex::Vector<amrex::Vector<amrex::Real>> zlevels_stag; // nominal height levels

// Terrain data structures
amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_phys_nd;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_phys_cc;

Expand All @@ -850,6 +852,7 @@ private:

amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_t_rk;

// Map scale factors
amrex::Vector<std::unique_ptr<amrex::MultiFab>> mapfac_m;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> mapfac_u;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> mapfac_v;
Expand Down
39 changes: 19 additions & 20 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand All @@ -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)
Expand Down Expand Up @@ -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--)
{
Expand All @@ -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<const Real> mapfac_arr = mapfac_m[lev]->const_array(mfi);
if (use_terrain) {
const Array4<const Real> 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<const Real> 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
Expand All @@ -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<const Real> mapfac_arr = mapfac_m[lev]->const_array(mfi);
if (use_terrain) {
const Array4<const Real> 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<const Real> 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
Expand Down Expand Up @@ -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));
Expand Down Expand Up @@ -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;
}
}

Expand Down Expand Up @@ -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],
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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);
}
}
Expand Down
13 changes: 8 additions & 5 deletions Source/ERF_MakeNewArrays.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<MultiFab>(ba,dm,1,1);

if (solverChoice.terrain_type == TerrainType::Moving)
Expand Down Expand Up @@ -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)
{
Expand Down Expand Up @@ -490,7 +492,7 @@ ERF::init_zphys (int lev, Real time)
void
ERF::remake_zphys (int lev, std::unique_ptr<MultiFab>& temp_zphys_nd)
{
if (lev > 0 && SolverChoice::terrain_type != TerrainType::None)
if (lev > 0 && SolverChoice::mesh_type == MeshType::VariableDz)
{
//
// First interpolate from coarser level
Expand All @@ -517,7 +519,8 @@ ERF::remake_zphys (int lev, std::unique_ptr<MultiFab>& 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]);
Expand Down Expand Up @@ -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);
}

Expand Down
4 changes: 2 additions & 2 deletions Source/ERF_MakeNewLevel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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));
Expand Down Expand Up @@ -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));
Expand Down
2 changes: 1 addition & 1 deletion Source/ERF_Tagging.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -161,7 +161,7 @@ ERF::refinement_criteria_setup ()
jlo = static_cast<int>((rbox_lo[1] - plo[1])/dx[1]);
ihi = static_cast<int>((rbox_hi[0] - plo[0])/dx[0]-1);
jhi = static_cast<int>((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();
Expand Down
4 changes: 2 additions & 2 deletions Source/IO/ERF_Checkpoint.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
19 changes: 10 additions & 9 deletions Source/IO/ERF_Plotfile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,8 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
//
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" &&
Expand Down Expand Up @@ -225,7 +226,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p

// Vector of MultiFabs for nodal data
Vector<MultiFab> 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);
Expand Down Expand Up @@ -578,7 +579,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
const Array4<Real>& derdat = mf[lev].array(mfi);
const Array4<Real> & p_arr = pres.array(mfi);

if (SolverChoice::terrain_type != TerrainType::None) {
if (SolverChoice::mesh_type != MeshType::ConstantDz) {
const Array4<Real const>& z_nd = z_phys_nd[lev]->const_array(mfi);

ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
Expand Down Expand Up @@ -672,7 +673,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
const Array4<Real>& derdat = mf[lev].array(mfi);
const Array4<Real> & p_arr = pres.array(mfi);

if (SolverChoice::terrain_type != TerrainType::None) {
if (SolverChoice::mesh_type != MeshType::ConstantDz) {
const Array4<Real const>& z_nd = z_phys_nd[lev]->const_array(mfi);

ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
Expand Down Expand Up @@ -843,7 +844,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> 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);
Expand Down Expand Up @@ -1350,7 +1351,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> 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];
Expand Down Expand Up @@ -1399,7 +1400,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> 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),
Expand Down Expand Up @@ -1510,7 +1511,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> 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),
Expand All @@ -1523,7 +1524,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
}

} else {
if (SolverChoice::terrain_type != TerrainType::None) {
if (SolverChoice::mesh_type != MeshType::ConstantDz) {
WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1,
GetVecOfConstPtrs(mf),
GetVecOfConstPtrs(mf_nd),
Expand Down
Loading

0 comments on commit 6997ffc

Please sign in to comment.