Fix MOST logic and correct issues with MOST average with terrain. (#1468

)

Docs/sphinx_doc/MOST.rst

@@ -136,6 +136,10 @@ In ERF, when the MOST boundary condition is applied, velocity and temperature in
 
       (\rho \theta)_{i,j,-n} = \rho_{i,j,-n} \left[ \frac{(\rho\theta)_{i,j,0}}{\rho_{i,j,0}} - \left. \frac{\tau_{\theta z}}{\rho} \right|_{i,j,0} \frac{\rho_{i,j,0}}{K_{\theta,v,(i,j,0)}} n \Delta z \right]
 
+   Finally, it must be noted that complex terrain will modify the surface normal and tangent vectors. Consequently, the MOST implentation with terrain will require local vector rotations. While the ERF dycore accounts for
+   terrain metrics when computing fluxes (e.g. for advection, diffusion, etc.), the impact of terrain metrics on MOST is still a work in progress. Therefore, running with terrain (``erf.use_terrain = true``) and with MOST
+   (``zlo.type = "Most"``) should be cautioned.
+
 MOST Inputs
 ~~~~~~~~~~~~~~~~~~~
 To evaluate the fluxes with MOST, the surface rougness parameter :math:`z_{0}` must be specified. This quantity may be considered a constant or may be parameterized through the friction velocity :math:`u_{\star}`. ERF supports three methods for parameterizing the surface roughness: ``constant``, ``charnock``, and ``modified_charnock``. The latter two methods parameterize :math:`z_{0} = f(u_{\star})` and are described in `Jimenez & Dudhia, American Meteorological Society, 2018 <https://doi.org/10.1175/JAMC-D-17-0137.1>`_. The rougness calculation method may be specified with

Exec/RegTests/Terrain3d_Hemisphere/inputs_most_test

@@ -18,7 +18,7 @@ zhi.type = "SlipWall"
 # MOST BOUNDARY
 #=================================================================
 zlo.type                   = "Most"
-erf.most.average_policy    = 1       # POLICY FOR AVERAGING
+erf.most.average_policy    = 0       # POLICY FOR AVERAGING
 erf.most.use_normal_vector = false   # USE NORMAL VECTOR W/ TERRAIN
 erf.most.use_interpolation = true    # USE INTERPOLATION ON DESTINATION   
 #erf.most.time_average      = true    # USE TIME AVERAGING
@@ -81,3 +81,7 @@ erf.alpha_T = 0.0 # [m^2/s]
 prob.T_0   = 300.0
 prob.U_0   = 1.0
 prob.rho_0 = 1.16
+prob.U_0_Pert_Mag = 0.01
+prob.V_0_Pert_Mag = 0.01
+prob.W_0_Pert_Mag = 0.0
+prob.pert_ref_height = 1.0

Exec/RegTests/Terrain3d_Hemisphere/prob.H

