WIP: Bump Parthenon version to current develop and new B Cleanup w/GMG

.gitmodules

@@ -1,7 +1,7 @@
 [submodule "external/parthenon"]
 	path = external/parthenon
-	url = https://github.com/AFD-Illinois/parthenon.git
-	branch = kharma
+	url = https://github.com/parthenon-hpc-lab/parthenon.git
+	branch = develop
 [submodule "external/variant"]
 	path = external/variant
 	url = https://github.com/mpark/variant.git

CMakeLists.txt

@@ -16,6 +16,8 @@ set(PARTHENON_ENABLE_CPP17 ON CACHE BOOL "KHARMA Override")
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
 
 # Parthenon options
+# Coordinates
+set(PARTHENON_COORDINATES "GRCoordinates" CACHE STRING "KHARMA Override")
 set(PARTHENON_DISABLE_EXAMPLES ON CACHE BOOL "KHARMA Override")
 set(PARTHENON_LINT_DEFAULT OFF CACHE BOOL "KHARMA Override")
 # Attempt HDF5 compression, requires recent/standard HDF5. YMMV
@@ -32,6 +34,7 @@ set(ENABLE_COMPILER_WARNINGS OFF CACHE BOOL "KHARMA Override")
 # Don't require static HDF5, as conda and other environments don't provide it
 set(HDF5_USE_STATIC_LIBRARIES OFF CACHE BOOL "KHARMA Override")
 
+
 # Kokkos options
 set(Kokkos_ENABLE_CUDA_LAMBDA ON CACHE BOOL "KHARMA Override")
 set(Kokkos_ENABLE_CUDA_CONSTEXPR ON CACHE BOOL "KHARMA Override")

external/kokkos-kernels/KokkosBatched_Dot_Internal.hpp

@@ -5,6 +5,8 @@
 
 #include "KokkosBatched_Util.hpp"
 
