Homme(xx)/SL: Finish C++/Kokkos for ETM; modify vertical discretization.

components/eamxx/cime_config/namelist_defaults_eamxx.xml

@@ -895,8 +895,11 @@ be lost if SCREAM_HACK_XML is not enabled.
     <!-- SL transport settings. SL defaults to on for pg2 configs. -->
     <dt_tracer_factor constraints="ge 1" hgrid=".*pg2">6</dt_tracer_factor>
     <hypervis_subcycle_q hgrid=".*pg2">6</hypervis_subcycle_q>
-    <transport_alg hgrid=".*pg2">12</transport_alg>
-    <semi_lagrange_trajectory_nsubstep>0</semi_lagrange_trajectory_nsubstep>
+    <transport_alg hgrid=".*pg2" doc="0 (Eulerian) or 12 (semi-Lagrangian) tracer transport">12</transport_alg>
+    <semi_lagrange_trajectory_nsubstep doc="Controls trajectory method. 0: original method. 1: New method. 2 or more: New method with this many trajectory substeps per tracer time step.">0</semi_lagrange_trajectory_nsubstep>
+    <semi_lagrange_trajectory_nvelocity doc="Number of velocity slices to use in new method. 2 or less maps to 2">-1</semi_lagrange_trajectory_nvelocity>
+    <semi_lagrange_halo doc="Max number of element halos available in communication. -1 triggers an automatic estimate.">-1</semi_lagrange_halo>
+    <semi_lagrange_diagnostics>0</semi_lagrange_diagnostics>
     <!-- Other settings that we'll trigger based on pg2 for convenience -->
     <se_ftype valid_values="0,2" hgrid=".*pg2">2</se_ftype>
     <mesh_file type="file">none</mesh_file>

components/eamxx/cime_config/testdefs/testmods_dirs/eamxx/sl_nsubstep2/shell_commands

@@ -1,3 +1,4 @@
 ATMCHANGE=$CIMEROOT/../components/eamxx/scripts/atmchange
 
 $ATMCHANGE semi_lagrange_trajectory_nsubstep=2 -b
+$ATMCHANGE semi_lagrange_trajectory_nvelocity=3 -b

components/homme/cmake/HommeMacros.cmake

@@ -829,6 +829,24 @@ MACRO(CREATE_CXX_VS_F90_TESTS_WITH_PROFILE TESTS_LIST testProfile)
   ENDFOREACH ()
 ENDMACRO(CREATE_CXX_VS_F90_TESTS_WITH_PROFILE)
 
