Updates for terrain / grid stretching (erf-model#1361)

* Enable use_terrain + input_sounding * Init zlevels_stag before amrex_probinit; warn if levels dont match problem domain It may make sense to update the problem domain if possible * Set the correct ztop from zlevels_stag init_terrain_grid() used to use ztop = geom.ProbHiArray()[2], which comes from the input file. For a stretched grid or user-specified "terrain_z_levels", there could have been a mismatch between ztop and the actual requested top of the domain. * Implement default init_custom_terrain() consistent with flat terrain defintions in most prob.cpps * Remove unnecessary override of init_custom_terrain Note: The regression tests DynamicRefinement, IsentropicVortex, and ScalarAdvDiff all explicitly initialized z_phys_nd to be k*dz and then called FillBoundary(). However, if use_terrain==true (and init_custom_terrain() was called), then the default behavior should be equivalent: init_zlevels() called without grid stretching, then init_terrain_grid() called to populate z_phys_nd and then FillBoundary() * Fix bug in advective flux in x when using higher-order + terrain * Cleanup * Update gold data after bug fix ed6ae8c * Remove zbot==0 assumption; inputs heights are absolute, not AGL * Clarify usage of zlevels_stag a.k.a. "terrain_z_levels" * Properly check that sounding is increasing; use zbot instead of 0; cleanup * Implement terrain-aware input-sounding initialization Notes: - Nominal z levels (e.g., from input terrain_z_levels) that are used for terrain smoothing should range from 0 (on the surface) to ztop - The same levels are used to calculate the base state by integrating the hydrostatic equation through a column of air - In the case of terrain or grid stretching, the base state profiles may have non-uniform vertical grid spacing - These 1-D profiles are then interpolated to the 3-D terrain grid; the state and hse fab arrays are interpolated at z_phys_cc whereas the velocity fab arrays are interpolated at face centers calculated from z_phys_nd --------- Co-authored-by: Ann Almgren <[email protected]>
ewquon · Jan 9, 2024 · d8f1aca · d8f1aca
1 parent 5786ab5
commit d8f1aca
Show file tree

Hide file tree

Showing 61 changed files with 178 additions and 750 deletions.
diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst
@@ -172,7 +172,8 @@ 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.
+be specified with the **erf.terrain_z_levels** parameter, which should vary
+from 0 (at the surface) to the domain height.
 
 .. _list-of-parameters-3:
 

diff --git a/Exec/ABL/prob.H b/Exec/ABL/prob.H
@@ -68,11 +68,6 @@ public:
         amrex::Array4<amrex::Real const> const& mf_v,
         const SolverChoice& sc) override;
 
-    void init_custom_terrain (
-        const amrex::Geometry& geom,
-        amrex::MultiFab& z_phys_nd,
-        const amrex::Real& time) override;
-
 protected:
     std::string name() override { return "ABL"; }
 

diff --git a/Exec/ABL/prob.cpp b/Exec/ABL/prob.cpp
@@ -164,28 +164,3 @@ Problem::init_custom_pert(
     }
   });
 }
-
-void
-Problem::init_custom_terrain(
-    const Geometry& /*geom*/,
-    MultiFab& z_phys_nd,
-    const Real& /*time*/)
-{
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        Array4<Real> z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Flat terrain with z = 0 at k = 0
-            z_arr(i,j,0) = 0.;
-        });
-    }
-}
diff --git a/Exec/ABL_input_sounding/prob.H b/Exec/ABL_input_sounding/prob.H
@@ -67,11 +67,6 @@ public:
         amrex::Array4<amrex::Real const> const& mf_v,
         const SolverChoice& sc) override;
 
-    void init_custom_terrain (
-        const amrex::Geometry& geom,
-        amrex::MultiFab& z_phys_nd,
-        const amrex::Real& time) override;
-
 protected:
     std::string name() override { return "ABL with input sounding"; }
 

diff --git a/Exec/ABL_input_sounding/prob.cpp b/Exec/ABL_input_sounding/prob.cpp
@@ -157,28 +157,3 @@ Problem::init_custom_pert(
     }
   });
 }
