CurvADM: add thorn

CurvADM/README

@@ -0,0 +1,9 @@
+Cactus Code Thorn CurvADM
+Author(s)    : Liwei Ji <jiliwei.phys@gmail.com>
+Maintainer(s): Liwei Ji <jiliwei.phys@gmail.com>
+Licence      : LGPL
+--------------------------------------------------------------------------
+
+1. Purpose
+
+Provide storage for the ADM variables

CurvADM/configuration.ccl

@@ -0,0 +1,3 @@
+# Configuration definition for thorn CurvADM
+
+REQUIRES Loop

CurvADM/interface.ccl

@@ -0,0 +1,29 @@
+# Interface definition for thorn CurvADM
+
+IMPLEMENTS: CurvADM
+
+USES INCLUDE HEADER: loop_device.hxx
+
+
+PUBLIC:
+
+CCTK_REAL metric TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { gxx gxy gxz gyy gyz gzz } "ADM 3-metric g_ij"
+
+CCTK_REAL excurv TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { kxx kxy kxz kyy kyz kzz } "ADM extrinsic curvature K_ij"
+
+CCTK_REAL lapse TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { alp } "ADM lapse function alpha"
+
+CCTK_REAL dtlapse TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dtalp } "Time derivative of ADM lapse function"
+
+CCTK_REAL shift TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { betax betay betaz } "ADM shift vector beta^i"
+
+CCTK_REAL dtshift TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dtbetax dtbetay dtbetaz } "Time derivative of ADM shift vector"
+
+
+# Variables necessary to calculate second time derivatives of the four-metric
+
+CCTK_REAL dtexcurv TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dtkxx dtkxy dtkxz dtkyy dtkyz dtkzz } "Time derivative of ADM extrinsic curvature K_ij"
+
+CCTK_REAL dt2lapse TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dt2alp } "Second time derivative of ADM lapse function"
+
+CCTK_REAL dt2shift TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dt2betax dt2betay dt2betaz } "Second time derivative of ADM shift vector"

CurvADM/param.ccl

@@ -0,0 +1,47 @@
+# Parameter definitions for thorn CurvADM
+
+RESTRICTED:
+
+KEYWORD initial_data "Initial metric and extrinsic curvature datasets"
+{
+  "none" :: "Do not set up initial conditions"
+  "Cartesian Minkowski" :: "Minkowski geometry in Cartesian coordinates"
+  "linear wave" :: "linear wave"
+} "Cartesian Minkowski"
+
+KEYWORD initial_lapse "Initial lapse value"
+{
+  "none" :: "Do not set up initial conditions"
+  "one" :: "Uniform lapse"
+} "one"
+
+KEYWORD initial_shift "Initial shift value"
+{
+  "none" :: "Do not set up initial conditions"
+  "zero" :: "Shift is zero"
+} "zero"
+
+KEYWORD initial_dtlapse "Initial dtlapse value"
+{
+  "none" :: "Do not set up initial conditions"
+  "zero" :: "Dtlapse is zero"
+} "zero"
+
+KEYWORD initial_dtshift "Initial dtshift value"
+{
+  "none" :: "Do not set up initial conditions"
+  "zero" :: "Dtshift is zero"
+} "zero"
+
+
+PRIVATE:
+
+CCTK_REAL linear_wave_amplitude "Linear wave amplitude"
+{
+  0.0:* :: ""
+} 1.0e-8
+
+CCTK_REAL linear_wave_wavelength "Linear wave wavelength"
+{
+  0.0:* :: ""
+} 1.0

CurvADM/schedule.ccl

