Z4c: Implement cell-centered spacetime

EinsteinToolkit · Nov 2, 2023 · 796f427 · 796f427
1 parent b725206
commit 796f427
Show file tree

Hide file tree

Showing 10 changed files with 45 additions and 44 deletions.
diff --git a/Z4c/interface.ccl b/Z4c/interface.ccl
@@ -14,6 +14,7 @@ USES INCLUDE HEADER: sum.hxx
 USES INCLUDE HEADER: vec.hxx
 USES INCLUDE HEADER: vect.hxx
 
+USES INCLUDE HEADER: newradx.hxx
 
 
 CCTK_INT FUNCTION GetCallFunctionCount()
@@ -23,29 +24,29 @@ REQUIRES FUNCTION GetCallFunctionCount
 
 # All variables have been shifted so that they tend to zero in flat space
 
-CCTK_REAL chi TYPE=gf TAGS='rhs="chi_rhs" dependents="ADMBaseX::metric"' "chi"
-CCTK_REAL gamma_tilde TYPE=gf TAGS='parities={+1 +1 +1   -1 -1 +1   -1 +1 -1   +1 +1 +1   +1 -1 -1   +1 +1 +1} rhs="gamma_tilde_rhs" dependents="ADMBaseX::metric"' { gammatxx gammatxy gammatxz gammatyy gammatyz gammatzz } "gamma-tilde"
-CCTK_REAL K_hat TYPE=gf TAGS='rhs="K_hat_rhs" dependents="ADMBaseX::curv"' { Kh } "K-hat"
-CCTK_REAL A_tilde TYPE=gf TAGS='parities={+1 +1 +1   -1 -1 +1   -1 +1 -1   +1 +1 +1   +1 -1 -1   +1 +1 +1} rhs="A_tilde_rhs" dependents="ADMBaseX::curv"' { Atxx Atxy Atxz Atyy Atyz Atzz } "A-tilde"
-CCTK_REAL Gam_tilde TYPE=gf TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} rhs="Gam_tilde_rhs"' { Gamtx Gamty Gamtz } "Gamma-tilde"
-CCTK_REAL Theta TYPE=gf TAGS='rhs="Theta_rhs" dependents="ADMBaseX::curv"' "Theta"
-CCTK_REAL alphaG TYPE=gf TAGS='rhs="alphaG_rhs" dependents="ADMBaseX::lapse ADMBaseX::dtlapse"' "alpha"
-CCTK_REAL betaG TYPE=gf TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} rhs="betaG_rhs" dependents="ADMBaseX::shift ADMBaseX::dtshift"' { betaGx betaGy betaGz } "beta"
+CCTK_REAL chi TYPE=gf CENTERING={ccc} TAGS='rhs="chi_rhs" dependents="ADMBaseX::metric"' "chi"
+CCTK_REAL gamma_tilde TYPE=gf CENTERING={ccc} TAGS='parities={+1 +1 +1   -1 -1 +1   -1 +1 -1   +1 +1 +1   +1 -1 -1   +1 +1 +1} rhs="gamma_tilde_rhs" dependents="ADMBaseX::metric"' { gammatxx gammatxy gammatxz gammatyy gammatyz gammatzz } "gamma-tilde"
+CCTK_REAL K_hat TYPE=gf CENTERING={ccc} TAGS='rhs="K_hat_rhs" dependents="ADMBaseX::curv"' { Kh } "K-hat"
+CCTK_REAL A_tilde TYPE=gf CENTERING={ccc} TAGS='parities={+1 +1 +1   -1 -1 +1   -1 +1 -1   +1 +1 +1   +1 -1 -1   +1 +1 +1} rhs="A_tilde_rhs" dependents="ADMBaseX::curv"' { Atxx Atxy Atxz Atyy Atyz Atzz } "A-tilde"
+CCTK_REAL Gam_tilde TYPE=gf CENTERING={ccc} TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} rhs="Gam_tilde_rhs"' { Gamtx Gamty Gamtz } "Gamma-tilde"
+CCTK_REAL Theta TYPE=gf CENTERING={ccc} TAGS='rhs="Theta_rhs" dependents="ADMBaseX::curv"' "Theta"
+CCTK_REAL alphaG TYPE=gf CENTERING={ccc} TAGS='rhs="alphaG_rhs" dependents="ADMBaseX::lapse ADMBaseX::dtlapse"' "alpha"
+CCTK_REAL betaG TYPE=gf CENTERING={ccc} TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} rhs="betaG_rhs" dependents="ADMBaseX::shift ADMBaseX::dtshift"' { betaGx betaGy betaGz } "beta"
 
 
 