+#define KOKKOS_IMPL_DO_NOT_USE_PRINTF(...) ::printf(__VA_ARGS__)
+
 namespace KokkosBatched {
 
 ///

external/kokkos-kernels/KokkosBatched_Util.hpp

@@ -772,8 +772,8 @@ KOKKOS_INLINE_FUNCTION auto subview_wrapper(ViewType v, IdxType1 i1,
 }
 template <class ViewType, class IdxType1>
 KOKKOS_INLINE_FUNCTION auto subview_wrapper(ViewType v, IdxType1 i1,
-                                            Kokkos::Impl::ALL_t i2,
-                                            Kokkos::Impl::ALL_t i3,
+                                            Kokkos::ALL_t i2,
+                                            Kokkos::ALL_t i3,
                                             const BatchLayout::Left &layout_tag,
                                             const Trans::Transpose) {
   auto sv_nt = subview_wrapper(v, i1, i3, i2, layout_tag);
@@ -805,7 +805,7 @@ KOKKOS_INLINE_FUNCTION auto subview_wrapper(
 }
 template <class ViewType, class IdxType1>
 KOKKOS_INLINE_FUNCTION auto subview_wrapper(
-    ViewType v, IdxType1 i1, Kokkos::Impl::ALL_t i2, Kokkos::Impl::ALL_t i3,
+    ViewType v, IdxType1 i1, Kokkos::ALL_t i2, Kokkos::ALL_t i3,
     const BatchLayout::Right &layout_tag, const Trans::Transpose &) {
   auto sv_nt = subview_wrapper(v, i1, i3, i2, layout_tag);
 

external/patches/parthenon-sycl.patch

@@ -0,0 +1,32 @@
+diff --git a/src/interface/swarm_device_context.hpp b/src/interface/swarm_device_context.hpp
+index e438dcf53..389a4df2b 100644
+--- a/src/interface/swarm_device_context.hpp
++++ b/src/interface/swarm_device_context.hpp
+@@ -65,10 +65,12 @@ class SwarmDeviceContext {
+     // Particle is on neither this block nor a neighboring block
+     if (i < 0 || i > 3 || ((j < 0 || j > 3) && ndim_ > 1) ||
+         ((k < 0 || k > 3) && ndim_ > 2)) {
++#ifndef KOKKOS_ENABLE_SYCL
+       printf("[%i] k = %i j = %i i = %i\n", n, k, j, i);
+       printf("x = %e [%e %e]\n", x, x_min_, x_max_);
+       printf("y = %e [%e %e]\n", y, y_min_, y_max_);
+       printf("z = %e [%e %e]\n", z, z_min_, z_max_);
++#endif
+       PARTHENON_FAIL("Particle neighbor indices out of bounds; particle has somehow "
+                      "moved beyond the halo of adjacent blocks which is not permitted.");
+     }
+diff --git a/src/utils/error_checking.hpp b/src/utils/error_checking.hpp
+index eec8dd5b7..890cd3c85 100644
+--- a/src/utils/error_checking.hpp
++++ b/src/utils/error_checking.hpp
+@@ -151,8 +151,10 @@ inline void require_throws(const char *const condition, std::stringstream const
+
+ [[noreturn]] KOKKOS_INLINE_FUNCTION void
+ fail(const char *const message, const char *const filename, int const linenumber) {
++#ifndef KOKKOS_ENABLE_SYCL
+   printf("### PARTHENON ERROR\n  Message:     %s\n  File:        %s\n  Line number: %i\n",
+          message, filename, linenumber);
++#endif
+   Kokkos::abort(message);
+   // Kokkos::abort ends control flow, but is not marked as `[[noreturn]]`, so we need this
+   // loop to supress a warning that the function does not return.

external/patches/parthenon-use-gr-coordinates.patch

@@ -1,17 +1,19 @@
 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
-index aaeabd8a..e354ef3d 100644
+index b5d275324..b648ea101 100644
 --- a/src/CMakeLists.txt
 +++ b/src/CMakeLists.txt
-@@ -90,7 +90,7 @@ set(COMPILED_WITH ${CMAKE_CXX_COMPILER})
+@@ -90,8 +90,8 @@ set(COMPILED_WITH ${CMAKE_CXX_COMPILER})
  set(COMPILER_COMMAND "<not-implemented>") # TODO: Put something more descriptive here
  set(COMPILER_FLAGS "<not-implemented>") # TODO: Put something more descriptive here
 
--set(COORDINATE_TYPE UniformCartesian) # TODO: Make this an option when more are available
-+set(COORDINATE_TYPE GRCoordinates) # TODO: Make this an option when more are available
-
- configure_file(config.hpp.in generated/config.hpp @ONLY)
-
-@@ -331,6 +331,8 @@ lint_target(parthenon)
+-set(COORDINATES_AVAILABLE UniformCartesian UniformCylindrical UniformSpherical)
+-set(PARTHENON_COORDINATES "UniformCartesian" CACHE STRING "Set the coordinate system from \{${COORDINATES_AVAILABLE}\} with UniformCartesian the default")
++set(COORDINATES_AVAILABLE UniformCartesian UniformCylindrical UniformSpherical GRCoordinates)
++set(PARTHENON_COORDINATES "GRCoordinates" CACHE STRING "Set the coordinate system from \{${COORDINATES_AVAILABLE}\} with UniformCartesian the default")
+ if (${PARTHENON_COORDINATES} IN_LIST COORDINATES_AVAILABLE)
+   set(COORDINATE_TYPE ${PARTHENON_COORDINATES})
+ else()
+@@ -386,6 +386,8 @@ lint_target(parthenon)
  target_include_directories(parthenon PUBLIC
    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
    $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/generated>
@@ -21,13 +23,13 @@ index aaeabd8a..e354ef3d 100644
    )
 
 diff --git a/src/coordinates/coordinates.hpp b/src/coordinates/coordinates.hpp
-index d1290dee..50bfc840 100644
+index 475232b86..88a422614 100644
 --- a/src/coordinates/coordinates.hpp
 +++ b/src/coordinates/coordinates.hpp
-@@ -16,6 +16,7 @@
- #include "config.hpp"
-
+@@ -18,6 +18,7 @@
  #include "uniform_cartesian.hpp"
+ #include "uniform_cylindrical.hpp"
+ #include "uniform_spherical.hpp"
 +#include "coordinates/gr_coordinates.hpp"
 
  namespace parthenon {

kharma/CMakeLists.txt

@@ -21,6 +21,7 @@ AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/flux EXE_NAME_SRC)
 
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/b_cd EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/b_cleanup EXE_NAME_SRC)
+AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/b_cleanup_gmg EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/b_ct EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/b_flux_ct EXE_NAME_SRC)
 AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR}/boundaries EXE_NAME_SRC)
@@ -48,6 +49,7 @@ include_directories(${CMAKE_CURRENT_SOURCE_DIR}/flux)
 
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/b_cd)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/b_cleanup)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/b_cleanup_gmg)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/b_ct)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/b_flux_ct)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/boundaries)

kharma/b_cd/b_cd.cpp

@@ -254,10 +254,10 @@ void FillOutput(MeshBlock *pmb, ParameterInput *pin)
     const int ndim = pmb->pmy_mesh->ndim;
     if (ndim < 2) return;
 
-    GridVector F1, F2, F3;
-    F1 = rc->Get("cons.B").flux[X1DIR];
-    F2 = rc->Get("cons.B").flux[X2DIR];
-    if (ndim > 2) F3 = rc->Get("cons.B").flux[X3DIR];
+    // GridVector F1, F2, F3;
+    // F1 = rc->Get("cons.B").flux[X1DIR];
+    // F2 = rc->Get("cons.B").flux[X2DIR];
+    // if (ndim > 2) F3 = rc->Get("cons.B").flux[X3DIR];
     auto& divB = rc->Get("divB").data;
 
     const IndexRange ib = rc->GetBoundsI(IndexDomain::interior);
