diff --git a/components/homme/src/share/cxx/HybridVCoord.cpp b/components/homme/src/share/cxx/HybridVCoord.cpp index 569eef1d6849..199b6bc4f739 100644 --- a/components/homme/src/share/cxx/HybridVCoord.cpp +++ b/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); diff --git a/components/homme/src/share/prim_driver_base.F90 b/components/homme/src/share/prim_driver_base.F90 index 27df23dba4d4..df83e1cc386b 100644 --- a/components/homme/src/share/prim_driver_base.F90 +++ b/components/homme/src/share/prim_driver_base.F90 @@ -1613,16 +1613,17 @@ 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 + +#if defined(CAM) && !defined(SCREAM) + adjust_ps=.true. ! Special case when CAM is defined, and SCREAM is not defined, + ! require forcing to stay on reference levels no matter dt_remap_factor #endif dp=elem%state%dp3d(:,:,:,np1) diff --git a/components/homme/src/theta-l_kokkos/cxx/ElementsState.cpp b/components/homme/src/theta-l_kokkos/cxx/ElementsState.cpp index f4ea57b71ca5..e7fe7f569c7d 100644 --- a/components/homme/src/theta-l_kokkos/cxx/ElementsState.cpp +++ b/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 random_dist(min_value, 1.0 / min_value); + std::uniform_real_distribution 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::HostMirror v)->bool { + Real* start = reinterpret_cast(v.data()); + Real* end = reinterpret_cast(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 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) { diff --git a/components/homme/src/theta-l_kokkos/cxx/ElementsState.hpp b/components/homme/src/theta-l_kokkos/cxx/ElementsState.hpp index 4fa353141402..aec7701661fa 100644 --- a/components/homme/src/theta-l_kokkos/cxx/ElementsState.hpp +++ b/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& phis); - + void randomize(const int seed, const HybridVCoord& hvcoord); + KOKKOS_INLINE_FUNCTION int num_elems() const { return m_num_elems; } diff --git a/components/homme/src/theta-l_kokkos/cxx/ForcingFunctor.hpp b/components/homme/src/theta-l_kokkos/cxx/ForcingFunctor.hpp index 80993d1d0f1d..28a702c1d273 100644 --- a/components/homme/src/theta-l_kokkos/cxx/ForcingFunctor.hpp +++ b/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; diff --git a/components/homme/src/theta-l_kokkos/cxx/LimiterFunctor.hpp b/components/homme/src/theta-l_kokkos/cxx/LimiterFunctor.hpp index bf93be710e94..cd3bf7c32526 100644 --- a/components/homme/src/theta-l_kokkos/cxx/LimiterFunctor.hpp +++ b/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(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(m_num_elems)) + , m_tu(m_policy_dp3d_lim) { m_np1 = -1; } diff --git a/components/homme/test_execs/thetal_kokkos_ut/forcing_ut.cpp b/components/homme/test_execs/thetal_kokkos_ut/forcing_ut.cpp index 7bfeeab0d748..5e4c51c7ca11 100644 --- a/components/homme/test_execs/thetal_kokkos_ut/forcing_ut.cpp +++ b/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(); - + p.dt_remap_factor = 1; + auto& hv = c.create(); 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