@@ -0,0 +1,105 @@
+# Schedule definitions for thorn CurvADM
+
+if (CCTK_IsThornActive("ODESolvers")) {
+
+  SCHEDULE GROUP CurvADM_InitialData IN ODESolvers_Initial
+  {
+  } "Schedule group for calculating ADM initial data"
+
+  SCHEDULE GROUP CurvADM_InitialGauge IN ODESolvers_Initial AFTER CurvADM_InitialData
+  {
+  } "Schedule group for the ADM initial gauge condition"
+
+  SCHEDULE GROUP CurvADM_PostInitial IN ODESolvers_Initial AFTER (CurvADM_InitialData CurvADM_InitialGauge)
+  {
+  } "Schedule group for modifying the ADM initial data, such as e.g. adding noise"
+
+  SCHEDULE GROUP CurvADM_SetADMVars IN ODESolvers_PostStep
+  {
+  } "Set ADM variables in this group"
+
+  SCHEDULE GROUP CurvADM_SetADMRHS IN ODESolvers_PostStep
+  {
+  } "Set ADM RHS variables in this group"
+
+} else {
+
+  SCHEDULE GROUP CurvADM_SetADMVars AT post_recover_variables
+  {
+  } "Set ADM variables in this group"
+
+  SCHEDULE GROUP CurvADM_SetADMRHS AT post_recover_variables
+  {
+  } "Set ADM RHS variables in this group"
+
+  SCHEDULE GROUP CurvADM_SetADMVars AT postregrid
+  {
+  } "Set ADM variables in this group"
+
+  SCHEDULE GROUP CurvADM_SetADMRHS AT postregrid
+  {
+  } "Set ADM RHS variables in this group"
+
+  SCHEDULE GROUP CurvADM_SetADMVars AT postrestrict
+  {
+  } "Set ADM variables in this group"
+
+  SCHEDULE GROUP CurvADM_SetADMRHS AT postrestrict
+  {
+  } "Set ADM RHS variables in this group"
+
+  SCHEDULE GROUP CurvADM_SetADMVars AT poststep
+  {
+  } "Set ADM variables in this group"
+
+  SCHEDULE GROUP CurvADM_SetADMRHS AT poststep
+  {
+  } "Set ADM RHS variables in this group"
+
+}
+
+if (CCTK_EQUALS(initial_data, "Cartesian Minkowski")) {
+  SCHEDULE CurvADM_initial_data IN CurvADM_InitialData
+  {
+    LANG: C
+    WRITES: metric(everywhere) curv(everywhere)
+  } "Set up Cartesian Minkowski initial data"
+} else if (CCTK_EQUALS(initial_data, "linear wave")) {
+  SCHEDULE CurvADM_linear_wave IN CurvADM_InitialData
+  {
+    LANG: C
+    WRITES: metric(everywhere) curv(everywhere)
+  } "Set up linear wave initial data"
+}
+
+if (CCTK_EQUALS(initial_lapse, "one")) {
+  SCHEDULE CurvADM_initial_lapse IN CurvADM_InitialGauge
+  {
+    LANG: C
+    WRITES: lapse(everywhere)
+  } "Set lapse to one"
+}
+
+if (CCTK_EQUALS(initial_dtlapse, "zero")) {
+  SCHEDULE CurvADM_initial_dtlapse IN CurvADM_InitialGauge
+  {
+    LANG: C
+    WRITES: dtlapse(everywhere)
+  } "Set dtlapse to zero"
+}
+
+if (CCTK_EQUALS(initial_shift, "zero")) {
+  SCHEDULE CurvADM_initial_shift IN CurvADM_InitialGauge
+  {
+    LANG: C
+    WRITES: shift(everywhere)
+  } "Set shift to zero"
+}
+
+if (CCTK_EQUALS(initial_dtshift, "zero")) {
+  SCHEDULE CurvADM_initial_dtshift IN CurvADM_InitialGauge
+  {
+    LANG: C
+    WRITES: dtshift(everywhere)
+  } "Set dtshift to zero"
+}

CurvADM/src/adm.cxx

