Properly initialize reference values for Rayleigh damping if using Ra…

…yleigh + input_sounding + terrain (#1429) * Remove duplicative print statement * Clarify comment * Properly initialize rayleigh from sounding when there's terrain * Remove out-of-date notification
erf-model · Feb 7, 2024 · e47a0f3 · e47a0f3
1 parent bfd928f
commit e47a0f3
Show file tree

Hide file tree

Showing 5 changed files with 54 additions and 28 deletions.
diff --git a/Source/ERF.cpp b/Source/ERF.cpp
@@ -506,10 +506,6 @@ ERF::InitData ()
         if (solverChoice.use_terrain) {
             if (init_type == "ideal") {
                 amrex::Abort("We do not currently support init_type = ideal with terrain");
-            } else if (init_type == "input_sounding") {
-                amrex::Print() << "Note: use_terrain==true with input_sounding"
-                    << " -- the expected use case is to enable grid stretching"
-                    << std::endl;
             }
         }
 
@@ -665,7 +661,9 @@ ERF::InitData ()
         initRayleigh();
         if (init_type == "input_sounding")
         {
-            // overwrite Ubar, Tbar, and thetabar with input profiles
+            // Overwrite ubar, vbar, and thetabar with input profiles; wbar is
+            // assumed to be 0. Note: the tau coefficient set by
+            // prob->erf_init_rayleigh() is still used
             bool restarting = (!restart_chkfile.empty());
             setRayleighRefFromSounding(restarting);
         }

diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp
@@ -5,6 +5,7 @@
 #include <ERF.H>
 #include <TileNoZ.H>
 #include <prob_common.H>
+#include <Utils/ParFunctions.H>
 
 using namespace amrex;
 
@@ -72,18 +73,29 @@ ERF::setRayleighRefFromSounding (bool restarting)
     for (int lev = 0; lev <= finest_level; lev++)
     {
         const int khi = geom[lev].Domain().bigEnd()[2];
-        const auto *const prob_lo = geom[lev].ProbLo();
-        const auto *const dx = geom[lev].CellSize();
+        Vector<Real> zcc(khi+1);
+
+        if (z_phys_cc[lev]) {
+            // use_terrain=1
+            // calculate the damping strength based on the max height at each k
+            reduce_to_max_per_level(zcc, z_phys_cc[lev]);
+        } else {
+            const auto *const prob_lo = geom[lev].ProbLo();
+            const auto *const dx = geom[lev].CellSize();
+            for (int k = 0; k <= khi; k++)
+            {
+                zcc[k] = prob_lo[2] + (k+0.5) * dx[2];
+            }
+        }
 
         for (int k = 0; k <= khi; k++)
         {
-            const Real z = prob_lo[2] + (k + 0.5) * dx[2];
-            h_rayleigh_ptrs[lev][Rayleigh::ubar][k]          = interpolate_1d(z_inp_sound, U_inp_sound, z, inp_sound_size);
-            h_rayleigh_ptrs[lev][Rayleigh::vbar][k]          = interpolate_1d(z_inp_sound, V_inp_sound, z, inp_sound_size);
-            h_rayleigh_ptrs[lev][Rayleigh::wbar][k]          = Real(0.0);
-            h_rayleigh_ptrs[lev][Rayleigh::thetabar][k] = interpolate_1d(z_inp_sound, theta_inp_sound, z, inp_sound_size);
+            h_rayleigh_ptrs[lev][Rayleigh::ubar][k]         = interpolate_1d(z_inp_sound, U_inp_sound, zcc[k], inp_sound_size);
+            h_rayleigh_ptrs[lev][Rayleigh::vbar][k]         = interpolate_1d(z_inp_sound, V_inp_sound, zcc[k], inp_sound_size);
+            h_rayleigh_ptrs[lev][Rayleigh::wbar][k]         = Real(0.0);
+            h_rayleigh_ptrs[lev][Rayleigh::thetabar][k] = interpolate_1d(z_inp_sound, theta_inp_sound, zcc[k], inp_sound_size);
             if (h_rayleigh_ptrs[lev][Rayleigh::tau][k] > 0) {
-                                                  amrex::Print() << z << ":" << " tau=" << h_rayleigh_ptrs[lev][Rayleigh::tau][k];
+                                                  amrex::Print() << zcc[k] << ":" << " tau=" << h_rayleigh_ptrs[lev][Rayleigh::tau][k];
                 if (solverChoice.rayleigh_damp_U) amrex::Print() << " ubar    = " << h_rayleigh_ptrs[lev][Rayleigh::ubar][k];
                 if (solverChoice.rayleigh_damp_V) amrex::Print() << " vbar    = " << h_rayleigh_ptrs[lev][Rayleigh::vbar][k];
                 if (solverChoice.rayleigh_damp_W) amrex::Print() << " wbar    = " << h_rayleigh_ptrs[lev][Rayleigh::wbar][k];

diff --git a/Source/Prob/init_rayleigh_damping.H b/Source/Prob/init_rayleigh_damping.H
@@ -1,3 +1,5 @@
+#include <Utils/ParFunctions.H>
+
 /**
  * Initialize a Rayleigh damping layer with the same structure as in WRF, based
  * on Durran and Klemp 1983
@@ -13,20 +15,7 @@ erf_init_rayleigh (amrex::Vector<amrex::Vector<amrex::Real> >& rayleigh_ptrs,
     if (z_phys_cc) {
         // use_terrain=1
         // calculate the damping strength based on the max height at each k
-        amrex::MultiArray4<const amrex::Real> const& ma_z_phys = z_phys_cc->const_arrays();
-        for (int k = 0; k <= khi; k++) {
-            zcc[k] = amrex::ParReduce(amrex::TypeList<amrex::ReduceOpMax>{},
-                                      amrex::TypeList<amrex::Real>{},
-                                      *z_phys_cc,
-                                      amrex::IntVect(0),
-            [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int) noexcept
-                -> amrex::GpuTuple<amrex::Real>
-            {
-                return { ma_z_phys[box_no](i,j,k) };
-            });
-        }
-        amrex::ParallelDescriptor::ReduceRealMax(zcc.data(), zcc.size());
-
+        reduce_to_max_per_level(zcc, z_phys_cc);
     } else {
         const amrex::Real* prob_lo = geom.ProbLo();
         const auto dx              = geom.CellSize();
@@ -52,7 +41,6 @@ erf_init_rayleigh (amrex::Vector<amrex::Vector<amrex::Real> >& rayleigh_ptrs,
             rayleigh_ptrs[Rayleigh::vbar][k]     = parms.V_0;
             rayleigh_ptrs[Rayleigh::wbar][k]     = parms.W_0;
             rayleigh_ptrs[Rayleigh::thetabar][k] = parms.T_0;
-            amrex::Print() << "  " << zcc[k] << " " << sinefac*sinefac << std::endl;
         }
         else
         {

diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package
@@ -13,6 +13,8 @@ CEXE_headers += Interpolation.H
 CEXE_headers += Interpolation_1D.H
 CEXE_sources += TerrainMetrics.cpp
 
+CEXE_headers += ParFunctions.H
+
 CEXE_headers += Sat_methods.H
 CEXE_headers += Water_vapor_saturation.H
 CEXE_headers += DirectionSelector.H

diff --git a/Source/Utils/ParFunctions.H b/Source/Utils/ParFunctions.H
@@ -0,0 +1,26 @@
+#ifndef ParFunctions_H
+#define ParFunctions_H
+
+/**
+ * Reduce a multifab to a vector of values at height levels
+ */
+AMREX_FORCE_INLINE
+void reduce_to_max_per_level (amrex::Vector<amrex::Real>& v,
+                              std::unique_ptr<amrex::MultiFab>& mf)
+{
+    amrex::MultiArray4<const amrex::Real> const& ma = mf->const_arrays();
+    for (int k = 0; k < v.size(); k++) {
+        v[k] = amrex::ParReduce(amrex::TypeList<amrex::ReduceOpMax>{},
+                                amrex::TypeList<amrex::Real>{},
+                                *mf,
+                                amrex::IntVect(0),
+        [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int) noexcept
+            -> amrex::GpuTuple<amrex::Real>
+        {
+            return { ma[box_no](i,j,k) };
+        });
+    }
+    amrex::ParallelDescriptor::ReduceRealMax(v.data(), v.size());
+}
+
+#endif /* ParFunctions_H */