Updates for WOA Problem (erf-model#1807)

* Handle domains with xlo != 0 * Fix indexing gotcha for nodal grid * Add back in factor of h_m; cleanup * Allow gamma_m to be set at runtime for STF terrain * Generalize WoA terrain * Use foextrap for conserved vars if inflow density is not set * Make WoA work with input_sounding and refinement boxes * Allow refinement boxes to be specified based on nominal zlevels_stag when there is terrain * cleanup * Move prob.zdamp and prob.dampceof to SolverChoice; get rid of 1D Rayleigh tau * Update source terms for rayleigh damping * Calculate Rayleigh damping layer using nominal z levels, accounting for refinement levels * Print out more helpful information for setting up refinement boxes * Update inputs * Cleanup * Update golds for updated STF terrain smoothing * Set previously hard-coded params * Handle case of trying to smooth flat terrain * Be extra safe with divide by 0
AMLattanzi · Sep 16, 2024 · db1faa5 · db1faa5
1 parent a9f72b4
commit db1faa5
Show file tree

Hide file tree

Showing 65 changed files with 413 additions and 824 deletions.
diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst
@@ -1195,7 +1195,11 @@ Examples of Usage
     BTF is used when generating the terrain following coordinate.
 
 -  **erf.terrain_smoothing**  = 1
-    STF is used when generating the terrain following coordinate.
+    STF is used when generating the terrain following coordinate. Additionally,
+    ``erf.terrain_gamma_m`` (default=0.5) may be used to set the minimum
+    allowable fractional grid spacing. From Klemp 2011, MWR: "Values of 0.5-0.6
+    seem to work best in 2D applications, while values about half this
+    magnitude appear better for 3D real-terrain simulations."
 
 -  **erf.terrain_smoothing**  = 2
     Sullivan TF is used when generating the terrain following coordinate.

diff --git a/Exec/ABL/ERF_prob.H b/Exec/ABL/ERF_prob.H
@@ -34,10 +34,6 @@ struct ProbParm : ProbParmDefaults {
   amrex::Real pert_periods_V = 5.0;
   amrex::Real pert_ref_height = 100.0;
 
-  // rayleigh damping
-  amrex::Real dampcoef = 0.2; // inverse time scale [1/s]
-  amrex::Real zdamp = 500.0; // damping depth [m] from model top
-
   // helper vars
   amrex::Real aval;
   amrex::Real bval;

diff --git a/Exec/ABL/ERF_prob.cpp b/Exec/ABL/ERF_prob.cpp
@@ -39,9 +39,6 @@ Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi)
   parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height;
   parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height;
 
-  pp.query("dampcoef", parms.dampcoef);
-  pp.query("zdamp", parms.zdamp);
-
   init_base_parms(parms.rho_0, parms.T_0);
 }
 

diff --git a/Exec/ABL/inputs_GABLS1_deardorff b/Exec/ABL/inputs_GABLS1_deardorff
@@ -82,9 +82,8 @@ erf.rayleigh_damp_U = true
 erf.rayleigh_damp_V = true
 erf.rayleigh_damp_W = true
 erf.rayleigh_damp_T = true
-
-prob.dampcoef = 0.2 # [1/s] following FastEddy
-prob.zdamp = 100. # from Beare et al. 2006, "most models applied gravity wave damping above 300 m"
+erf.rayleigh_dampcoef = 0.2 # [1/s] following FastEddy
+erf.rayleigh_zdamp = 100. # from Beare et al. 2006, "most models applied gravity wave damping above 300 m"
 
 # Initial conditions from Beare et al. 2006
 prob.KE_0            = 0.4 # [m2/s2]

diff --git a/Exec/ABL_input_sounding/ERF_prob.H b/Exec/ABL_input_sounding/ERF_prob.H
@@ -29,10 +29,6 @@ struct ProbParm : ProbParmDefaults {
   amrex::Real pert_periods_V = 5.0;
   amrex::Real pert_ref_height = 100.0;
 
-  // rayleigh damping
-  amrex::Real dampcoef = 0.003; // inverse time scale [1/s]
-  amrex::Real zdamp = 5000.0; // damping depth [m] from model top
-
   // helper vars
   amrex::Real aval;
   amrex::Real bval;

diff --git a/Exec/ABL_input_sounding/ERF_prob.cpp b/Exec/ABL_input_sounding/ERF_prob.cpp
@@ -35,9 +35,6 @@ Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi)
   parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height;
   parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height;
 
-  pp.query("dampcoef", parms.dampcoef);
-  pp.query("zdamp", parms.zdamp);
-
   init_base_parms(parms.rho_0, parms.T_0);
 }
 