@@ -12,12 +12,30 @@ struct ProbParm : ProbParmDefaults {
   amrex::Real T_0 = 300.0; // surface temperature == mean potential temperature
   amrex::Real U_0 = 0.0;
   amrex::Real V_0 = 0.0;
+
+  // random initial perturbations (legacy code)
+  amrex::Real U_0_Pert_Mag = 0.0;
+  amrex::Real V_0_Pert_Mag = 0.0;
+  amrex::Real W_0_Pert_Mag = 0.0;
+
+  // divergence-free initial perturbations
+  amrex::Real pert_deltaU = 0.0;
+  amrex::Real pert_deltaV = 0.0;
+  amrex::Real pert_periods_U  = 5.0;
+  amrex::Real pert_periods_V  = 5.0;
+  amrex::Real pert_ref_height = 1.0;
+
+  // helper vars
+  amrex::Real aval;
+  amrex::Real bval;
+  amrex::Real ufac;
+  amrex::Real vfac;
 }; // namespace ProbParm
 
 class Problem : public ProblemBase
 {
 public:
-    Problem();
+    Problem (const amrex::Real* problo, const amrex::Real* probhi);
 
 #include "Prob/init_density_hse_dry_terrain.H"
 

Exec/RegTests/Terrain3d_Hemisphere/prob.cpp

@@ -4,25 +4,33 @@
 using namespace amrex;
 
 std::unique_ptr<ProblemBase>
-amrex_probinit(
-    const amrex_real* /*problo*/,
-    const amrex_real* /*probhi*/)
+amrex_probinit (const amrex_real* problo,
+                const amrex_real* probhi)
 {
-    return std::make_unique<Problem>();
+    return std::make_unique<Problem>(problo, probhi);
 }
 
-Problem::Problem()
+Problem::Problem (const amrex::Real* problo,
+                  const amrex::Real* probhi)
 {
   // Parse params
   amrex::ParmParse pp("prob");
   pp.query("rho_0", parms.rho_0);
   pp.query("U_0", parms.U_0);
+  pp.query("U_0_Pert_Mag", parms.U_0_Pert_Mag);
+  pp.query("V_0_Pert_Mag", parms.V_0_Pert_Mag);
+  pp.query("W_0_Pert_Mag", parms.W_0_Pert_Mag);
+  pp.query("pert_ref_height", parms.pert_ref_height);
+  parms.aval = parms.pert_periods_U * 2.0 * PI / (probhi[1] - problo[1]);
+  parms.bval = parms.pert_periods_V * 2.0 * PI / (probhi[0] - problo[0]);
+  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;
 
   init_base_parms(parms.rho_0, parms.T_0);
 }
 
 void
-Problem::init_custom_pert(
+Problem::init_custom_pert (
     const Box& bx,
     const Box& xbx,
     const Box& ybx,
@@ -34,66 +42,107 @@ Problem::init_custom_pert(
     Array4<Real      > const&,
     Array4<Real      > const&,
     Array4<Real const> const& z_nd,
-    Array4<Real const> const& z_cc,
+    Array4<Real const> const& /*z_cc*/,
     GeometryData const& geomdata,
     Array4<Real const> const& /*mf_m*/,
     Array4<Real const> const& /*mf_u*/,
     Array4<Real const> const& /*mf_v*/,
     const SolverChoice& sc)
 {
-  const int khi = geomdata.Domain().bigEnd()[2];
+    const int khi = geomdata.Domain().bigEnd()[2];
 
     const bool use_moisture = (sc.moisture_type != MoistureType::None);
+    const bool use_terrain  = sc.use_terrain;
 
-  AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);
+    AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);
 
-  // Geometry (note we must include these here to get the data on device)
-  amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-  {
-    const auto prob_lo  = geomdata.ProbLo();
-    const auto dx       = geomdata.CellSize();
-    const Real x = prob_lo[0] + (i + 0.5) * dx[0];
-    const Real z = z_cc(i,j,k);
-
-    // Set scalar = 0 everywhere
-    state(i, j, k, RhoScalar_comp) = 0.0;
+    // Geometry (note we must include these here to get the data on device)
+    ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    {
+        // Set scalar = 0 everywhere
+        state(i, j, k, RhoScalar_comp) = 0.0;
 
-    if (use_moisture) {
-        state(i, j, k, RhoQ1_comp) = 0.0;
-        state(i, j, k, RhoQ2_comp) = 0.0;
-    }
-  });
-
-  // Set the x-velocity
-  amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-  {
-      x_vel(i, j, k) = 10.0;
-  });
-
-  // Set the y-velocity
-  amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-  {
-      y_vel(i, j, k) = 0.0;
-  });
-
-  const auto dx = geomdata.CellSize();
-  amrex::GpuArray<Real, AMREX_SPACEDIM> dxInv;
-  dxInv[0] = 1. / dx[0];
-  dxInv[1] = 1. / dx[1];
-  dxInv[2] = 1. / dx[2];
-
-  // Set the z-velocity from impenetrable condition
-  amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-  {
-      z_vel(i, j, k) = 0.0;
-  });
-
-  amrex::Gpu::streamSynchronize();
+        if (use_moisture) {
+            state(i, j, k, RhoQ1_comp) = 0.0;
+            state(i, j, k, RhoQ2_comp) = 0.0;
+        }
+    });
 