-CCTK_REAL ZtC TYPE=gf TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} checkpoint="no"' { ZtCx ZtCy ZtCz } "Z-tilde"
-CCTK_REAL HC TYPE=gf TAGS='checkpoint="no"' "H"
-CCTK_REAL MtC TYPE=gf TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} checkpoint="no"' { MtCx MtCy MtCz } "M-tilde"
-CCTK_REAL allC TYPE=gf TAGS='checkpoint="no"' "constraint monitor"
+CCTK_REAL ZtC TYPE=gf CENTERING={ccc} TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} checkpoint="no"' { ZtCx ZtCy ZtCz } "Z-tilde"
+CCTK_REAL HC TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "H"
+CCTK_REAL MtC TYPE=gf CENTERING={ccc} TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} checkpoint="no"' { MtCx MtCy MtCz } "M-tilde"
+CCTK_REAL allC TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "constraint monitor"
 
 
 
-CCTK_REAL chi_rhs TYPE=gf TAGS='checkpoint="no"' "chi"
-CCTK_REAL gamma_tilde_rhs TYPE=gf TAGS='parities={+1 +1 +1   -1 -1 +1   -1 +1 -1   +1 +1 +1   +1 -1 -1   +1 +1 +1} checkpoint="no"' { gammatxx_rhs gammatxy_rhs gammatxz_rhs gammatyy_rhs gammatyz_rhs gammatzz_rhs } "gamma-tilde"
-CCTK_REAL K_hat_rhs TYPE=gf TAGS='checkpoint="no"' { Kh_rhs } "K-hat"
-CCTK_REAL A_tilde_rhs TYPE=gf TAGS='parities={+1 +1 +1   -1 -1 +1   -1 +1 -1   +1 +1 +1   +1 -1 -1   +1 +1 +1} checkpoint="no"' { Atxx_rhs Atxy_rhs Atxz_rhs Atyy_rhs Atyz_rhs Atzz_rhs } "A-tilde"
-CCTK_REAL Gam_tilde_rhs TYPE=gf TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} checkpoint="no"' { Gamtx_rhs Gamty_rhs Gamtz_rhs } "Gamma-tilde"
-CCTK_REAL Theta_rhs TYPE=gf TAGS='checkpoint="no"' "Theta"
-CCTK_REAL alphaG_rhs TYPE=gf TAGS='checkpoint="no"' "alpha"
-CCTK_REAL betaG_rhs TYPE=gf TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} checkpoint="no"' { betaGx_rhs betaGy_rhs betaGz_rhs } "beta"
+CCTK_REAL chi_rhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "chi"
+CCTK_REAL gamma_tilde_rhs CENTERING={ccc} TYPE=gf TAGS='parities={+1 +1 +1   -1 -1 +1   -1 +1 -1   +1 +1 +1   +1 -1 -1   +1 +1 +1} checkpoint="no"' { gammatxx_rhs gammatxy_rhs gammatxz_rhs gammatyy_rhs gammatyz_rhs gammatzz_rhs } "gamma-tilde"
+CCTK_REAL K_hat_rhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { Kh_rhs } "K-hat"
+CCTK_REAL A_tilde_rhs TYPE=gf CENTERING={ccc} TAGS='parities={+1 +1 +1   -1 -1 +1   -1 +1 -1   +1 +1 +1   +1 -1 -1   +1 +1 +1} checkpoint="no"' { Atxx_rhs Atxy_rhs Atxz_rhs Atyy_rhs Atyz_rhs Atzz_rhs } "A-tilde"
+CCTK_REAL Gam_tilde_rhs TYPE=gf CENTERING={ccc} TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} checkpoint="no"' { Gamtx_rhs Gamty_rhs Gamtz_rhs } "Gamma-tilde"
+CCTK_REAL Theta_rhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "Theta"
+CCTK_REAL alphaG_rhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "alpha"
+CCTK_REAL betaG_rhs TYPE=gf CENTERING={ccc} TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} checkpoint="no"' { betaGx_rhs betaGy_rhs betaGz_rhs } "beta"
diff --git a/Z4c/src/adm.cxx b/Z4c/src/adm.cxx
@@ -24,7 +24,7 @@ extern "C" void Z4c_ADM(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTS_Z4c_ADM;
   DECLARE_CCTK_PARAMETERS;
 
-  const array<int, dim> indextype = {0, 0, 0};
+  const array<int, dim> indextype = {1, 1, 1};
   const GF3D2layout layout1(cctkGH, indextype);
 
   const GF3D2<const CCTK_REAL> gf_chi1(layout1, chi);
@@ -91,7 +91,7 @@ extern "C" void Z4c_ADM(CCTK_ARGUMENTS) {
 #ifdef __CUDACC__
   const nvtxRangeId_t range = nvtxRangeStartA("Z4c_ADM::adm");
 #endif
-  grid.loop_all_device<0, 0, 0, vsize>(
+  grid.loop_all_device<1, 1, 1, vsize>(
       grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
         const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
         const GF3D2index index1(layout1, p.I);

diff --git a/Z4c/src/adm2.cxx b/Z4c/src/adm2.cxx
@@ -42,11 +42,11 @@ extern "C" void Z4c_ADM2(CCTK_ARGUMENTS) {
 
   //
 
-  const array<int, dim> indextype = {0, 0, 0};
+  const array<int, dim> indextype = {1, 1, 1};
   const array<int, dim> nghostzones = {cctk_nghostzones[0], cctk_nghostzones[1],
                                        cctk_nghostzones[2]};
   vect<int, dim> imin, imax;
-  GridDescBase(cctkGH).box_int<0, 0, 0>(nghostzones, imin, imax);
+  GridDescBase(cctkGH).box_int<1, 1, 1>(nghostzones, imin, imax);
   // Suffix 1: with ghost zones, suffix 0: without ghost zones
   const GF3D2layout layout1(cctkGH, indextype);
   const GF3D5layout layout0(imin, imax);
@@ -189,7 +189,7 @@ extern "C" void Z4c_ADM2(CCTK_ARGUMENTS) {
 #ifdef __CUDACC__
   const nvtxRangeId_t range = nvtxRangeStartA("Z4c_ADM2::adm2");
 #endif
-  grid.loop_int_device<0, 0, 0, vsize>(
+  grid.loop_int_device<1, 1, 1, vsize>(
       grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
         const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
         const GF3D2index index1(layout1, p.I);

diff --git a/Z4c/src/constraints.cxx b/Z4c/src/constraints.cxx
@@ -40,11 +40,11 @@ extern "C" void Z4c_Constraints(CCTK_ARGUMENTS) {
 
   //
 
-  const array<int, dim> indextype = {0, 0, 0};
+  const array<int, dim> indextype = {1, 1, 1};
   const array<int, dim> nghostzones = {cctk_nghostzones[0], cctk_nghostzones[1],
                                        cctk_nghostzones[2]};
   vect<int, dim> imin, imax;
-  GridDescBase(cctkGH).box_int<0, 0, 0>(nghostzones, imin, imax);
+  GridDescBase(cctkGH).box_int<1, 1, 1>(nghostzones, imin, imax);
   // Suffix 1: with ghost zones, suffix 0: without ghost zones
   const GF3D2layout layout1(cctkGH, indextype);
   const GF3D5layout layout0(imin, imax);
@@ -187,7 +187,7 @@ extern "C" void Z4c_Constraints(CCTK_ARGUMENTS) {
 #ifdef __CUDACC__
   const nvtxRangeId_t range = nvtxRangeStartA("Z4c_Constraints::constraints");
 #endif
-  grid.loop_int_device<0, 0, 0, vsize>(
+  grid.loop_int_device<1, 1, 1, vsize>(
       grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
         const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
         const GF3D2index index1(layout1, p.I);

diff --git a/Z4c/src/derivs.hxx b/Z4c/src/derivs.hxx
@@ -328,7 +328,7 @@ calc_derivs(const cGH *restrict const cctkGH, const GF3D2<const T> &gf1,
   const vec<CCTK_REAL, dim> dx([&](int a) { return CCTK_DELTA_SPACE(a); });
 
   const Loop::GridDescBaseDevice grid(cctkGH);
-  grid.loop_int_device<0, 0, 0, vsize>(
+  grid.loop_int_device<1, 1, 1, vsize>(
       grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
         const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
         const GF3D5index index0(layout0, p.I);
@@ -353,7 +353,7 @@ calc_derivs2(const cGH *restrict const cctkGH, const GF3D2<const T> &gf1,
   const vec<CCTK_REAL, dim> dx([&](int a) { return CCTK_DELTA_SPACE(a); });
 
   const Loop::GridDescBaseDevice grid(cctkGH);
-  grid.loop_int_device<0, 0, 0, vsize>(
+  grid.loop_int_device<1, 1, 1, vsize>(
       grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
         const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
         const int vavail = p.imax - p.i;
@@ -424,7 +424,7 @@ apply_upwind_diss(const cGH *restrict const cctkGH, const GF3D2<const T> &gf_,
   if (epsdiss == 0) {
 
     const Loop::GridDescBaseDevice grid(cctkGH);
-    grid.loop_int_device<0, 0, 0, vsize>(
+    grid.loop_int_device<1, 1, 1, vsize>(
         grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
           const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
           const vec<vreal, dim> betaG = gf_betaG_(mask, p.I);
@@ -437,7 +437,7 @@ apply_upwind_diss(const cGH *restrict const cctkGH, const GF3D2<const T> &gf_,
   } else {
 
     const Loop::GridDescBaseDevice grid(cctkGH);
-    grid.loop_int_device<0, 0, 0, vsize>(
+    grid.loop_int_device<1, 1, 1, vsize>(
         grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
           const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
           const vec<vreal, dim> betaG = gf_betaG_(mask, p.I);

diff --git a/Z4c/src/enforce.cxx b/Z4c/src/enforce.cxx
@@ -25,7 +25,7 @@ extern "C" void Z4c_Enforce(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTSX_Z4c_Enforce;
   DECLARE_CCTK_PARAMETERS;
 
-  const array<int, dim> indextype = {0, 0, 0};
+  const array<int, dim> indextype = {1, 1, 1};
   const GF3D2layout layout1(cctkGH, indextype);
 
   const GF3D2<CCTK_REAL> &gf_chi = chi;
@@ -49,7 +49,7 @@ extern "C" void Z4c_Enforce(CCTK_ARGUMENTS) {
 #ifdef __CUDACC__
   const nvtxRangeId_t range = nvtxRangeStartA("Z4c_Enforce::enforce");
 #endif
-  grid.loop_int_device<0, 0, 0, vsize>(
+  grid.loop_int_device<1, 1, 1, vsize>(
       grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
         const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
         const GF3D2index index1(layout1, p.I);

diff --git a/Z4c/src/field.hxx b/Z4c/src/field.hxx
@@ -36,7 +36,7 @@ template <typename T> class field {
 
 public:
   field(const GridDescBase &grid) : grid(grid) {
-    grid.box_all<0, 0, 0>(grid.nghostzones, imin, imax);
+    grid.box_all<1, 1, 1>(grid.nghostzones, imin, imax);
     ptrdiff_t str = 1;
     for (int d = 0; d < dim; ++d) {
       istr[d] = str;

diff --git a/Z4c/src/initial1.cxx b/Z4c/src/initial1.cxx
@@ -22,7 +22,7 @@ using namespace std;
 extern "C" void Z4c_Initial1(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTS_Z4c_Initial1;
 
-  const array<int, dim> indextype = {0, 0, 0};
+  const array<int, dim> indextype = {1, 1, 1};
   const GF3D2layout layout1(cctkGH, indextype);
 
   const smat<GF3D2<const CCTK_REAL>, 3> gf_g1{
@@ -80,7 +80,7 @@ extern "C" void Z4c_Initial1(CCTK_ARGUMENTS) {
 #ifdef __CUDACC__
   const nvtxRangeId_t range = nvtxRangeStartA("Z4c_Initial1::initial1");
 #endif
-  grid.loop_int_device<0, 0, 0, vsize>(
+  grid.loop_int_device<1, 1, 1, vsize>(
       grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
         const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
         const GF3D2index index1(layout1, p.I);

diff --git a/Z4c/src/initial2.cxx b/Z4c/src/initial2.cxx
@@ -25,7 +25,7 @@ extern "C" void Z4c_Initial2(CCTK_ARGUMENTS) {
       CCTK_DELTA_SPACE(2),
   };
 
-  const array<int, dim> indextype = {0, 0, 0};
+  const array<int, dim> indextype = {1, 1, 1};
   const GF3D2layout layout1(cctkGH, indextype);
 
   const smat<GF3D2<const CCTK_REAL>, 3> gf_gammat1{
@@ -50,7 +50,7 @@ extern "C" void Z4c_Initial2(CCTK_ARGUMENTS) {
 #ifdef __CUDACC__
   const nvtxRangeId_t range = nvtxRangeStartA("Z4c_Initial2::initial2");
 #endif
-  grid.loop_int_device<0, 0, 0, vsize>(
+  grid.loop_int_device<1, 1, 1, vsize>(
       grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
         const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
         const GF3D2index index1(layout1, p.I);

diff --git a/Z4c/src/rhs.cxx b/Z4c/src/rhs.cxx
@@ -42,11 +42,11 @@ extern "C" void Z4c_RHS(CCTK_ARGUMENTS) {
 
   //
 
-  const array<int, dim> indextype = {0, 0, 0};
+  const array<int, dim> indextype = {1, 1, 1};
   const array<int, dim> nghostzones = {cctk_nghostzones[0], cctk_nghostzones[1],
                                        cctk_nghostzones[2]};
   vect<int, dim> imin, imax;
-  GridDescBase(cctkGH).box_int<0, 0, 0>(nghostzones, imin, imax);
+  GridDescBase(cctkGH).box_int<1, 1, 1>(nghostzones, imin, imax);
   // Suffix 1: with ghost zones, suffix 0: without ghost zones
   const GF3D2layout layout1(cctkGH, indextype);
   const GF3D5layout layout0(imin, imax);
@@ -218,7 +218,7 @@ extern "C" void Z4c_RHS(CCTK_ARGUMENTS) {
   const nvtxRangeId_t range = nvtxRangeStartA("Z4c_RHS::rhs");
 #endif
   noinline([&]() __attribute__((__flatten__, __hot__)) {
-    grid.loop_int_device<0, 0, 0, vsize>(
+    grid.loop_int_device<1, 1, 1, vsize>(
         grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
           const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
           const GF3D2index index1(layout1, p.I);
@@ -259,7 +259,7 @@ extern "C" void Z4c_RHS(CCTK_ARGUMENTS) {
 #else
 
   noinline([&]() __attribute__((__flatten__, __hot__)) {
-    grid.loop_int_device<0, 0, 0, vsize>(
+    grid.loop_int_device<1, 1, 1, vsize>(
         grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
           const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
           const GF3D2index index1(layout1, p.I);
@@ -292,7 +292,7 @@ extern "C" void Z4c_RHS(CCTK_ARGUMENTS) {
   });
 
   noinline([&]() __attribute__((__flatten__, __hot__)) {
-    grid.loop_int_device<0, 0, 0, vsize>(
+    grid.loop_int_device<1, 1, 1, vsize>(
         grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
           const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
           const GF3D2index index1(layout1, p.I);