Merge 'CONRAD/investigate_adjust_ps' (PR #6063)

This is an incremental change to introduce adjustment of dp3d
instead of adjustment of ps.

adjust_ps is still used for preqx (forcing was never rewritten),
rsplit=0 (only adjust_ps is possible), or EAM (to not touch every test).

Fixed bugs in forcing_ut init, testing it only
with adjust_ps=false option.

[non-BFB] for HOMME tests with baselines (with moist forcing
and some without nontrivial forcing, because even zero forcing
is processed differently for adjustment of dp3d).

components/homme/src/share/cxx/HybridVCoord.cpp

@@ -158,8 +158,8 @@ void HybridVCoord::random_init(int seed) {
 
     Errors::runtime_check(curr>prev,"Error! hybrid_a+hybrid_b is not increasing.\n", -1);
 
-    host_hybrid_am_real(i-1) = (host_hybrid_ai(i) + host_hybrid_ai(i))/2.0;
-    host_hybrid_bm_real(i-1) = (host_hybrid_bi(i) + host_hybrid_bi(i))/2.0;
+    host_hybrid_am_real(i-1) = (host_hybrid_ai(i) + host_hybrid_ai(i-1))/2.0;
+    host_hybrid_bm_real(i-1) = (host_hybrid_bi(i) + host_hybrid_bi(i-1))/2.0;
   }
 
   Kokkos::deep_copy(hybrid_ai, host_hybrid_ai);

components/homme/src/share/prim_driver_base.F90

@@ -1613,16 +1613,16 @@ subroutine applyCAMforcing_tracers(elem,hvcoord,np1,np1_qdp,dt,adjustment)
 
 #ifdef MODEL_THETA_L
   if (dt_remap_factor==0) then
-     adjust_ps=.true.   ! stay on reference levels for Eulerian case
+    adjust_ps=.true.   ! stay on reference levels for Eulerian case
   else
-#ifdef SCREAM
-     adjust_ps=.false.  ! Lagrangian case can support adjusting dp3d or ps
-#else
-     adjust_ps=.true.   ! Lagrangian case can support adjusting dp3d or ps
-#endif
+    adjust_ps=.false.  ! Lagrangian case can support adjusting dp3d or ps
   endif
 #else
-  adjust_ps=.true.      ! preqx requires forcing to stay on reference levels
+  adjust_ps=.true.     ! preqx requires forcing to stay on reference levels
+#endif
+
+#ifdef CAM
+  adjust_ps=.true.     ! For CAM runs, require forcing to stay on reference levels
 #endif
 
   dp=elem%state%dp3d(:,:,:,np1)

components/homme/src/theta-l_kokkos/cxx/ElementsState.cpp

@@ -238,6 +238,80 @@ void ElementsState::randomize(const int seed,
   Kokkos::fence();
 }
 
+
+void ElementsState::randomize(const int seed,
+                              const HybridVCoord& hvcoord) {
+  // Check elements were inited
+  assert (m_num_elems>0);
+
+  // Check data makes sense
+  assert (hvcoord.ps0>0);
+  assert (hvcoord.hybrid_ai0>=0);
+
+  // Arbitrary minimum value to generate
+  constexpr const Real min_value = 0.015625;
+
+  std::mt19937_64 engine(seed);
+  std::uniform_real_distribution<Real> random_dist(min_value, 1.0 / min_value);
+  std::uniform_real_distribution<Real> pdf_vtheta_dp(100.0, 1000.0);
+
+  genRandArray(m_v,         engine, random_dist);
+  genRandArray(m_w_i,       engine, random_dist);
+  genRandArray(m_vtheta_dp, engine, pdf_vtheta_dp);
+  // Note: to avoid errors in the equation of state, we need phi to be increasing.
+  //       Rather than using a constraint (which may call the function many times,
+  //       we simply ask that there are no duplicates, then we sort it later.
+  auto sort_and_chek = [](const ExecViewManaged<Scalar[NUM_LEV_P]>::HostMirror v)->bool {
+    Real* start = reinterpret_cast<Real*>(v.data());
+    Real* end   = reinterpret_cast<Real*>(v.data()) + NUM_LEV_P*VECTOR_SIZE;
+    std::sort(start,end);
+    std::reverse(start,end);
+    auto it = std::unique(start,end);
+    return it==end;
+  };
+  for (int ie=0; ie<m_num_elems; ++ie) {
+    for (int itl=0; itl<NUM_TIME_LEVELS; ++itl) {
+      for (int igp=0; igp<NP; ++igp) {
+        for (int jgp=0; jgp<NP; ++ jgp) {
+          genRandArray(Homme::subview(m_phinh_i,ie,itl,igp,jgp),engine,random_dist,sort_and_chek);
+        }
+      }
+    }
+  }
+
+  std::uniform_real_distribution<Real> pressure_pdf(800, 1200);
+  genRandArray(m_ps_v, engine, pressure_pdf);
+
+  auto dp = m_dp3d;
+  auto ps = m_ps_v;
+  auto ps0 = hvcoord.ps0;
+  auto hyai = hvcoord.hybrid_ai_packed;
+  auto hybi = hvcoord.hybrid_bi_packed;
+  auto hyai_delta = hvcoord.hybrid_ai_delta;
+  auto hybi_delta = hvcoord.hybrid_bi_delta;
+  const auto tu = m_tu;
+  Kokkos::parallel_for(m_policy, KOKKOS_LAMBDA(const TeamMember& team) {
+    KernelVariables kv(team, tu);
+    const int ie = kv.ie / NUM_TIME_LEVELS;
+    const int tl = kv.ie % NUM_TIME_LEVELS;
+
+    Kokkos::parallel_for(Kokkos::TeamThreadRange(kv.team,NP*NP),
+                         [&](const int idx) {
+      const int igp  = idx / NP;
+      const int jgp  = idx % NP;
+      ColumnOps::compute_midpoint_delta(kv,hyai,hyai_delta);
+      ColumnOps::compute_midpoint_delta(kv,hybi,hybi_delta);
+      team.team_barrier();
+      Kokkos::parallel_for(Kokkos::ThreadVectorRange(kv.team,NUM_LEV),
+                           [&](const int ilev) {
+        dp(ie,tl,igp,jgp,ilev) = ps0*hyai_delta(ilev)
+                               + ps(ie,tl,igp,jgp)*hybi_delta(ilev);
+      });
+    });
+  });
+  Kokkos::fence();
+}
+
 void ElementsState::pull_from_f90_pointers (CF90Ptr& state_v,         CF90Ptr& state_w_i,
                                             CF90Ptr& state_vtheta_dp, CF90Ptr& state_phinh_i,
                                             CF90Ptr& state_dp3d,      CF90Ptr& state_ps_v) {

components/homme/src/theta-l_kokkos/cxx/ElementsState.hpp

@@ -63,10 +63,10 @@ class ElementsState {
   void randomize(const int seed);
   void randomize(const int seed, const Real max_pressure);
   void randomize(const int seed, const Real max_pressure, const Real ps0, const Real hyai0);
-
   void randomize(const int seed, const Real max_pressure, const Real ps0, const Real hyai0,
                  const ExecViewUnmanaged<const Real*[NP][NP]>& phis);
-
+  void randomize(const int seed, const HybridVCoord& hvcoord);
+
   KOKKOS_INLINE_FUNCTION
   int num_elems() const { return m_num_elems; }
 

components/homme/src/theta-l_kokkos/cxx/ForcingFunctor.hpp

@@ -99,8 +99,7 @@ class ForcingFunctor
     m_hydrostatic = p.theta_hydrostatic_mode;
     m_qsize = p.qsize;
 
-    // TODO: this may change, depending on the simulation params
-    m_adjust_ps = true;
+    m_adjust_ps = (p.dt_remap_factor == 0);
 
     m_eos.init(m_hydrostatic,m_hvcoord);
     m_elem_ops.init(m_hvcoord);
@@ -153,7 +152,7 @@ class ForcingFunctor
     constexpr int int_size = NP*NP*NUM_LEV_P*VECTOR_SIZE;
 
     // 3 persistent midlayers, 2 non-persistent midlayer, and 1 non-persistent interface
-    return mid_size*(nelems*4+nslots) + (m_hydrostatic ? int_size*nslots : 0);
+    return mid_size*(nelems*4+nslots) + int_size*nslots;
   }
 
   void init_buffers (const FunctorsBuffersManager& fbm) {
@@ -370,7 +369,7 @@ class ForcingFunctor
           });
           if (!m_adjust_ps) {
             Kokkos::parallel_for(Kokkos::ThreadVectorRange(kv.team,NUM_LEV),
-                                 [&](const int ilev) {
+                                [&](const int ilev) {
               dp_adj(ilev) = dp(ilev) + dp(ilev)*(fq(ilev)-q(ilev));
             });
           }
@@ -386,7 +385,7 @@ class ForcingFunctor
           });
           if (!m_adjust_ps) {
             Kokkos::parallel_for(Kokkos::ThreadVectorRange(kv.team,NUM_LEV),
-                                 [&](const int& ilev) {
+                                [&](const int& ilev) {
               dp_adj(ilev) = dp(ilev) + compute_fqdt_pack(ilev,fq,qdp);
             });
           }
@@ -472,14 +471,14 @@ class ForcingFunctor
         } else {
           // Compute hydrostatic p from dp. Store in exner, then add to pnh
           auto p_i = Homme::subview(m_pi_i,kv.team_idx,igp,jgp);
-          m_elem_ops.compute_hydrostatic_p(kv,dp,p_i,exner);
+          m_elem_ops.compute_hydrostatic_p(kv,dp_adj,p_i,exner);
           Kokkos::parallel_for(Kokkos::ThreadVectorRange(kv.team,NUM_LEV),
-                               [&](const int ilev) {
+                              [&](const int ilev) {
             pnh(ilev) += exner(ilev);
             dp(ilev) = dp_adj(ilev);
           });
         }
-        
+
         Kokkos::parallel_for(Kokkos::ThreadVectorRange(kv.team,NUM_LEV),
                              [&](const int ilev) {
           using namespace PhysicalConstants;

components/homme/src/theta-l_kokkos/cxx/LimiterFunctor.hpp

@@ -67,10 +67,10 @@ struct LimiterFunctor {
       , m_hvcoord(hvcoord)
       , m_state(elements.m_state)
       , m_geometry(elements.m_geometry)
-      , m_policy_dp3d_lim (Homme::get_default_team_policy<ExecSpace,TagDp3dLimiter>(m_num_elems))
-      , m_tu(m_policy_dp3d_lim)
       , m_dp3d_thresh(params.dp3d_thresh)
       , m_vtheta_thresh(params.vtheta_thresh)
+      , m_policy_dp3d_lim (Homme::get_default_team_policy<ExecSpace,TagDp3dLimiter>(m_num_elems))
+      , m_tu(m_policy_dp3d_lim)
   {
     m_np1 = -1;
   }

components/homme/test_execs/thetal_kokkos_ut/forcing_ut.cpp

@@ -60,7 +60,8 @@ TEST_CASE("forcing", "forcing") {
   // Init everything through singleton, which is what happens in normal runs
   auto& c = Context::singleton();
   auto& p = c.create<SimulationParams>();
-
+  p.dt_remap_factor = 1;
+
   auto& hv = c.create<HybridVCoord>();
   hv.random_init(seed);
 
@@ -165,7 +166,7 @@ TEST_CASE("forcing", "forcing") {
           std::cout << "     -> adjustment: " << (adjustment ? "true" : "false") << "\n";
 
           // Reset state, tracers, and forcing to the original random values
-          state.randomize(seed, 10*hv.ps0, hv.ps0, hv.hybrid_ai0);
+          state.randomize(seed, hv);
           tracers.randomize(seed,-1e-3,1e-3);
           forcing.randomize(seed);
 
@@ -289,7 +290,7 @@ TEST_CASE("forcing", "forcing") {
     // Reset state and forcing to the original random values
     std::cout << "Testing dynamics forcing.\n";
 
-    state.randomize(seed, 10*hv.ps0, hv.ps0, hv.hybrid_ai0);
+    state.randomize(seed, hv);
     forcing.randomize(seed);
 
     // Sync views