+function(check_transport_error_norms
+    TEST_NAME TCEN_ERROR_ANCHOR TCEN_FILENAME TCEN_UPPER_BOUNDS)
+  # Check prescribed-wind tracer transport test error-norm output, which is
+  # available with some of these tests. This function encapsulates the setup
+  # steps to add a check after a tracer transport test runs. See
+  # TransportCheckErrorNorms.cmake.in for details.
+  configure_file(
+    ${HOMME_SOURCE_DIR}/cmake/TransportCheckErrorNorms.cmake.in
+    ${HOMME_BINARY_DIR}/tests/${TEST_NAME}/check.cmake
+    @ONLY)
+  add_test(
+    NAME "${TEST_NAME}_l2err"
+    COMMAND ${CMAKE_COMMAND} -P check.cmake
+    WORKING_DIRECTORY ${HOMME_BINARY_DIR}/tests/${TEST_NAME})
+  set_tests_properties(
+    "${TEST_NAME}_l2err" PROPERTIES DEPENDS "${TEST_NAME}")
+endfunction()
+
 macro(testQuadPrec HOMME_QUAD_PREC)
 
   TRY_COMPILE(COMPILE_RESULT_VAR

components/homme/cmake/TransportCheckErrorNorms.cmake.in

@@ -0,0 +1,39 @@
+# This utility checks output of the form
+#   planar_conv>                          l2                    linf
+#   planar_conv> Q 1  3.7612654153016604E-03  9.8620728873592828E-03
+#   planar_conv> Q 2  3.9526251584022527E-03  8.2639849437505641E-03
+# against caller-provided upper bounds on the l2 norms. Output of this sort
+# appears on stdout at the end of some tracer transport tests. 'planar_conv' can
+# be something else; it is caller-configurable. Example usage:
+#   set(TCEN_ERROR_ANCHOR "planar_conv")
+#   set(TCEN_FILENAME "out.txt")
+#   set(TCEN_UPPER_BOUNDS "3.9e-3;4.1e-3")
+#   configure_file(TransportCheckErrorNorms.cmake.in check.cmake @ONLY)
+#   add_test(
+#     NAME "${TEST_NAME}_l2err"
+#     COMMAND ${CMAKE_COMMAND} -P check.cmake
+#     WORKING_DIRECTORY ${HOMME_BINARY_DIR}/tests/${TEST_NAME})
+#   set_tests_properties(
+#     "${TEST_NAME}_l2err" PROPERTIES DEPENDS "${TEST_NAME}")
+# Function check_transport_error_norms in HommeMacros.cmake encapsulates the
+# above commands.
+
+function(check_error_norms test_name filename upper_bounds)
+  execute_process(COMMAND grep "${test_name}> Q" ${filename} OUTPUT_VARIABLE grep_stdout)
+  separate_arguments(tokens NATIVE_COMMAND "${grep_stdout}")
+  list(LENGTH upper_bounds nchk)
+  math(EXPR nchk "${nchk} - 1")
+  foreach(idx RANGE 0 ${nchk})
+    math(EXPR tidx "3 + 5 * ${idx}")
+    list(GET tokens ${tidx} e)
+    list(GET upper_bounds ${idx} ub)
+    if(e GREATER ub)
+      message(FATAL_ERROR "${e} should be <= ${ub} but is not")
+    endif()
+  endforeach()
+endfunction()
+
+check_error_norms(
+  "@TCEN_ERROR_ANCHOR@"
+  "@TCEN_FILENAME@"
+  "@TCEN_UPPER_BOUNDS@")

components/homme/cmake/machineFiles/aurora-aot.cmake

@@ -44,7 +44,8 @@ SET(CMAKE_Fortran_COMPILER "mpifort" CACHE STRING "")
 SET(CMAKE_CXX_COMPILER "mpicxx" CACHE STRING "")
 
 #AOT flags
-SET(SYCL_COMPILE_FLAGS "-std=c++17 -fsycl -fsycl-device-code-split=per_kernel -fno-sycl-id-queries-fit-in-int -fsycl-unnamed-lambda -Xclang -fsycl-allow-virtual-functions")
+#SET(SYCL_COMPILE_FLAGS "-std=c++17 -fsycl -fsycl-device-code-split=per_kernel -fno-sycl-id-queries-fit-in-int -fsycl-unnamed-lambda -Xclang -fsycl-allow-virtual-functions")
+SET(SYCL_COMPILE_FLAGS "-std=c++17 -fsycl -fsycl-device-code-split=per_kernel -fno-sycl-id-queries-fit-in-int -fsycl-unnamed-lambda")
 SET(SYCL_LINK_FLAGS "-Wl,--no-relax -flink-huge-device-code -fsycl-max-parallel-link-jobs=32 -fsycl -fsycl-device-code-split=per_kernel -fsycl-targets=intel_gpu_pvc")
 
 SET(ADD_Fortran_FLAGS "-fc=ifx -fpscomp logicals -O3 -DNDEBUG -DCPRINTEL -g" CACHE STRING "")

components/homme/cmake/machineFiles/chrysalis.cmake

@@ -54,6 +54,7 @@ IF (${IFORT_RESULT} EQUAL 0)
   SET (ADD_Fortran_FLAGS "-traceback" CACHE STRING "")
   SET (ADD_C_FLAGS       "-traceback" CACHE STRING "")
   SET (ADD_CXX_FLAGS     "-traceback" CACHE STRING "")
+  SET (BUILD_HOMME_THETA_KOKKOS TRUE CACHE BOOL "")
 ELSE()
   SET (MKLROOT $ENV{MKLROOT} CACHE FILEPATH "")
   SET (HOMME_FIND_BLASLAPACK TRUE CACHE BOOL "")

components/homme/src/preqx_kokkos/CMakeLists.txt

@@ -97,6 +97,7 @@ MACRO(PREQX_KOKKOS_SETUP)
     ${SRC_SHARE_DIR}/cxx/prim_cxx_driver_base.F90
     ${SRC_SHARE_DIR}/cxx/utilities/bfb_mod.F90
     ${SRC_SHARE_DIR}/planar_mod.F90
+    ${TEST_SRC_DIR}/planar_transport.F90
     ${SRC_SHARE_DIR}/geometry_mod.F90
     ${SRC_SHARE_DIR}/planar_mesh_mod.F90
   )