+    // Set the x-velocity
+    ParallelForRNG(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept
+    {
+        const Real* prob_lo = geomdata.ProbLo();
+        const Real* dx = geomdata.CellSize();
+        const Real z = use_terrain ? 0.25*( z_nd(i,j  ,k) + z_nd(i,j  ,k+1)
+                                          + z_nd(i,j+1,k) + z_nd(i,j+1,k+1) )
+            : prob_lo[2] + (k + 0.5) * dx[2];
+
+        x_vel(i, j, k) = 10.0;
+        if ((z <= parms.pert_ref_height) && (parms.U_0_Pert_Mag != 0.0))
+        {
+            Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
+            Real x_vel_prime = (rand_double*2.0 - 1.0)*parms.U_0_Pert_Mag;
+            x_vel(i, j, k) += x_vel_prime;
+        }
+    });
+
+    // Set the y-velocity
+    ParallelForRNG(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept
+    {
+        const Real* prob_lo = geomdata.ProbLo();
+        const Real* dx = geomdata.CellSize();
+        const Real x = prob_lo[0] + (i + 0.5) * dx[0];
+        const Real z = use_terrain ? 0.25*( z_nd(i  ,j,k) + z_nd(i  ,j,k+1)
+                                          + z_nd(i+1,j,k) + z_nd(i+1,j,k+1) )
+            : prob_lo[2] + (k + 0.5) * dx[2];
+
+        // Set the y-velocity
+        y_vel(i, j, k) = 0.0;
+        if ((z <= parms.pert_ref_height) && (parms.V_0_Pert_Mag != 0.0))
+        {
+            Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
+            Real y_vel_prime = (rand_double*2.0 - 1.0)*parms.V_0_Pert_Mag;
+            y_vel(i, j, k) += y_vel_prime;
+        }
+        if (parms.pert_deltaV != 0.0)
+        {
+            const amrex::Real xl = x - prob_lo[0];
+            const amrex::Real zl = z / parms.pert_ref_height;
+            const amrex::Real damp = std::exp(-0.5 * zl * zl);
+            y_vel(i, j, k) += parms.vfac * damp * z * std::cos(parms.bval * xl);
+        }
+    });
+
+    const auto dx = geomdata.CellSize();
+    amrex::GpuArray<Real, AMREX_SPACEDIM> dxInv;
+    dxInv[0] = 1. / dx[0];
+    dxInv[1] = 1. / dx[1];
+    dxInv[2] = 1. / dx[2];
+
+    // Set the z-velocity from impenetrable condition
+    ParallelForRNG(zbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept
+    {
+        const int dom_lo_z = geomdata.Domain().smallEnd()[2];
+        const int dom_hi_z = geomdata.Domain().bigEnd()[2];
+
+        // Set the z-velocity
+        if (k == dom_lo_z || k == dom_hi_z+1)
+        {
+            z_vel(i, j, k) = 0.0;
+        }
+        else if (parms.W_0_Pert_Mag != 0.0)
+        {
+            Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
+            Real z_vel_prime = (rand_double*2.0 - 1.0)*parms.W_0_Pert_Mag;
+            z_vel(i, j, k) = z_vel_prime;
+        }
+    });
+
+    amrex::Gpu::streamSynchronize();
 }
 
 void