-
-void
-Problem::init_custom_terrain(
-    const Geometry& /*geom*/,
-    MultiFab& z_phys_nd,
-    const Real& /*time*/)
-{
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        Array4<Real> z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Flat terrain with z = 0 at k = 0
-            z_arr(i,j,0) = 0.;
-        });
-    }
-}
diff --git a/Exec/DevTests/MiguelDev/prob.H b/Exec/DevTests/MiguelDev/prob.H
@@ -57,11 +57,6 @@ public:
         amrex::Array4<amrex::Real const> const& mf_v,
         const SolverChoice& sc) override;
 
-    void init_custom_terrain (
-        const amrex::Geometry& geom,
-        amrex::MultiFab& z_phys_nd,
-        const amrex::Real& time) override;
-
 protected:
     std::string name() override { return "ABL test"; }
 

diff --git a/Exec/DevTests/MiguelDev/prob.cpp b/Exec/DevTests/MiguelDev/prob.cpp
@@ -140,31 +140,3 @@ Problem::init_custom_pert(
   });
 
 }
-
-void
-Problem::init_custom_terrain(
-    const Geometry& /*geom*/,
-    MultiFab& z_phys_nd,
-    const Real& /*time*/)
-{
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Bottom of domain
-    int k0 = 0;
-
-    for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Flat terrain with z = 0 at k = 0
-            z_arr(i,j,k0) = 0.0;
-        });
-    }
-}
diff --git a/Exec/DevTests/MultiBlock/prob.H b/Exec/DevTests/MultiBlock/prob.H
@@ -46,11 +46,6 @@ public:
         amrex::Array4<amrex::Real const> const& mf_v,
         const SolverChoice& sc) override;
 
-    void init_custom_terrain (
-        const amrex::Geometry& geom,
-        amrex::MultiFab& z_phys_nd,
-        const amrex::Real& time) override;
-
 protected:
     std::string name() override { return "MultiBlock"; }
 

diff --git a/Exec/DevTests/MultiBlock/prob.cpp b/Exec/DevTests/MultiBlock/prob.cpp
@@ -111,31 +111,3 @@ Problem::init_custom_pert(
 
   amrex::Gpu::streamSynchronize();
 }
-
-void
-Problem::init_custom_terrain(
-    const Geometry& /*geom*/,
-    MultiFab& z_phys_nd,
-    const Real& /*time*/)
-{
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Bottom of domain
-    int k0 = 0;
-
-    for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Flat terrain with z = 0 at k = 0
-            z_arr(i,j,k0) = 0.0;
-        });
-    }
-}
diff --git a/Exec/LLJ/prob.H b/Exec/LLJ/prob.H
@@ -39,11 +39,6 @@ public:
         amrex::Array4<amrex::Real const> const& mf_v,
         const SolverChoice& sc) override;
 
-    void init_custom_terrain (
-        const amrex::Geometry& geom,
-        amrex::MultiFab& z_phys_nd,
-        const amrex::Real& time) override;
-
 protected:
     std::string name() override { return "Low-Level Jet"; }
 

diff --git a/Exec/LLJ/prob.cpp b/Exec/LLJ/prob.cpp
@@ -42,28 +42,3 @@ Problem::init_custom_pert(
 {
   amrex::Print() << "Dummy function..Needed for linking" << std::endl;
 }
-
-void
-Problem::init_custom_terrain(
-    const Geometry& /*geom*/,
-    MultiFab& z_phys_nd,
-    const Real& /*time*/)
-{
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        Array4<Real> z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Flat terrain with z = 0 at k = 0
-            z_arr(i,j,0) = 0.;
-        });
-    }
-}
diff --git a/Exec/Radiation/prob.H b/Exec/Radiation/prob.H
@@ -46,11 +46,6 @@ public:
         amrex::Array4<amrex::Real const> const& mf_v,
         const SolverChoice& sc) override;
 
