Add thorn ADMconstraints (#50)

* ADMconstraints: add admconstraints.cxx

* ADMconstraints: add ADM_vars.wl

* ADMconstraints: add ADM_rhs.wl

* ADMconstraints: add type_in_rhs.nb

* ADMconstraints: nv type_in_rhs.nb outside wl

* ADMconstraints: update type_in_rhs.nb

* ADMconstraints: update ADM_vars.wl

* ADMconstraints: update type_in_rhs.nb

* ADMconstraints: update type_in_rhs.nb

* ADMconstraints: update ADM_rhs.wl

* ADMconstraints: add check_rhs.ipynb

* ADMconstraints: update check_rhs.ipynb

* ADMconstraints: add ADM_set_constraint.wl

* ADMconstraints: add ADM_set_constraint.hxx

* ADMconstraints: add interface.ccl

* ADMconstraints: rename MtC to MC

* ADMconstraints: add param.ccl

* ADMconstraints: update admconstraints.cxx

* ADMconstraints: update admconstraints.cxx

* ADMconstraints: remove prefix ADM for variables

* ADMconstraints: update admconstraints.cxx

* ADMconstraints: update admconstraints.cxx

* ADMconstraints: add make.code.defn

* ADMconstraints: add configuration.ccl

* ADMconstraints: add schedule.ccl

* ADMconstraints: make it compile

* ADMconstraints: format ADM_rhs.wl

* ADMconstraints: treat DexK separatedly to remove warnings

* ADMconstraints: treat dexK separatedly to remove warnings

* ADMconstraints: add lapse mask

* ADMconstraints: add extra SYNCs

* ADMconstraints: fix comments in schedule

* ADMconstraints: add par lapse_mask_outer_radius

* ADMconstraints: mask out region which is larger than a radius

* ADMconstraints: fix type of default value for lapse_mask_outer_radius

* ADMconstraints: use or instead of and when apply mask

* ADMconstraints: sync HC and MC after calculation

ADMconstraints/configuration.ccl

@@ -0,0 +1,3 @@
+# Configuration definitions for thorn ADMconstraints
+
+REQUIRES Arith Loop Derivs

ADMconstraints/interface.ccl

@@ -0,0 +1,20 @@
+# Interface definition for thorn ADMconstraints
+
+IMPLEMENTS: ADMconstraints
+
+INHERITS: ADMBaseX TmunuBaseX
+
+USES INCLUDE HEADER: defs.hxx
+USES INCLUDE HEADER: derivs.hxx
+USES INCLUDE HEADER: div.hxx
+USES INCLUDE HEADER: dual.hxx
+USES INCLUDE HEADER: loop_device.hxx
+USES INCLUDE HEADER: mat.hxx
+USES INCLUDE HEADER: simd.hxx
+USES INCLUDE HEADER: sum.hxx
+USES INCLUDE HEADER: vec.hxx
+USES INCLUDE HEADER: vect.hxx
+
+CCTK_REAL HC TYPE=gf TAGS='checkpoint="no"' "H"
+
+CCTK_REAL MC TYPE=gf TAGS='parities={-1 +1 +1   +1 -1 +1   +1 +1 -1} checkpoint="no"' { MCx MCy MCz } "M"

ADMconstraints/param.ccl

@@ -0,0 +1,27 @@
+# Parameter definitions for thorn ADMconstraints
+
+BOOLEAN calc_constraints "Calculate constraints" STEERABLE=recover
+{
+} yes
+
+CCTK_INT deriv_order "Order of spatial finite differencing" STEERABLE=never
+{
+  2 :: "Second order finite difference"
+  4 :: "Fourth order finite difference"
+  6 :: "Sixth order finite difference"
+  8 :: "Eighth order finite difference"
+} 4
+
+BOOLEAN use_lapse_mask "Exclude the region when lapse is small" STEERABLE=recover
+{
+} yes
+
+CCTK_REAL lapse_mask_cutoff "Exclude the region when lapse is smaller than the cutoff" STEERABLE=always
+{
+  *:* :: ""
+} 0.3
+
+CCTK_REAL lapse_mask_outer_radius "Exclude the region when the distance to the origin is larger than the radius" STEERABLE=always
+{
+  *:* :: ""
+} 30.0

ADMconstraints/schedule.ccl

@@ -0,0 +1,52 @@
+# Schedule definitions for thorn ADMconstraints
+
+STORAGE: HC
+STORAGE: MC
+
+
+
+
+SCHEDULE GROUP ADMconstraints_AnalysisGroup AT analysis
+{
+} "Analyse ADM variables"
+
+
+if (calc_constraints) {
+  SCHEDULE ADMconstraints_Sync IN ADMconstraints_AnalysisGroup
+  {
+    LANG: C
+    OPTIONS: global
+    SYNC: ADMBaseX::metric
+    SYNC: ADMBaseX::curv
+    SYNC: ADMBaseX::lapse
+    SYNC: ADMBaseX::shift
+  } "Synchronize"
+
+  SCHEDULE ADMconstraints_CalcConstraints IN ADMconstraints_AnalysisGroup AFTER ADMconstraints_Sync
+  {
+    LANG: C
+    READS: ADMBaseX::metric(everywhere)
+    READS: ADMBaseX::curv(everywhere)
+    READS: ADMBaseX::lapse(everywhere)
+    READS: ADMBaseX::shift(everywhere)
+    READS: TmunuBaseX::eTtt(interior)
+    READS: TmunuBaseX::eTti(interior)
+    READS: TmunuBaseX::eTij(interior)
+    WRITES: HC(interior)
+    WRITES: MC(interior)
+    SYNC: HC MC
+  } "Calculate ADM Constraints"
+
+  if (use_lapse_mask) {
+    SCHEDULE ADMconstraints_LapseMask IN ADMconstraints_AnalysisGroup AFTER ADMconstraints_CalcConstraints
+    {
+      LANG: C
+      READS: ADMBaseX::lapse(interior)
+      READS: HC(interior)
+      READS: MC(interior)
+      WRITES: HC(interior)
+      WRITES: MC(interior)
+      SYNC: HC MC
+    } "Apply lapse mask"
+  }
+}

ADMconstraints/src/admconstraints.cxx

@@ -0,0 +1,145 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#ifdef __CUDACC__
+// Disable CCTK_DEBUG since the debug information takes too much
+// parameter space to launch the kernels
+#ifdef CCTK_DEBUG
+#undef CCTK_DEBUG
+#endif
+#endif
+
+#include <derivs.hxx>
+#include <loop_device.hxx>
+#include <simd.hxx>
+
+#ifdef __CUDACC__
+#include <nvtx3/nvToolsExt.h>
+#endif
+
+#include <cmath>
+
+namespace ADMconstraints {
+using namespace Arith;
+using namespace Loop;
+
+template <typename T>
+CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T Power(T x, int y) {
+  return (y == 2) ? Arith::pow2(x) : Arith::pown(x, y);
+}
+
+extern "C" void ADMconstraints_Sync(CCTK_ARGUMENTS) {
+  // do nothing
+}
+
+extern "C" void ADMconstraints_CalcConstraints(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_ADMconstraints_CalcConstraints;
+  DECLARE_CCTK_PARAMETERS;
+
+  for (int d = 0; d < 3; ++d)
+    if (cctk_nghostzones[d] < deriv_order / 2)
+      CCTK_VERROR("Need at least %d ghost zones", deriv_order / 2);
+
+  const vect<CCTK_REAL, dim> dx{
+      CCTK_DELTA_SPACE(0),
+      CCTK_DELTA_SPACE(1),
+      CCTK_DELTA_SPACE(2),
+  };
+
+  const array<int, dim> indextype = {0, 0, 0};
+  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);
+  // suffix 2: with ghost zones, suffix 5: without ghost zones
+  const GF3D2layout layout2(cctkGH, indextype);
+  const GF3D5layout layout5(imin, imax);
+
+  // Input grid functions
+  const smat<GF3D2<const CCTK_REAL>, 3> gf_gam{gxx, gxy, gxz, gyy, gyz, gzz};
+  const smat<GF3D2<const CCTK_REAL>, 3> gf_exK{kxx, kxy, kxz, kyy, kyz, kzz};
+  const GF3D2<const CCTK_REAL> &gf_alpha = alp;
+  const vec<GF3D2<const CCTK_REAL>, 3> gf_beta{betax, betay, betaz};
+
+  // More input grid functions
+  const GF3D2<const CCTK_REAL> &gf_eTtt = eTtt;
+  const vec<GF3D2<const CCTK_REAL>, 3> gf_eTt{eTtx, eTty, eTtz};
+  const smat<GF3D2<const CCTK_REAL>, 3> gf_eT{eTxx, eTxy, eTxz,
+                                              eTyy, eTyz, eTzz};
+
+  // Output grid functions
+  const GF3D2<CCTK_REAL> &gf_HC = HC;
+  const vec<GF3D2<CCTK_REAL>, 3> gf_MC{MCx, MCy, MCz};
+
+  // Define derivs lambdas
+  const auto calccopy = [&](const auto &gf, const auto &gf0) {
+    Derivs::calc_copy<0, 0, 0>(gf, layout5, grid, gf0);
+  };
+  const auto calcderivs = [&](const auto &gf, const auto &dgf,
+                              const auto &gf0) {
+    Derivs::calc_derivs<0, 0, 0>(gf, dgf, layout5, grid, gf0, dx, deriv_order);
+  };
+  const auto calcderivs2 = [&](const auto &gf, const auto &dgf,
+                               const auto &ddgf, const auto &gf0) {
+    Derivs::calc_derivs2<0, 0, 0>(gf, dgf, ddgf, layout5, grid, gf0, dx,
+                                  deriv_order);
+  };
+
+  // Tile variables for derivatives and so on
+  const int ntmps = 88;
+  GF3D5vector<CCTK_REAL> tmps(layout5, ntmps);
+  int itmp = 0;
+
+  const auto make_gf = [&]() { return GF3D5<CCTK_REAL>(tmps(itmp++)); };
+  const auto make_vec = [&](const auto &f) {
+    return vec<result_of_t<decltype(f)()>, 3>([&](int) { return f(); });
+  };
+  const auto make_mat = [&](const auto &f) {
+    return smat<result_of_t<decltype(f)()>, 3>([&](int, int) { return f(); });
+  };
+  const auto make_vec_gf = [&]() { return make_vec(make_gf); };
+  const auto make_mat_gf = [&]() { return make_mat(make_gf); };
+  const auto make_mat_vec_gf = [&]() { return make_mat(make_vec_gf); };
+  const auto make_mat_mat_gf = [&]() { return make_mat(make_mat_gf); };
+
+  const smat<GF3D5<CCTK_REAL>, 3> tl_gam(make_mat_gf());
+  const smat<vec<GF3D5<CCTK_REAL>, 3>, 3> tl_dgam(make_mat_vec_gf());
+  const smat<smat<GF3D5<CCTK_REAL>, 3>, 3> tl_ddgam(make_mat_mat_gf());
+  calcderivs2(tl_gam, tl_dgam, tl_ddgam, gf_gam);
+
+  const smat<GF3D5<CCTK_REAL>, 3> tl_exK(make_mat_gf());
+  const smat<vec<GF3D5<CCTK_REAL>, 3>, 3> tl_dexK(make_mat_vec_gf());
+  calcderivs(tl_exK, tl_dexK, gf_exK);
+
+  const GF3D5<CCTK_REAL> tl_alpha(make_gf());
+  calccopy(tl_alpha, gf_alpha);
+
+  const vec<GF3D5<CCTK_REAL>, 3> tl_beta(make_vec_gf());
+  calccopy(tl_beta, gf_beta);
+
+  if (itmp != ntmps)
+    CCTK_VERROR("Wrong number of temporary variables: ntmps=%d itmp=%d", ntmps,
+                itmp);
+  itmp = -1;
+
+  // simd types
+  typedef simd<CCTK_REAL> vreal;
+  typedef simdl<CCTK_REAL> vbool;
+  constexpr size_t vsize = tuple_size_v<vreal>;
+
+  // parameters
+  const CCTK_REAL cpi = M_PI;
+
+#ifdef __CUDACC__
+  const nvtxRangeId_t range = nvtxRangeStartA("ADMconstraints::constraints");
+#endif
+
+#include "../wolfram/ADM_set_constraint.hxx"
+
+#ifdef __CUDACC__
+  nvtxRangeEnd(range);
+#endif
+}
+
+} // namespace ADMconstraints

ADMconstraints/src/lapsemask.cxx

@@ -0,0 +1,29 @@
+#include <loop_device.hxx>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+
+namespace ADMconstraints {
+using namespace Arith;
+using namespace Loop;
+
+extern "C" void ADMconstraints_LapseMask(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_ADMconstraints_LapseMask;
+  DECLARE_CCTK_PARAMETERS;
+
+  const CCTK_REAL local_cutoff = lapse_mask_cutoff;
+  const CCTK_REAL local_outer_radius = lapse_mask_outer_radius;
+
+  grid.loop_int_device<0, 0, 0>(
+      grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
+        const CCTK_REAL rad = sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
+        if (alp(p.I) < local_cutoff || rad > local_outer_radius) {
+          HC(p.I) = 0.0;
+          MCx(p.I) = 0.0;
+          MCy(p.I) = 0.0;
+          MCz(p.I) = 0.0;
+        }
+      });
+}
+
+} // namespace ADMconstraints

ADMconstraints/src/make.code.defn

@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn Z4c
+
+# Source files in this directory
+SRCS = \
+	admconstraints.cxx \
+  lapsemask.cxx
+
+# Subdirectories containing source files
+SUBDIRS =