-Problem::init_custom_terrain(
+Problem::init_custom_terrain (
     const Geometry& geom,
     MultiFab& z_phys_nd,
     const Real& /*time*/)
@@ -111,7 +160,6 @@ Problem::init_custom_terrain(
 
     // User function parameters
     Real a    = 0.5;
-    Real num  = 8 * a * a * a;
     Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
     Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
 

Source/BoundaryConditions/ABLMost.H

@@ -286,9 +286,9 @@ public:
         MODIFIED_CHARNOCK
     };
 
-    FluxCalcType  flux_type;
-    ThetaCalcType theta_type;
-    RoughCalcType rough_type;
+    FluxCalcType   flux_type{FluxCalcType::MOENG};
+    ThetaCalcType theta_type{ThetaCalcType::ADIABATIC};
+    RoughCalcType rough_type{RoughCalcType::CONSTANT};
 
 private:
     bool use_moisture;

Source/BoundaryConditions/ABLMost.cpp

@@ -174,7 +174,8 @@ ABLMost::compute_most_bcs (const int& lev,
         const auto *const u_mean     = m_ma.get_average(lev,0);
         const auto *const v_mean     = m_ma.get_average(lev,1);
         const auto *const t_mean     = m_ma.get_average(lev,2);
-        const auto *const u_mag_mean = m_ma.get_average(lev,3);
+        const auto *const q_mean     = m_ma.get_average(lev,3);
+        const auto *const u_mag_mean = m_ma.get_average(lev,4);
 
         const auto um_arr  = u_mean->array(mfi);
         const auto vm_arr  = v_mean->array(mfi);
@@ -202,7 +203,6 @@ ABLMost::compute_most_bcs (const int& lev,
                 Box b2d = bx;
                 b2d.setBig(2,zlo-1);
                 int n = RhoTheta_comp;
-
                 ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                 {
                     Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain;

Source/BoundaryConditions/MOSTAverage.cpp

@@ -317,13 +317,15 @@ MOSTAverage::set_k_indices_T()
         for (int lev(0); lev < m_maxlev; lev++) {
             int kmax = m_geom[lev].Domain().bigEnd(2);
             for (MFIter mfi(*m_k_indx[lev], TileNoZ()); mfi.isValid(); ++mfi) {
-                Box npbx  = mfi.tilebox(); npbx.convert({1,1,0});
+                Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
                 const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
                 auto k_arr = m_k_indx[lev]->array(mfi);
                 ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     k_arr(i,j,k) = 0;
-                    Real z_target = d_zref + z_phys_arr(i,j,k);
+                    Real z_bot_face  = 0.25 * ( z_phys_arr(i  ,j  ,k) + z_phys_arr(i+1,j  ,k)
+                                              + z_phys_arr(i  ,j+1,k) + z_phys_arr(i+1,j+1,k) );
+                    Real z_target    = z_bot_face + d_zref;
                     for (int lk(0); lk<=kmax; ++lk) {
                         Real z_lo = 0.25 * ( z_phys_arr(i,j  ,lk  ) + z_phys_arr(i+1,j  ,lk  )
                                            + z_phys_arr(i,j+1,lk  ) + z_phys_arr(i+1,j+1,lk  ) );
@@ -392,7 +394,9 @@ MOSTAverage::set_norm_indices_T()
                 j_arr(i,j,k) = j_new;
 
                 // Search for k (grid is stretched in z)
-                Real z_target = delta_z + z_phys_arr(i,j,k);
+                Real z_bot_face  = 0.25 * ( z_phys_arr(i  ,j  ,k) + z_phys_arr(i+1,j  ,k)
+                                          + z_phys_arr(i  ,j+1,k) + z_phys_arr(i+1,j+1,k) );
+                Real z_target    = z_bot_face + delta_z;
                 for (int lk(0); lk<=kmax; ++lk) {
                     Real z_lo = 0.25 * ( z_phys_arr(i_new,j_new  ,lk  ) + z_phys_arr(i_new+1,j_new  ,lk  )
                                        + z_phys_arr(i_new,j_new+1,lk  ) + z_phys_arr(i_new+1,j_new+1,lk  ) );
@@ -479,7 +483,7 @@ MOSTAverage::set_norm_positions_T()
         const auto dxInv  = m_geom[lev].InvCellSizeArray();
         IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
         for (MFIter mfi(*m_x_pos[lev], TileNoZ()); mfi.isValid(); ++mfi) {
-            Box npbx  = mfi.tilebox(); npbx.convert({1,1,0});
+            Box npbx  = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
             Box gpbx  = mfi.growntilebox(ng);
             RealBox grb{gpbx,dx.data(),base.dataPtr()};
 
@@ -506,7 +510,9 @@ MOSTAverage::set_norm_positions_T()
                 // Final position at end of vector
                 x_pos_arr(i,j,k) = x0 + delta_x;
                 y_pos_arr(i,j,k) = y0 + delta_y;
-                z_pos_arr(i,j,k) = z_phys_arr(i,j,k) + delta_z;
+                Real z_bot_face  = 0.25 * ( z_phys_arr(i  ,j  ,k) + z_phys_arr(i+1,j  ,k)
+                                          + z_phys_arr(i  ,j+1,k) + z_phys_arr(i+1,j+1,k) );
+                z_pos_arr(i,j,k) = z_bot_face + delta_z;
 
                 // Destination position must be contained on the current process!
                 Real pos[] = {x_pos_arr(i,j,k),y_pos_arr(i,j,k),0.5*dx[2]};