@@ -0,0 +1,80 @@
+#include <loop_device.hxx>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <array>
+
+namespace CurvADM {
+using namespace std;
+using namespace Loop;
+
+extern "C" void CurvADM_initial_data(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_CurvADM_initial_data;
+  DECLARE_CCTK_PARAMETERS;
+
+  grid.loop_all_device<1, 1, 1>(grid.nghostzones,
+                                [=] CCTK_DEVICE(const PointDesc &p)
+                                    CCTK_ATTRIBUTE_ALWAYS_INLINE {
+                                      gxx(p.I) = 1;
+                                      gxy(p.I) = 0;
+                                      gxz(p.I) = 0;
+                                      gyy(p.I) = 1;
+                                      gyz(p.I) = 0;
+                                      gzz(p.I) = 1;
+
+                                      kxx(p.I) = 0;
+                                      kxy(p.I) = 0;
+                                      kxz(p.I) = 0;
+                                      kyy(p.I) = 0;
+                                      kyz(p.I) = 0;
+                                      kzz(p.I) = 0;
+                                    });
+}
+
+extern "C" void CurvADM_initial_lapse(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_CurvADM_initial_lapse;
+  DECLARE_CCTK_PARAMETERS;
+
+  grid.loop_all_device<1, 1, 1>(
+      grid.nghostzones, [=] CCTK_DEVICE(const PointDesc &p)
+                            CCTK_ATTRIBUTE_ALWAYS_INLINE { alp(p.I) = 1; });
+}
+
+extern "C" void CurvADM_initial_dtlapse(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_CurvADM_initial_dtlapse;
+  DECLARE_CCTK_PARAMETERS;
+
+  grid.loop_all_device<1, 1, 1>(
+      grid.nghostzones, [=] CCTK_DEVICE(const PointDesc &p)
+                            CCTK_ATTRIBUTE_ALWAYS_INLINE { dtalp(p.I) = 0; });
+}
+
+extern "C" void CurvADM_initial_shift(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_CurvADM_initial_shift;
+  DECLARE_CCTK_PARAMETERS;
+
+  grid.loop_all_device<1, 1, 1>(grid.nghostzones,
+                                [=] CCTK_DEVICE(const PointDesc &p)
+                                    CCTK_ATTRIBUTE_ALWAYS_INLINE {
+                                      betax(p.I) = 0;
+                                      betay(p.I) = 0;
+                                      betaz(p.I) = 0;
+                                    });
+}
+
+extern "C" void CurvADM_initial_dtshift(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_CurvADM_initial_dtshift;
+  DECLARE_CCTK_PARAMETERS;
+
+  grid.loop_all_device<1, 1, 1>(grid.nghostzones,
+                                [=] CCTK_DEVICE(const PointDesc &p)
+                                    CCTK_ATTRIBUTE_ALWAYS_INLINE {
+                                      dtbetax(p.I) = 0;
+                                      dtbetay(p.I) = 0;
+                                      dtbetaz(p.I) = 0;
+                                    });
+}
+
+} // namespace CurvADM

CurvADM/src/linear_wave.cxx

@@ -0,0 +1,47 @@
+#include <loop.hxx>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <cmath>
+
+namespace CurvADM {
+using namespace Loop;
+using namespace std;
+
+extern "C" void CurvADM_linear_wave(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_CurvADM_linear_wave;
+  DECLARE_CCTK_PARAMETERS;
+
+  const CCTK_REAL t = cctk_time;
+
+  // See arXiv:1111.2177 [gr-qc], (74-75)
+
+  const auto b = [&](const PointDesc &p) {
+    return linear_wave_amplitude *
+           sin(2 * CCTK_REAL(M_PI) * (p.x - t) / linear_wave_wavelength);
+  };
+
+  const auto bt = [&](const PointDesc &p) {
+    return -2 * CCTK_REAL(M_PI) * linear_wave_amplitude /
+           linear_wave_wavelength *
+           cos(2 * CCTK_REAL(M_PI) * (p.x - t) / linear_wave_wavelength);
+  };
+
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gxx(p.I) = 1; });
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gxy(p.I) = 0; });
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gxz(p.I) = 0; });
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gyy(p.I) = 1 + b(p); });
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gyz(p.I) = 0; });
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gzz(p.I) = 1 - b(p); });
+
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kxx(p.I) = 0; });
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kxy(p.I) = 0; });
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kxz(p.I) = 0; });
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kyy(p.I) = bt(p) / 2; });
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kyz(p.I) = 0; });
+  loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kzz(p.I) = -bt(p) / 2; });
+}
+
+} // namespace CurvADM

CurvADM/src/make.code.defn

@@ -0,0 +1,7 @@
+# Main make.code.defn file for thorn CurvADM
+
+# Source files in this directory
+SRCS = adm.cxx linear_wave.cxx
+
+# Subdirectories containing source files
+SUBDIRS =