MPI support for A_phi calculation in GR

dev/nix/kokkos.nix

@@ -10,6 +10,7 @@ let
   pversion = "4.6.01";
   compilerPkgs = {
     "HIP" = with pkgs.rocmPackages; [
+      llvm.rocm-merged-llvm
       rocm-core
       clr
       rocthrust
@@ -18,11 +19,13 @@ let
       rocm-smi
     ];
     "CUDA" = with pkgs.cudaPackages; [
+      llvmPackages_18.clang-tools
       cudatoolkit
       cuda_cudart
       pkgs.gcc13
     ];
     "NONE" = [
+      pkgs.llvmPackages_18.clang-tools
       pkgs.gcc13
     ];
   };
@@ -36,9 +39,7 @@ let
     "HIP" = [
       "-D Kokkos_ENABLE_HIP=ON"
       "-D Kokkos_ARCH_${getArch { }}=ON"
-      # remove leading AMD_
       "-D AMDGPU_TARGETS=${builtins.replaceStrings [ "amd_" ] [ "" ] (pkgs.lib.toLower (getArch { }))}"
-      "-D CMAKE_C_COMPILER=hipcc"
       "-D CMAKE_CXX_COMPILER=hipcc"
     ];
     "CUDA" = [

dev/nix/shell.nix

@@ -42,9 +42,6 @@ pkgs.mkShell {
     zlib
     cmake
 
-    llvmPackages_18.clang-tools
-    libgcc
-
     adios2Pkg
     kokkosPkg
 

pgens/wald/pgen.hpp

@@ -7,6 +7,7 @@
 #include "arch/kokkos_aliases.h"
 #include "arch/traits.h"
 #include "utils/comparators.h"
+#include "utils/error.h"
 #include "utils/formatting.h"
 #include "utils/log.h"
 #include "utils/numeric.h"
@@ -17,14 +18,31 @@
 #include "framework/domain/domain.h"
 #include "framework/domain/metadomain.h"
 
+#include <string>
 #include <vector>
 
+enum InitFieldGeometry {
+  Wald,
+  Vertical,
+};
+
 namespace user {
   using namespace ntt;
 
   template <class M, Dimension D>
   struct InitFields {
-    InitFields(M metric_) : metric { metric_ } {}
+    InitFields(M metric_, const std::string& init_field_geometry)
+      : metric { metric_ } {
+      if (init_field_geometry == "wald") {
+        field_geometry = InitFieldGeometry::Wald;
+      } else if (init_field_geometry == "vertical") {
+        field_geometry = InitFieldGeometry::Vertical;
+      } else {
+        raise::Error(fmt::format("Unrecognized field geometry: %s",
+                                 init_field_geometry.c_str()),
+                     HERE);
+      }
+    }
 
     Inline auto A_3(const coord_t<D>& x_Cd) const -> real_t {
       return HALF * (metric.template h_<3, 3>(x_Cd) +
@@ -83,102 +101,132 @@ namespace user {
 
     Inline auto bx3(
       const coord_t<D>& x_Ph) const -> real_t { // at ( i + HALF , j + HALF )
-      coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
-      metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
-
-      x0m[0] = xi[0];
-      x0m[1] = xi[1] - HALF;
-      x0p[0] = xi[0];
-      x0p[1] = xi[1] + HALF;
-
-      real_t inv_sqrt_detH_iPjP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
-      return -(A_1(x0p) - A_1(x0m)) * inv_sqrt_detH_iPjP;
+      if (field_geometry == InitFieldGeometry::Wald) {
+        coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
+        metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
+
+        x0m[0] = xi[0];
+        x0m[1] = xi[1] - HALF;
+        x0p[0] = xi[0];
+        x0p[1] = xi[1] + HALF;
+
+        real_t inv_sqrt_detH_iPjP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
+        return -(A_1(x0p) - A_1(x0m)) * inv_sqrt_detH_iPjP;
+      } else if (field_geometry == InitFieldGeometry::Vertical) {
+        return ZERO;
+      } else {
+        raise::KernelError(HERE, "Unrecognized field geometry");
+        return ZERO;
+      }
     }
 
     Inline auto dx1(const coord_t<D>& x_Ph) const -> real_t { // at ( i + HALF , j )
-      coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
-      metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
-
-      real_t alpha_iPj { metric.alpha({ xi[0], xi[1] }) };
-      real_t inv_sqrt_detH_ij { ONE / metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) };
-      real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) };
-      real_t beta_ij { metric.beta1({ xi[0] - HALF, xi[1] }) };
-      real_t alpha_ij { metric.alpha({ xi[0] - HALF, xi[1] }) };
-
-      // D1 at ( i + HALF , j )
-      x0m[0] = xi[0] - HALF;
-      x0m[1] = xi[1];
-      x0p[0] = xi[0] + HALF;
-      x0p[1] = xi[1];
-      real_t E1d { (A_0(x0p) - A_0(x0m)) };
-      real_t D1d { E1d / alpha_iPj };
-
-      // D3 at ( i , j )
-      x0m[0] = xi[0] - HALF - HALF;
-      x0m[1] = xi[1];
-      x0p[0] = xi[0] - HALF + HALF;
-      x0p[1] = xi[1];
-      real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij };
-
-      real_t D1u { metric.template h<1, 1>({ xi[0], xi[1] }) * D1d +
-                   metric.template h<1, 3>({ xi[0], xi[1] }) * D3d };
-
-      return D1u;
+      if (field_geometry == InitFieldGeometry::Wald) {
+        coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
+        metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
+
+        real_t alpha_iPj { metric.alpha({ xi[0], xi[1] }) };
+        real_t inv_sqrt_detH_ij { ONE / metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) };
+        real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) };
+        real_t beta_ij { metric.beta1({ xi[0] - HALF, xi[1] }) };
+        real_t alpha_ij { metric.alpha({ xi[0] - HALF, xi[1] }) };
+
+        // D1 at ( i + HALF , j )
+        x0m[0] = xi[0] - HALF;
+        x0m[1] = xi[1];
+        x0p[0] = xi[0] + HALF;
+        x0p[1] = xi[1];
+        real_t E1d { (A_0(x0p) - A_0(x0m)) };
+        real_t D1d { E1d / alpha_iPj };
+
+        // D3 at ( i , j )
+        x0m[0] = xi[0] - HALF - HALF;
+        x0m[1] = xi[1];
+        x0p[0] = xi[0] - HALF + HALF;
+        x0p[1] = xi[1];
+        real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij };
+
+        real_t D1u { metric.template h<1, 1>({ xi[0], xi[1] }) * D1d +
+                     metric.template h<1, 3>({ xi[0], xi[1] }) * D3d };
+
+        return D1u;
+      } else if (field_geometry == InitFieldGeometry::Vertical) {
+        return ZERO;
+      } else {
+        raise::KernelError(HERE, "Unrecognized field geometry");
+        return ZERO;
+      }
     }
 
     Inline auto dx2(const coord_t<D>& x_Ph) const -> real_t { // at ( i , j + HALF )
-      coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
-      metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
-      x0m[0] = xi[0];
-      x0m[1] = xi[1] - HALF;
-      x0p[0] = xi[0];
-      x0p[1] = xi[1] + HALF;
-      real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
-      real_t sqrt_detH_ijP { metric.sqrt_det_h({ xi[0], xi[1] }) };
-      real_t alpha_ijP { metric.alpha({ xi[0], xi[1] }) };
-      real_t beta_ijP { metric.beta1({ xi[0], xi[1] }) };
-
-      real_t E2d { (A_0(x0p) - A_0(x0m)) };
-      real_t D2d { E2d / alpha_ijP - (A_1(x0p) - A_1(x0m)) * beta_ijP / alpha_ijP };
-      real_t D2u { metric.template h<2, 2>({ xi[0], xi[1] }) * D2d };
-
-      return D2u;
+      if (field_geometry == InitFieldGeometry::Wald) {
+        coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
+        metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
+        x0m[0] = xi[0];
+        x0m[1] = xi[1] - HALF;
+        x0p[0] = xi[0];
+        x0p[1] = xi[1] + HALF;
+        real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
+        real_t sqrt_detH_ijP { metric.sqrt_det_h({ xi[0], xi[1] }) };
+        real_t alpha_ijP { metric.alpha({ xi[0], xi[1] }) };
+        real_t beta_ijP { metric.beta1({ xi[0], xi[1] }) };
+
+        real_t E2d { (A_0(x0p) - A_0(x0m)) };
+        real_t D2d { E2d / alpha_ijP -
+                     (A_1(x0p) - A_1(x0m)) * beta_ijP / alpha_ijP };
+        real_t D2u { metric.template h<2, 2>({ xi[0], xi[1] }) * D2d };
+
+        return D2u;
+      } else if (field_geometry == InitFieldGeometry::Vertical) {
+        return ZERO;
+      } else {
+        raise::KernelError(HERE, "Unrecognized field geometry");
+        return ZERO;
+      }
     }
 
     Inline auto dx3(const coord_t<D>& x_Ph) const -> real_t { // at ( i , j )
-      coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
-      metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
-      real_t inv_sqrt_detH_ij { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
-      real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0], xi[1] }) };
-      real_t beta_ij { metric.beta1({ xi[0], xi[1] }) };
-      real_t alpha_ij { metric.alpha({ xi[0], xi[1] }) };
-      real_t alpha_iPj { metric.alpha({ xi[0] + HALF, xi[1] }) };
-
-      // D3 at ( i , j )
-      x0m[0] = xi[0] - HALF;
-      x0m[1] = xi[1];
-      x0p[0] = xi[0] + HALF;
-      x0p[1] = xi[1];
-      real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij };
-
-      // D1 at ( i + HALF , j )
-      x0m[0] = xi[0] + HALF - HALF;
-      x0m[1] = xi[1];
-      x0p[0] = xi[0] + HALF + HALF;
-      x0p[1] = xi[1];
-      real_t E1d { (A_0(x0p) - A_0(x0m)) };
-      real_t D1d { E1d / alpha_iPj };
-
-      if (cmp::AlmostZero(x_Ph[1])) {
-        return metric.template h<1, 3>({ xi[0], xi[1] }) * D1d;
+      if (field_geometry == InitFieldGeometry::Wald) {
+        coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
+        metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
+        real_t inv_sqrt_detH_ij { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
+        real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0], xi[1] }) };
+        real_t beta_ij { metric.beta1({ xi[0], xi[1] }) };
+        real_t alpha_ij { metric.alpha({ xi[0], xi[1] }) };
+        real_t alpha_iPj { metric.alpha({ xi[0] + HALF, xi[1] }) };
+
+        // D3 at ( i , j )
+        x0m[0] = xi[0] - HALF;
+        x0m[1] = xi[1];
+        x0p[0] = xi[0] + HALF;
+        x0p[1] = xi[1];
+        real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij };
+
+        // D1 at ( i + HALF , j )
+        x0m[0] = xi[0] + HALF - HALF;
+        x0m[1] = xi[1];
+        x0p[0] = xi[0] + HALF + HALF;
+        x0p[1] = xi[1];
+        real_t E1d { (A_0(x0p) - A_0(x0m)) };
+        real_t D1d { E1d / alpha_iPj };
+
+        if (cmp::AlmostZero(x_Ph[1])) {
+          return metric.template h<1, 3>({ xi[0], xi[1] }) * D1d;
+        } else {
+          return metric.template h<3, 3>({ xi[0], xi[1] }) * D3d +
+                 metric.template h<1, 3>({ xi[0], xi[1] }) * D1d;
+        }
+      } else if (field_geometry == InitFieldGeometry::Vertical) {
+        return ZERO;
       } else {
-        return metric.template h<3, 3>({ xi[0], xi[1] }) * D3d +
-               metric.template h<1, 3>({ xi[0], xi[1] }) * D1d;
+        raise::KernelError(HERE, "Unrecognized field geometry");
+        return ZERO;
       }
     }
 
   private:
-    const M metric;
+    const M           metric;
+    InitFieldGeometry field_geometry;
   };
 
   template <SimEngine::type S, class M>
@@ -201,7 +249,8 @@ namespace user {
     inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
       : arch::ProblemGenerator<S, M> { p }
       , global_domain { m }
-      , init_flds { m.mesh().metric } {}
+      , init_flds { m.mesh().metric,
+                    p.template get<std::string>("setup.init_field", "wald") } {}
   };
 
 } // namespace user

pgens/wald/wald.toml

@@ -4,11 +4,11 @@
   runtime = 100.0
 
 [grid]
-  resolution = [128, 128]
+  resolution = [512, 512]
   extent     = [[1.0, 10.0]]
 
   [grid.metric]
-    metric  = "kerr_schild"
+    metric  = "qkerr_schild"
     qsph_r0 = 0.0
     qsph_h  = 0.0
     ks_a    = 0.95
@@ -38,10 +38,10 @@
   ppc0 = 2.0
 
 [setup]
+  init_field = "wald" # or "vertical"
 
 [output]
-  format         = "hdf5"
-  separate_files = false
+  format = "hdf5"
 
   [output.fields]
     interval_time = 1.0

src/framework/domain/metadomain.h

@@ -261,6 +261,9 @@ namespace ntt {
 #if defined(OUTPUT_ENABLED)
     out::Writer        g_writer;
     checkpoint::Writer g_checkpoint_writer;
+  #if defined(MPI_ENABLED)
+    void CommunicateVectorPotential(unsigned short);
+  #endif
 #endif
 
 #if defined(MPI_ENABLED)