@@ -272,11 +272,11 @@ void FillOutput(MeshBlock *pmb, ParameterInput *pin)
 
     pmb->par_for("B_field_bsqmax", kl.s, kl.e, jl.s, jl.e, il.s, il.e,
         KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
-            double divb_local = ((F1(V1, k, j, i+1) - F1(V1, k, j, i)) / G.Dxc<1>(i) +
-                                 (F2(V2, k, j+1, i) - F2(V2, k, j, i)) / G.Dxc<2>(j));
-            if (ndim > 2) divb_local += (F3(V3, k+1, j, i) - F3(V3, k, j, i)) / G.Dxc<3>(k);
+            // double divb_local = ((F1(V1, k, j, i+1) - F1(V1, k, j, i)) / G.Dxc<1>(i) +
+            //                      (F2(V2, k, j+1, i) - F2(V2, k, j, i)) / G.Dxc<2>(j));
+            // if (ndim > 2) divb_local += (F3(V3, k+1, j, i) - F3(V3, k, j, i)) / G.Dxc<3>(k);
 
-            divB(k, j, i) = divb_local;
+            divB(k, j, i) = 0.;
         }
     );
 }

kharma/b_cleanup/b_cleanup.cpp

@@ -50,11 +50,19 @@
 
 // The package should never be loaded if there is not a global solve to be done.
 // Therefore we yell at load time rather than waiting for the first solve
+// But we still need some stubs to compile
 std::shared_ptr<KHARMAPackage> B_Cleanup::Initialize(ParameterInput *pin, std::shared_ptr<Packages_t>& packages)
-{throw std::runtime_error("KHARMA was compiled without global solvers!  Cannot clean B Field!");}
-// We still need a stub for CleanupDivergence() in order to compile, but it will never be called
-TaskStatus B_Cleanup::CleanupDivergence(std::shared_ptr<MeshData<Real>>& md) {}
-bool B_Cleanup::CleanupThisStep(Mesh* pmesh, int nstep) {}
+{
+    throw std::runtime_error("KHARMA was compiled without global solvers!  Cannot clean B Field!");
+}
+TaskStatus B_Cleanup::CleanupDivergence(std::shared_ptr<MeshData<Real>>& md)
+{
+    throw std::runtime_error("KHARMA was compiled without global solvers!  Cannot clean B Field!");
+}
+bool B_Cleanup::CleanupThisStep(Mesh* pmesh, int nstep)
+{
+    return false;
+}
 
 #else
 
@@ -222,7 +230,7 @@ TaskStatus B_Cleanup::CleanupDivergence(std::shared_ptr<MeshData<Real>>& md)
     // Add a solver container as a shallow copy on the default MeshData
     // msolve is just a sub-set of vars we need from md, making MPI syncs etc faster
     std::vector<std::string> names = KHARMA::GetVariableNames(&pmesh->packages, {Metadata::GetUserFlag("B_Cleanup"), Metadata::GetUserFlag("StartupOnly")});
-    auto &msolve = pmesh->mesh_data.AddShallow("solve", names);
+    auto &msolve = pmesh->mesh_data.AddShallow("solve", md, names);
 
     // Initialize the divB variable, which we'll be solving against.
     if (use_b_ct) {
@@ -257,9 +265,9 @@ TaskStatus B_Cleanup::CleanupDivergence(std::shared_ptr<MeshData<Real>>& md)
     // execute it to perform BiCGStab iteration
     TaskID t_none(0);
     TaskCollection tc;
-    auto tr = tc.AddRegion(1);
+    auto& tr = tc.AddRegion(1);
     auto t_solve_step = solver.CreateTaskList(t_none, 0, tr, msolve, msolve);
-    while (!tr.Execute());
+    tc.Execute();
     // Make sure solution's ghost zones are sync'd
     KHARMADriver::SyncAllBounds(msolve);
 
@@ -541,4 +549,4 @@ TaskStatus B_Cleanup::CenterLaplacian(MeshData<Real>* md, const std::string& p_v
     return TaskStatus::complete;
 }
 
-#endif
+#endif // DISABLE_CLEANUP

kharma/b_cleanup/solvers.cpp

@@ -11,6 +11,10 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+#if DISABLE_CLEANUP
+
+#else
+
 #include "bicgstab_solver.hpp"
 
 namespace parthenon {
@@ -20,3 +24,5 @@ int BiCGStabCounter::global_num_bicgstab_solvers = 0;
 
 } // namespace solvers
 } // namespace parthenon
+
+#endif // DISABLE_CLEANUP