diff --git a/Exec/DevTests/ABL_Rayleigh_check/inputs_rayleigh_check b/Exec/DevTests/ABL_Rayleigh_check/inputs_rayleigh_check
@@ -54,9 +54,9 @@ erf.molec_diff_type = "None"
 
 # RAYLEIGH DAMPING
 erf.rayleigh_damp_U = true
+erf.rayleigh_dampcoef = 0.02 # timescale is 50 s
+erf.rayleigh_zdamp = 9e99 # no vertical variation in damping layer strength
 
 # PROBLEM CHOICE
-prob.zdamp = 9e99 # no vertical variation in damping layer strength
-prob.dampcoef = 0.02 # timescale is 50 s
 prob.U_0 = -10.0 # subtract out mean wind from profile to start with U = 0
 
diff --git a/Exec/DevTests/ABL_with_WW3/ERF_prob.H b/Exec/DevTests/ABL_with_WW3/ERF_prob.H
@@ -35,10 +35,6 @@ struct ProbParm : ProbParmDefaults {
   amrex::Real pert_periods_V = 5.0;
   amrex::Real pert_ref_height = 100.0;
 
-  // rayleigh damping
-  amrex::Real dampcoef = 0.2; // inverse time scale [1/s]
-  amrex::Real zdamp = 500.0; // damping depth [m] from model top
-
   // helper vars
   amrex::Real aval;
   amrex::Real bval;

diff --git a/Exec/DevTests/ABL_with_WW3/ERF_prob.cpp b/Exec/DevTests/ABL_with_WW3/ERF_prob.cpp
@@ -40,9 +40,6 @@ Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi)
   parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height;
   parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height;
 
-  pp.query("dampcoef", parms.dampcoef);
-  pp.query("zdamp", parms.zdamp);
-
   init_base_parms(parms.rho_0, parms.T_0);
 }
 

diff --git a/Exec/DevTests/EB_Test/ERF_prob.H b/Exec/DevTests/EB_Test/ERF_prob.H
@@ -60,8 +60,9 @@ public:
 
     void erf_init_rayleigh (
         amrex::Vector<amrex::Vector<amrex::Real> >& rayleigh_ptrs,
-        amrex::Geometry      const& geom,
-        std::unique_ptr<amrex::MultiFab>& z_phys_cc) override;
+        amrex::Geometry const& geom,
+        std::unique_ptr<amrex::MultiFab>& z_phys_nd,
+        amrex::Real zdamp) override;
 
 protected:
     std::string name() override { return "Scalar Advection/Diffusion"; }

