diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index a972f7b68..672ef6a4c 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -164,6 +164,61 @@ Examples of Usage level-1 grids will be created every 2 level-0 time steps, and new level-2 grids will be created every 2 level-1 time steps. + +Grid Stretching +=============== + +This automatically activates **erf.use_terrain**. By default, the +problem-specific terrain is initialized to be flat at an elevation of z=0. +These inputs are used to automatically generate the staggered z levels used to +calculate the grid metric transformation. Alternatively, arbitrary z levels may +be specified with the **erf.terrain_z_levels** parameter. + +.. _list-of-parameters-3: + +List of Parameters +------------------ + ++-------------------------------+-----------------+-----------------+-------------+ +| Parameter | Definition | Acceptable | Default | +| | | Values | | ++===============================+=================+=================+=============+ +| **erf.grid_stretching_ratio** | scaling factor | Real > 0 | 0 (no grid | +| | applied to | | stretching) | +| | delta z at each | | | +| | level | | | ++-------------------------------+-----------------+-----------------+-------------+ +| **erf.initial_dz** | vertical grid | Real > 0 | must be set | +| | spacing for the | | if grid | +| | cell above the | | stretching | +| | bottom surface | | ratio is set| ++-------------------------------+-----------------+-----------------+-------------+ +| **erf.terrain_z_levels** | nominal | List of Real | NONE | +| | staggered | | | +| | z levels | | | ++-------------------------------+-----------------+-----------------+-------------+ + +.. _notes-3: + +Notes +----- + +- If both **erf.terrain_z_levels** and **erf.grid_stretching_ratio** are + specified, the simple grid stretching will be ignored. +- The number of input **erf.terrain_z_levels** must be equal **amr.n_cell** in + the z direction + 1. + +.. _examples-of-usage-3: + +Examples of Usage +----------------- + +- **erf.grid_stretching_ratio** = 1.025 + +- | **erf.initial_dz** = 5.0 + | the first cell center would be at z=2.5 + + Regridding ========== diff --git a/Exec/ABL/inputs_most_pbl b/Exec/ABL/inputs_most_pbl new file mode 100644 index 000000000..151c1e8fe --- /dev/null +++ b/Exec/ABL/inputs_most_pbl @@ -0,0 +1,65 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 3000 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 1024 1024 1024 +amr.n_cell = 4 4 64 + +#erf.grid_stretching_ratio = 1.025 +#erf.initial_dz = 16.0 + +geometry.is_periodic = 1 1 0 + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +zlo.type = "Most" +erf.most.z0 = 0.1 +erf.most.zref = 8.0 + +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 0.1 # fixed time step depending on grid resolution + +# 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 = -1 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta + +# 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.pbl_type = "MYNN2.5" + +erf.init_type = "uniform" + +# 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 + diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index b52ea0556..6a2abcb1d 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -37,6 +37,18 @@ struct SolverChoice { // Do we have terrain (or grid stretching)? pp.query("use_terrain", use_terrain); + pp.query("grid_stretching_ratio", grid_stretching_ratio); + AMREX_ASSERT_WITH_MESSAGE((grid_stretching_ratio >= 0.), + "The grid stretching ratio must be greater than 0"); + if (grid_stretching_ratio > 0) { + if (!use_terrain) { + amrex::Print() << "Turning terrain on to enable grid stretching" << std::endl; + use_terrain = true; + } + pp.query("zsurface", zsurf); + pp.get("initial_dz", dz0); + } + // Do we set map scale factors to 0.5 instead of 1 for testing? pp.query("test_mapfactor", test_mapfactor); @@ -266,6 +278,11 @@ struct SolverChoice { amrex::Real c_p = Cp_d; // specific heat at constant pressure for dry air [J/(kg-K)] amrex::Real rdOcp; + // Staggered z levels for vertical grid stretching + amrex::Real grid_stretching_ratio = 0; + amrex::Real zsurf = 0.0; + amrex::Real dz0; + #if defined(ERF_USE_WARM_NO_PRECIP) amrex::Real tau_cond = 1.0; // Default time of 1 sec -- this is somewhat arbitray #endif diff --git a/Source/ERF.H b/Source/ERF.H index ca240c279..456a7a8f5 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -550,6 +550,8 @@ private: amrex::Vector> SFS_hfx1_lev, SFS_hfx2_lev, SFS_hfx3_lev; amrex::Vector> SFS_diss_lev; + // Terrain / grid stretching + amrex::Vector zlevels_stag; amrex::Vector> z_phys_nd; amrex::Vector> z_phys_cc; amrex::Vector> detJ_cc; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 2e54eded7..7dc8b4440 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -122,6 +122,14 @@ ERF::ERF () // Geometry on all levels has been defined already. + if (solverChoice.use_terrain) { + init_zlevels(zlevels_stag, + geom[0], + solverChoice.grid_stretching_ratio, + solverChoice.zsurf, + solverChoice.dz0); + } + // No valid BoxArray and DistributionMapping have been defined. // But the arrays for them have been resized. diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index bfe5d0667..bb03cae21 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -410,7 +410,7 @@ ERF::update_terrain_arrays (int lev, Real time) if (solverChoice.use_terrain) { if (init_type != "real" && init_type != "metgrid") { prob->init_custom_terrain(geom[lev],*z_phys_nd[lev],time); - init_terrain_grid(geom[lev],*z_phys_nd[lev]); + init_terrain_grid(geom[lev],*z_phys_nd[lev],zlevels_stag); } if (lev>0 && (init_type == "real" || init_type == "metgrid")) { PhysBCFunctNoOp empty_bc; diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index d9d037115..f410aa4bb 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -103,7 +103,7 @@ ERF::init_from_metgrid(int lev) } // mf // This defines all the z(i,j,k) values given z(i,j,0) from above. - init_terrain_grid(geom[lev], *z_phys); + init_terrain_grid(geom[lev], *z_phys, zlevels_stag); // This makes the Jacobian. make_J (geom[lev],*z_phys, *detJ_cc[lev]); diff --git a/Source/TimeIntegration/TI_fast_rhs_fun.H b/Source/TimeIntegration/TI_fast_rhs_fun.H index 9d28b4dde..b5658e2c4 100644 --- a/Source/TimeIntegration/TI_fast_rhs_fun.H +++ b/Source/TimeIntegration/TI_fast_rhs_fun.H @@ -29,13 +29,13 @@ auto fast_rhs_fun = [&](int fast_step, int n_sub, // Make "old" fast geom -- store in z_phys_nd for convenience if (verbose) Print() << "Making geometry at start of substep time: " << old_substep_time << std::endl; prob->init_custom_terrain(fine_geom,*z_phys_nd[level],old_substep_time); - init_terrain_grid (fine_geom,*z_phys_nd[level]); + init_terrain_grid (fine_geom,*z_phys_nd[level], zlevels_stag); make_J (fine_geom,*z_phys_nd[level], *detJ_cc[level]); // Make "new" fast geom if (verbose) Print() << "Making geometry for end of substep time :" << new_substep_time << std::endl; prob->init_custom_terrain(fine_geom,*z_phys_nd_new[level],new_substep_time); - init_terrain_grid (fine_geom,*z_phys_nd_new[level]); + init_terrain_grid (fine_geom,*z_phys_nd_new[level], zlevels_stag); make_J (fine_geom,*z_phys_nd_new[level], *detJ_cc_new[level]); Real inv_dt = 1./dtau; diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 08fc4ce2e..f30077ddc 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -29,17 +29,17 @@ if (verbose) Print() << "Re-making old geometry at old time : " << old_step_time << std::endl; prob->init_custom_terrain(fine_geom,*z_phys_nd[level],old_step_time); - init_terrain_grid (fine_geom,*z_phys_nd[level]); + init_terrain_grid (fine_geom,*z_phys_nd[level], zlevels_stag); make_J (fine_geom,*z_phys_nd[level], *detJ_cc[level]); if (verbose) Print() << "Making src geometry at old_stage_time: " << old_stage_time << std::endl; prob->init_custom_terrain(fine_geom,*z_phys_nd_src[level],old_stage_time); - init_terrain_grid (fine_geom,*z_phys_nd_src[level]); + init_terrain_grid (fine_geom,*z_phys_nd_src[level], zlevels_stag); make_J (fine_geom,*z_phys_nd_src[level], *detJ_cc_src[level]); if (verbose) Print() << "Making new geometry at new_stage_time: " << new_stage_time << std::endl; prob->init_custom_terrain(fine_geom,*z_phys_nd_new[level],new_stage_time); - init_terrain_grid (fine_geom,*z_phys_nd_new[level]); + init_terrain_grid (fine_geom,*z_phys_nd_new[level], zlevels_stag); make_J (fine_geom,*z_phys_nd_new[level], *detJ_cc_new[level]); Real inv_dt = 1./slow_dt; diff --git a/Source/Utils/TerrainMetrics.H b/Source/Utils/TerrainMetrics.H index a95c19278..0d8faa9af 100644 --- a/Source/Utils/TerrainMetrics.H +++ b/Source/Utils/TerrainMetrics.H @@ -10,8 +10,16 @@ * Utility routines for constructing terrain metric terms */ -// Declare function for ERF.cpp -void init_terrain_grid( const amrex::Geometry& geom, amrex::MultiFab& z_phys_nd); +// Declare functions for ERF.cpp +void init_zlevels (amrex::Vector& zlevels_stag, + const amrex::Geometry& geom, + const amrex::Real grid_stretching_ratio, + const amrex::Real zsurf, + const amrex::Real dz0); + +void init_terrain_grid (const amrex::Geometry& geom, + amrex::MultiFab& z_phys_nd, + amrex::Vector const& z_levels_h); //***************************************************************************************** // Compute terrain metric terms at cell-center diff --git a/Source/Utils/TerrainMetrics.cpp b/Source/Utils/TerrainMetrics.cpp index be9bfe89b..4fd2ed9ba 100644 --- a/Source/Utils/TerrainMetrics.cpp +++ b/Source/Utils/TerrainMetrics.cpp @@ -5,6 +5,59 @@ using namespace amrex; +void +init_zlevels (amrex::Vector& zlevels_stag, + const amrex::Geometry& geom, + const amrex::Real grid_stretching_ratio, + const amrex::Real zsurf, + const amrex::Real dz0) +{ + auto dx = geom.CellSizeArray(); + const amrex::Box& domain = geom.Domain(); + int nz = domain.length(2)+1; // staggered + zlevels_stag.resize(nz); + + if (grid_stretching_ratio == 0) { + // This is the default for z_levels + for (int k = 0; k < nz; k++) + { + zlevels_stag[k] = k * dx[2]; + } + } else { + // Create stretched grid based on initial dz and stretching ratio + zlevels_stag[0] = zsurf; + amrex::Real dz = dz0; + amrex::Print() << "Stretched grid levels: " << zsurf; + for (int k = 1; k < nz; k++) + { + zlevels_stag[k] = zlevels_stag[k-1] + dz; + amrex::Print() << " " << zlevels_stag[k]; + dz *= grid_stretching_ratio; + } + amrex::Print() << std::endl; + } + + // Try reading in terrain_z_levels, which allows arbitrarily spaced grid + // levels to be specified and will take precedence over the + // grid_stretching_ratio parameter + ParmParse pp("erf"); + int n_zlevels = pp.countval("terrain_z_levels"); + if (n_zlevels > 0) + { + if (n_zlevels != nz) { + amrex::Print() << "You supplied " << n_zlevels << " staggered terrain_z_levels " << std::endl; + amrex::Print() << "but n_cell+1 in the z-direction is " << nz << std::endl; + amrex::Abort("You must specify a z_level for every value of k"); + } + + if (grid_stretching_ratio > 0) { + amrex::Print() << "Note: Found terrain_z_levels, ignoring grid_stretching_ratio" << std::endl; + } + + pp.queryarr("terrain_z_levels", zlevels_stag, 0, nz); + } +} + /** * Computation of the terrain grid from BTF, STF, or Sullivan TF model * @@ -14,7 +67,7 @@ using namespace amrex; */ void -init_terrain_grid (const Geometry& geom, MultiFab& z_phys_nd) +init_terrain_grid (const Geometry& geom, MultiFab& z_phys_nd, amrex::Vector const& z_levels_h) { auto dx = geom.CellSizeArray(); auto ProbHiArr = geom.ProbHiArray(); @@ -24,6 +77,7 @@ init_terrain_grid (const Geometry& geom, MultiFab& z_phys_nd) int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1; int domlo_z = domain.smallEnd(2); int domhi_z = domain.bigEnd(2) + 1; + int nz = domain.length(2)+1; // User-selected method from inputs file (BTF default) ParmParse pp("erf"); @@ -31,28 +85,7 @@ init_terrain_grid (const Geometry& geom, MultiFab& z_phys_nd) pp.query("terrain_smoothing", terrain_smoothing); // "custom" refers to terrain analytically specified, such as WitchOfAgnesi - std::string terrain_type = "custom"; - - int nz = domain.length(2)+1; - amrex::Vector z_levels_h; - z_levels_h.resize(nz); - - // This is the default for z_levels - for (int k = 0; k < nz; k++) - { - z_levels_h[k] = k * dx[2]; - } - - // But we can read them in from the inputs file as an alternative - int n_zlevels = pp.countval("terrain_z_levels"); - - pp.queryarr("terrain_z_levels", z_levels_h, 0, nz); - - if (n_zlevels > 0 && n_zlevels != nz) { - amrex::Print() << "You supplied " << n_zlevels << " z_levels " << std::endl; - amrex::Print() << "but n_cell in the z-direction is " << nz << std::endl; - amrex::Abort("You must specify a z_level for every value of k"); - } + //std::string terrain_type = "custom"; amrex::Gpu::DeviceVector z_levels_d; z_levels_d.resize(nz); diff --git a/Source/prob_common.H b/Source/prob_common.H index 9f6c874eb..9e501ca88 100644 --- a/Source/prob_common.H +++ b/Source/prob_common.H @@ -107,16 +107,19 @@ public: /** * Function to perform custom initialization of terrain * + * Note: Terrain functionality can also be used to provide grid stretching. + * * @param[in] geom container for geometric information * @param[out] z_phys_nd height coordinate at nodes * @param[in] time current time */ virtual void init_custom_terrain (const amrex::Geometry& /*geom*/, - amrex::MultiFab& /*z_phys_nd*/, + amrex::MultiFab& z_phys_nd, const amrex::Real& /*time*/) { - amrex::Error("Should never call init_custom_terrain for "+name()+" problem"); + amrex::Print() << "Initialized flat terrain at z=0" << std::endl; + z_phys_nd.setVal(0.0); } #ifdef ERF_USE_TERRAIN_VELOCITY