components/homme/src/share/compose/compose_hommexx.cpp

@@ -24,8 +24,8 @@ using View = typename TracerArrays<ko::MachineTraits>::View<DataType>;
 void set_views (const SetView<HommexxReal***>& spheremp,
                 const SetView<HommexxReal****>& dp, const SetView5& dp3d,
                 const SetView<HommexxReal******>& qdp, const SetView5& q,
-                const SetView5& dep_points, const SetView5& vnode,
-                const SetView5& vdep, const Int ndim) {
+                const SetView5& dep_points, const SetView5& vnode, const Int ndim,
+                const SetView5& vdep, const Int vdep_ndim) {
   static_assert(std::is_same<Real, HommexxReal>::value,
                 "Hommexx and Compose real types must be the same.");
 #ifdef COMPOSE_PORT
@@ -42,23 +42,22 @@ void set_views (const SetView<HommexxReal***>& spheremp,
   if (vnode.data())
     ta.vnode = View<Real****>(vnode.data(), nel, vnode.extent_int(1), np2, ndim);
   if (vdep.data())
-    ta.vdep  = View<Real****>(vdep.data(),  nel, vdep .extent_int(1), np2, ndim);
+    ta.vdep  = View<Real****>(vdep.data(),  nel, vdep .extent_int(1), np2, vdep_ndim);
 #else
   slmm_throw_if(true, "Running a Hommexx code path with the non-Hommexx build"
                 " is not supported.\n");
 #endif
 }
 
-void set_hvcoord (const HommexxReal etai_beg, const HommexxReal etai_end,
-                  const HommexxReal* etam) {
+void set_hvcoord (const HommexxReal* etai, const HommexxReal* etam) {
   auto& cm = *get_isl_mpi_singleton();
-  islmpi::set_hvcoord(cm, etai_beg, etai_end, etam);
+  islmpi::set_hvcoord(cm, etai, etam);
 }
 
-void calc_v_departure (const int step, const HommexxReal dtsub) {
+void interp_v_update (const int step, const HommexxReal dtsub) {
   auto& cm = *get_isl_mpi_singleton();
-  islmpi::calc_v_departure<>(cm, 0, cm.nelemd - 1, step, dtsub,
-                             nullptr, nullptr, nullptr);
+  islmpi::interp_v_update<>(cm, 0, cm.nelemd - 1, step, dtsub,
+                            nullptr, nullptr, nullptr);
 }
 
 void advect (const int np1, const int n0_qdp, const int np1_qdp) {

components/homme/src/share/compose/compose_hommexx.hpp

@@ -16,12 +16,11 @@ typedef SetView<HommexxReal*****> SetView5;
 void set_views(const SetView<HommexxReal***>& spheremp,
                const SetView<HommexxReal****>& dp, const SetView5& dp3d,
                const SetView<HommexxReal******>& qdp, const SetView5& q,
-               const SetView5& dep_points, const SetView5& vnode,
-               const SetView5& vdep, const int trajectory_ndim);
+               const SetView5& dep_points, const SetView5& vnode, const int ndim,
+               const SetView5& vdep, const int vdep_ndim);
 
-void set_hvcoord(const HommexxReal etai_beg, const HommexxReal etai_end,
-                 const HommexxReal* etam);
-void calc_v_departure(const int step, const HommexxReal dtsub);
+void set_hvcoord(const HommexxReal* etai, const HommexxReal* etam);
+void interp_v_update(const int step, const HommexxReal dtsub);
 
 void advect(const int np1, const int n0_qdp, const int np1_qdp);
 

components/homme/src/share/compose/compose_slmm.cpp

@@ -124,41 +124,47 @@ void finalize_init_phase (IslMpi<MT>& cm, typename IslMpi<MT>::Advecter& advecte
 }
 
 template <typename MT>
-void set_hvcoord (IslMpi<MT>& cm, const Real etai_beg, const Real etai_end,
-                  const Real* etam) {
+void set_hvcoord (IslMpi<MT>& cm, const Real* etai, const Real* etam) {
   if (cm.etam.size() > 0) return;
 #if defined COMPOSE_HORIZ_OPENMP
 # pragma omp barrier
 # pragma omp master
 #endif
   {
     slmm_assert(cm.nlev > 0);
-    cm.etai_beg = etai_beg;
-    cm.etai_end = etai_end;
+    cm.etai_beg = etai[0];
+    cm.etai_end = etai[cm.nlev];
+    cm.etai = typename IslMpi<MT>::template ArrayD<Real*>("etai", cm.nlev+1);
     cm.etam = typename IslMpi<MT>::template ArrayD<Real*>("etam", cm.nlev);
-    const auto h = ko::create_mirror_view(cm.etam);
+    const auto hi = ko::create_mirror_view(cm.etai);
+    const auto hm = ko::create_mirror_view(cm.etam);
+    for (int k = 0; k <= cm.nlev; ++k) {
+      hi(k) = etai[k];
+      slmm_assert(k == 0 or hi(k) > hi(k-1));
+      slmm_assert(hi(k) >= 0 && hi(k) <= 1);
+    }
     for (int k = 0; k < cm.nlev; ++k) {
-      h(k) = etam[k];
-      slmm_assert(k == 0 || h(k) > h(k-1));
-      slmm_assert(h(k) > 0 && h(k) < 1);
+      hm(k) = etam[k];
+      slmm_assert(k == 0 or hm(k) > hm(k-1));
+      slmm_assert(hm(k) > 0 && hm(k) < 1);
     }
-    ko::deep_copy(cm.etam, h);
+    ko::deep_copy(cm.etai, hi);
+    ko::deep_copy(cm.etam, hm);
   }
 #if defined COMPOSE_HORIZ_OPENMP
 # pragma omp barrier
 #endif
 }
 
 template void set_hvcoord(
-  IslMpi<ko::MachineTraits>& cm, const Real etai_beg, const Real etai_end,
-  const Real* etam);
+  IslMpi<ko::MachineTraits>& cm, const Real* etai, const Real* etam);
 
 // Set pointers to HOMME data arrays.
 template <typename MT>
 void set_elem_data (IslMpi<MT>& cm, const Int ie, Real* qdp, const Int n0_qdp,
                     const Real* dp, Real* q, const Int nelem_in_patch) {
   slmm_assert(ie < cm.ed_h.size());
-  slmm_assert(cm.halo > 1 || cm.ed_h(ie).nbrs.size() == nelem_in_patch);
+  slmm_assert(cm.halo > 1 or cm.ed_h(ie).nbrs.size() == nelem_in_patch);
   auto& e = cm.ed_h(ie);
 #if defined COMPOSE_PORT
   cm.tracer_arrays->pqdp.set_ie_ptr(ie, qdp);
@@ -395,15 +401,14 @@ void slmm_check_ref2sphere (homme::Int ie, homme::Cartesian3D* p) {
   amb::dev_fin_threads();
 }
 
-void slmm_set_hvcoord (const homme::Real etai_beg, const homme::Real etai_end,
-                       const homme::Real* etam) {
+void slmm_set_hvcoord (const homme::Real* etai, const homme::Real* etam) {
   amb::dev_init_threads();
   slmm_assert(homme::g_csl_mpi);
-  homme::islmpi::set_hvcoord(*homme::g_csl_mpi, etai_beg, etai_end, etam);
+  homme::islmpi::set_hvcoord(*homme::g_csl_mpi, etai, etam);
   amb::dev_fin_threads();
 }
 
-void slmm_calc_v_departure (
+void slmm_interp_v_update (
   homme::Int nets, homme::Int nete, homme::Int step, homme::Real dtsub,
   homme::Real* dep_points, homme::Int dep_points_ndim, homme::Real* vnode,
   homme::Real* vdep, homme::Int* info)
@@ -419,8 +424,8 @@ void slmm_calc_v_departure (
     homme::sl_traj_h2d(*cm.tracer_arrays, dep_points, vnode, vdep,
                        cm.dep_points_ndim);
   }
-  homme::islmpi::calc_v_departure(cm, nets - 1, nete - 1, step - 1,
-                                  dtsub, dep_points, vnode, vdep);
+  homme::islmpi::interp_v_update(cm, nets - 1, nete - 1, step - 1,
+                                 dtsub, dep_points, vnode, vdep);
   *info = 0;
   {
     slmm::Timer timer("d2h");

components/homme/src/share/compose/compose_slmm_islmpi.hpp

@@ -579,10 +579,10 @@ struct IslMpi {
   const typename Advecter::ConstPtr advecter;
   const Int np, np2, nlev, qsize, qsized, nelemd, halo;
   const bool traj_3d;
-  const Int traj_nsubstep, dep_points_ndim;
+  const Int traj_nsubstep, dep_points_ndim, traj_msg_sz;
 
   Real etai_beg, etai_end;
-  ArrayD<Real*> etam;
+  ArrayD<Real*> etai, etam;
 
   ElemDataListH ed_h; // this rank's owned cells, indexed by LID
   ElemDataListD ed_d;
@@ -637,8 +637,10 @@ struct IslMpi {
           Int itraj_3d, Int itraj_nsubstep)
     : p(ip), advecter(advecter),
       np(inp), np2(np*np), nlev(inlev), qsize(iqsize), qsized(iqsized), nelemd(inelemd),
-      halo(ihalo), traj_3d(itraj_3d), traj_nsubstep(itraj_nsubstep),
-      dep_points_ndim(traj_3d && traj_nsubstep > 0 ? 4 : 3),
+      halo(ihalo), traj_3d(itraj_3d),
+      traj_nsubstep(itraj_nsubstep),
+      dep_points_ndim(traj_3d and traj_nsubstep > 0 ? 4 : 3),
+      traj_msg_sz(traj_3d ? 5 : dep_points_ndim),
       tracer_arrays(itracer_arrays)
   {}
 
@@ -768,11 +770,10 @@ void step(
   Real* q_min_r, Real* q_max_r);
 
 template <typename MT = ko::MachineTraits>
-void set_hvcoord(IslMpi<MT>& cm, const Real etai_beg, const Real etai_end,
-                 const Real* etam);
+void set_hvcoord(IslMpi<MT>& cm, const Real* etai, const Real* etam);
 
 template <typename MT = ko::MachineTraits>
-void calc_v_departure(
+void interp_v_update(
   IslMpi<MT>& cm, const Int nets, const Int nete, const Int step, const Real dtsub,
   Real* dep_points_r, const Real* vnode, Real* vdep);