-    void init_custom_terrain (
-        const amrex::Geometry& geom,
-        amrex::MultiFab& z_phys_nd,
-        const amrex::Real& time) override;
-
     void erf_init_rayleigh (
         amrex::Vector<amrex::Real>& tau,
         amrex::Vector<amrex::Real>& ubar,

diff --git a/Exec/Radiation/prob.cpp b/Exec/Radiation/prob.cpp
@@ -212,34 +212,6 @@ Problem::init_custom_pert(
   amrex::Gpu::streamSynchronize();
 }
 
-void
-Problem::init_custom_terrain(
-    const Geometry& /*geom*/,
-    MultiFab& z_phys_nd,
-    const Real& /*time*/)
-{
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Bottom of domain
-    int k0 = 0;
-
-    for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Flat terrain with z = 0 at k = 0
-            z_arr(i,j,k0) = 0.0;
-        });
-    }
-}
-
 void
 Problem::erf_init_rayleigh(
     amrex::Vector<amrex::Real>& tau,

diff --git a/Exec/RegTests/Bubble/prob.H b/Exec/RegTests/Bubble/prob.H
@@ -72,11 +72,6 @@ public:
                                   std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                                   amrex::Geometry const& geom) override;
 
-    void init_custom_terrain (
-        const amrex::Geometry& geom,
-        amrex::MultiFab& z_phys_nd,
-        const amrex::Real& time) override;
-
     AMREX_FORCE_INLINE
     AMREX_GPU_HOST_DEVICE
     amrex::Real compute_theta (const amrex::Real T_b, const amrex::Real p_b);

diff --git a/Exec/RegTests/Bubble/prob.cpp b/Exec/RegTests/Bubble/prob.cpp
@@ -457,31 +457,3 @@ Problem::init_custom_pert(
 
     amrex::Gpu::streamSynchronize();
 }
-
-void
-Problem::init_custom_terrain(
-    const Geometry& /*geom*/,
-    MultiFab& z_phys_nd,
-    const Real& /*time*/)
-{
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Bottom of domain
-    int k0 = 0;
-
-    for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Flat terrain with z = 0 at k = 0
-            z_arr(i,j,k0) = 0.0;
-        });
-    }
-}
diff --git a/Exec/RegTests/CouetteFlow/prob.H b/Exec/RegTests/CouetteFlow/prob.H
@@ -41,11 +41,6 @@ public:
         amrex::Array4<amrex::Real const> const& mf_v,
         const SolverChoice& sc) override;
 
-    void init_custom_terrain (
-        const amrex::Geometry& geom,
-        amrex::MultiFab& z_phys_nd,
-        const amrex::Real& time) override;
-
 protected:
     std::string name() override { return "Couette Flow"; }
 

diff --git a/Exec/RegTests/CouetteFlow/prob.cpp b/Exec/RegTests/CouetteFlow/prob.cpp
@@ -75,28 +75,3 @@ Problem::init_custom_pert(
     z_vel(i, j, k) = parms.w_0;
   });
 }
-
-void
-Problem::init_custom_terrain(
-    const Geometry& /*geom*/,
-    MultiFab& z_phys_nd,
-    const Real& /*time*/)
-{
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        Array4<Real> z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Flat terrain with z = 0 at k = 0
-            z_arr(i,j,0) = 0.;
-        });
-    }
-}
diff --git a/Exec/RegTests/DensityCurrent/prob.H b/Exec/RegTests/DensityCurrent/prob.H
@@ -44,11 +44,6 @@ public:
         amrex::Array4<amrex::Real const> const& mf_v,
         const SolverChoice& sc) override;
 
-    void init_custom_terrain (
-        const amrex::Geometry& geom,
-        amrex::MultiFab& z_phys_nd,
-        const amrex::Real& time) override;
-
 protected:
     std::string name() override { return "Density Current"; }
 

diff --git a/Exec/RegTests/DensityCurrent/prob.cpp b/Exec/RegTests/DensityCurrent/prob.cpp
@@ -148,32 +148,3 @@ Problem::init_custom_pert(
 
   amrex::Gpu::streamSynchronize();
 }
-
-void
-Problem::init_custom_terrain(
-    const Geometry& /*geom*/,
-    MultiFab& z_phys_nd,
-    const Real& /*time*/)
-{
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Bottom of domain
-    int k0 = 0;
-
-    for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Flat terrain with z = 0 at k = 0
-            z_arr(i,j,k0) = 0.0;
-        });
-    }
-}
-