diff --git a/Exec/DevTests/EB_Test/ERF_prob.cpp b/Exec/DevTests/EB_Test/ERF_prob.cpp
@@ -38,15 +38,15 @@ Problem::Problem()
 void
 Problem::erf_init_rayleigh(
     amrex::Vector<amrex::Vector<amrex::Real> >& rayleigh_ptrs,
-    amrex::Geometry      const& geom,
-    std::unique_ptr<MultiFab>& /*z_phys_cc*/)
+    amrex::Geometry const& geom,
+    std::unique_ptr<MultiFab>& /*z_phys_nd*/,
+    amrex::Real /*zdamp*/)
 {
   const int khi = geom.Domain().bigEnd()[2];
 
   // We just use these values to test the Rayleigh damping
   for (int k = 0; k <= khi; k++)
   {
-      rayleigh_ptrs[Rayleigh::tau][k]      = 1.0;
       rayleigh_ptrs[Rayleigh::ubar][k]     = 2.0;
       rayleigh_ptrs[Rayleigh::vbar][k]     = 1.0;
       rayleigh_ptrs[Rayleigh::wbar][k]     = 0.0;

diff --git a/Exec/DevTests/MovingTerrain/ERF_prob.H b/Exec/DevTests/MovingTerrain/ERF_prob.H
@@ -46,8 +46,9 @@ public:
 
     void erf_init_rayleigh (
         amrex::Vector<amrex::Vector<amrex::Real> >& rayleigh_ptrs,
-        amrex::Geometry      const& geom,
-        std::unique_ptr<amrex::MultiFab>& z_phys_cc) override;
+        amrex::Geometry const& geom,
+        std::unique_ptr<amrex::MultiFab>& z_phys_nd,
+        amrex::Real zdamp) override;
 
 protected:
     std::string name() override { return "MovingTerrain"; }

diff --git a/Exec/DevTests/MovingTerrain/ERF_prob.cpp b/Exec/DevTests/MovingTerrain/ERF_prob.cpp
@@ -123,12 +123,12 @@ Problem::init_custom_pert(
 void
 Problem::erf_init_rayleigh (Vector<amrex::Vector<amrex::Real> >& rayleigh_ptrs,
                             Geometry const& geom,
-                            std::unique_ptr<MultiFab>& /*z_phys_cc*/)
+                            std::unique_ptr<MultiFab>& /*z_phys_nd*/,
+                            amrex::Real /*zdamp*/)
 {
     const int khi = geom.Domain().bigEnd()[2];
     for (int k = 0; k <= khi; k++)
     {
-        rayleigh_ptrs[Rayleigh::tau][k]      =   0.0;
         rayleigh_ptrs[Rayleigh::ubar][k]     =   0.0;
         rayleigh_ptrs[Rayleigh::vbar][k]     =   0.0;
         rayleigh_ptrs[Rayleigh::wbar][k]     =   0.0;
@@ -138,7 +138,6 @@ Problem::erf_init_rayleigh (Vector<amrex::Vector<amrex::Real> >& rayleigh_ptrs,
     // Damping above k = 60
     for (int k = 60; k <= khi; k++)
     {
-        rayleigh_ptrs[Rayleigh::tau][k]      =  10.0;  // Remember that this gets multiplied by dt
         rayleigh_ptrs[Rayleigh::ubar][k]     =   2.0;
         rayleigh_ptrs[Rayleigh::vbar][k]     =   1.0;
         rayleigh_ptrs[Rayleigh::wbar][k]     =   0.0;

diff --git a/Exec/DevTests/TemperatureSource/ERF_prob.H b/Exec/DevTests/TemperatureSource/ERF_prob.H
@@ -30,10 +30,6 @@ struct ProbParm : ProbParmDefaults {
   amrex::Real pert_periods_V = 5.0;
   amrex::Real pert_ref_height = 100.0;
 
-  // rayleigh damping
-  amrex::Real dampcoef = 0.2; // inverse time scale [1/s]
-  amrex::Real zdamp = 500.0; // damping depth [m] from model top
-
   //==============================================
   // USER-DEFINED INPUTS
   // source terms

diff --git a/Exec/DevTests/TemperatureSource/ERF_prob.cpp b/Exec/DevTests/TemperatureSource/ERF_prob.cpp
@@ -37,9 +37,6 @@ Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi)
   parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height;
   parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height;
 
-  pp.query("dampcoef", parms.dampcoef);
-  pp.query("zdamp", parms.zdamp);
-
   //===========================================================================
   // READ USER-DEFINED INPUTS
   pp.get("advection_heating_rate", parms.advection_heating_rate);

diff --git a/Exec/DevTests/TropicalCyclone/ERF_prob.H b/Exec/DevTests/TropicalCyclone/ERF_prob.H
@@ -17,12 +17,9 @@ struct ProbParm : ProbParmDefaults {
     amrex::Real RZERO = 800.0;  // Radius of zero wind speed
     amrex::Real ZZERO = 2000.0; // Height of zero wind speed
 
-    // rayleigh damping
     amrex::Real U_0 = 0.0;
     amrex::Real V_0 = 0.0;
     amrex::Real W_0 = 0.0;
-    amrex::Real dampcoef = 0.2; // inverse time scale [1/s]
-    amrex::Real zdamp = 500.0; // damping depth [m] from model top
 
 }; // namespace ProbParm
 

diff --git a/Exec/DevTests/TropicalCyclone/ERF_prob.cpp b/Exec/DevTests/TropicalCyclone/ERF_prob.cpp
@@ -28,9 +28,6 @@ Problem::Problem()
     pp.query("RZERO", parms.RZERO);
     pp.query("ZZERO", parms.ZZERO);
 
-    pp.query("dampcoef", parms.dampcoef);
-    pp.query("zdamp", parms.zdamp);
-
     init_base_parms(parms.rho_0, parms.T_0);
 }
 

diff --git a/Exec/DevTests/TropicalCyclone/inputs b/Exec/DevTests/TropicalCyclone/inputs
@@ -83,8 +83,8 @@ erf.rotational_time_period = 86164.0900027328
 
 # RAYLEIGH DAMPING
 erf.rayleigh_damp_W = true
-prob.dampcoef = 0.2
-prob.zdamp = 5000.
+erf.rayleigh_dampcoef = 0.2
+erf.rayleigh_zdamp = 5000.
 
 # PROBLEM PARAMETERS
 # initial vortex

diff --git a/Exec/EWP/ERF_prob.H b/Exec/EWP/ERF_prob.H
@@ -30,10 +30,6 @@ struct ProbParm : ProbParmDefaults {
   amrex::Real pert_periods_V = 5.0;
   amrex::Real pert_ref_height = 100.0;
 
-  // rayleigh damping
-  amrex::Real dampcoef = 0.2; // inverse time scale [1/s]
-  amrex::Real zdamp = 500.0; // damping depth [m] from model top
-
   // helper vars
   amrex::Real aval;
   amrex::Real bval;

diff --git a/Exec/EWP/ERF_prob.cpp b/Exec/EWP/ERF_prob.cpp
@@ -36,9 +36,6 @@ Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi)
   parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height;
   parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height;
 
-  pp.query("dampcoef", parms.dampcoef);
-  pp.query("zdamp", parms.zdamp);
-
   init_base_parms(parms.rho_0, parms.T_0);
 }
 

diff --git a/Exec/Radiation/ERF_prob.H b/Exec/Radiation/ERF_prob.H
@@ -49,8 +49,9 @@ public:
 
     void erf_init_rayleigh (
         amrex::Vector<amrex::Vector<amrex::Real> >& rayleigh_ptrs,
-        amrex::Geometry      const& geom,
-        std::unique_ptr<amrex::MultiFab>& z_phys_cc) override;
+        amrex::Geometry const& geom,
+        std::unique_ptr<amrex::MultiFab>& z_phys_nd,
+        amrex::Real zdamp) override;
 
 protected:
     std::string name() override { return "Supercell"; }

diff --git a/Exec/Radiation/ERF_prob.cpp b/Exec/Radiation/ERF_prob.cpp
@@ -214,15 +214,15 @@ Problem::init_custom_pert(
 void
 Problem::erf_init_rayleigh(
     amrex::Vector<amrex::Vector<Real> >& rayleigh_ptrs,
-    amrex::Geometry      const& geom,
-    std::unique_ptr<MultiFab>& /*z_phys_cc*/)
+    amrex::Geometry const& geom,
+    std::unique_ptr<MultiFab>& /*z_phys_nd*/,
+    amrex::Real /*zdamp*/)
 {
     const int khi = geom.Domain().bigEnd()[2];
 
     // We just use these values to test the Rayleigh damping
     for (int k = 0; k <= khi; k++)
     {
-        rayleigh_ptrs[Rayleigh::tau][k]      = 1.0;
         rayleigh_ptrs[Rayleigh::ubar][k]     = 2.0;
         rayleigh_ptrs[Rayleigh::vbar][k]     = 1.0;
         rayleigh_ptrs[Rayleigh::wbar][k]     = 0.0;

diff --git a/Exec/RegTests/Bomex/ERF_prob.H b/Exec/RegTests/Bomex/ERF_prob.H
@@ -31,10 +31,6 @@ struct ProbParm : ProbParmDefaults {
     amrex::Real pert_periods_V = 5.0;
     amrex::Real pert_ref_height = 100.0;
 
-    // rayleigh damping
-    amrex::Real dampcoef = 0.2; // inverse time scale [1/s]
-    amrex::Real zdamp = 500.0; // damping depth [m] from model top
-
     //==============================================
     // USER-DEFINED INPUTS
     // source terms

diff --git a/Exec/RegTests/Bomex/ERF_prob.cpp b/Exec/RegTests/Bomex/ERF_prob.cpp
@@ -38,9 +38,6 @@ Problem::Problem (const Real* problo, const Real* probhi)
     parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height;
     parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height;
 
-    pp.query("dampcoef", parms.dampcoef);
-    pp.query("zdamp", parms.zdamp);
-
     //===========================================================================
     // READ USER-DEFINED INPUTS
     pp.query("advection_heating_rate", parms.advection_heating_rate);

diff --git a/Exec/RegTests/Bubble/ERF_prob.H b/Exec/RegTests/Bubble/ERF_prob.H
@@ -31,10 +31,6 @@ struct ProbParm : ProbParmDefaults {
     bool T_pert_is_airtemp = true; // T_pert input is air temperature
     bool perturb_rho = true; // not rho*theta (i.e., p is constant); otherwise perturb rho*theta
 
-    // rayleigh damping
-    amrex::Real dampcoef = 0.0;    // inverse time scale [1/s]
-    amrex::Real zdamp    = 5000.0; // damping depth [m] from model top
-
     // Moist bubble params
     bool do_moist_bubble = false;
     amrex::Real theta_pert  = 2.0;

diff --git a/Exec/RegTests/Bubble/ERF_prob.cpp b/Exec/RegTests/Bubble/ERF_prob.cpp
@@ -29,8 +29,6 @@ Problem::Problem()
   pp.query("T_pert", parms.T_pert);
   pp.query("T_pert_is_airtemp", parms.T_pert_is_airtemp);
   pp.query("perturb_rho", parms.perturb_rho);
-  pp.query("dampcoef", parms.dampcoef);
-  pp.query("zdamp", parms.zdamp);
 
   pp.query("do_moist_bubble", parms.do_moist_bubble);
   pp.query("theta_pert", parms.theta_pert);

diff --git a/Exec/RegTests/ScalarAdvDiff/ERF_prob.H b/Exec/RegTests/ScalarAdvDiff/ERF_prob.H
@@ -56,8 +56,9 @@ public:
 
     void erf_init_rayleigh (
         amrex::Vector<amrex::Vector<amrex::Real> >& rayleigh_ptrs,
-        amrex::Geometry      const& geom,
-        std::unique_ptr<amrex::MultiFab>& z_phys_cc) override;
+        amrex::Geometry const& geom,
+        std::unique_ptr<amrex::MultiFab>& z_phys_nd,
+        amrex::Real zdamp) override;
 
 protected:
     std::string name() override { return "Scalar Advection/Diffusion"; }

diff --git a/Exec/RegTests/ScalarAdvDiff/ERF_prob.cpp b/Exec/RegTests/ScalarAdvDiff/ERF_prob.cpp
@@ -38,15 +38,15 @@ Problem::Problem()
 void
 Problem::erf_init_rayleigh(
     amrex::Vector<amrex::Vector<amrex::Real> >& rayleigh_ptrs,
-    amrex::Geometry      const& geom,
-    std::unique_ptr<MultiFab>& /*z_phys_cc*/)
+    amrex::Geometry const& geom,
+    std::unique_ptr<MultiFab>& /*z_phys_nd*/,
+    amrex::Real /*zdamp*/)
 {
   const int khi = geom.Domain().bigEnd()[2];
 
   // We just use these values to test the Rayleigh damping
   for (int k = 0; k <= khi; k++)
   {
-      rayleigh_ptrs[Rayleigh::tau][k]      = 1.0;
       rayleigh_ptrs[Rayleigh::ubar][k]     = 2.0;
       rayleigh_ptrs[Rayleigh::vbar][k]     = 1.0;
       rayleigh_ptrs[Rayleigh::wbar][k]     = 0.0;

diff --git a/Exec/RegTests/TurbulentInflow/ERF_prob.H b/Exec/RegTests/TurbulentInflow/ERF_prob.H
@@ -35,10 +35,6 @@ struct ProbParm : ProbParmDefaults {
   amrex::Real pert_periods_V = 5.0;
   amrex::Real pert_ref_height = 100.0;
 
-  // rayleigh damping
-  amrex::Real dampcoef = 0.2; // inverse time scale [1/s]
-  amrex::Real zdamp = 500.0; // damping depth [m] from model top
-
   // helper vars
   amrex::Real aval;
   amrex::Real bval;

diff --git a/Exec/RegTests/TurbulentInflow/ERF_prob.cpp b/Exec/RegTests/TurbulentInflow/ERF_prob.cpp
@@ -40,9 +40,6 @@ Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi)
   parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height;
   parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height;
 
-  pp.query("dampcoef", parms.dampcoef);
-  pp.query("zdamp", parms.zdamp);
-
   init_base_parms(parms.rho_0, parms.T_0);
 }
 

diff --git a/Exec/RegTests/WitchOfAgnesi/ERF_prob.H b/Exec/RegTests/WitchOfAgnesi/ERF_prob.H
@@ -8,11 +8,16 @@
 #include "ERF_prob_common.H"
 
 struct ProbParm : ProbParmDefaults {
+  // perturbations to initial conditions
+  amrex::Real rho_0 = 0.0;
+  amrex::Real T_0 = 0.0;
   amrex::Real U_0 = 0.0;
   amrex::Real V_0 = 0.0;
-  amrex::Real W_0 = 0.0;
-  amrex::Real dampcoef = 0.003;
-  amrex::Real zdamp = 1.0;
+  amrex::Real W_0 = 0.0; // needed for rayleigh damping
+
+  // hill parameters
+  amrex::Real hmax = 0.0;   // full hill height, can be negative to simulate valleys
+  amrex::Real L = 0.0;      // hill length at half-height
 }; // namespace ProbParm
 
 class Problem : public ProblemBase