Initial commit.

6 files changed, 237 insertions, 0 deletions

diff --git a/configuration.ccl b/configuration.ccl
new file mode 100644
index 0000000..e0b0aba
--- /dev/null
+++ b/configuration.ccl
@@ -0,0 +1,3 @@
+# Configuration definition for thorn Trumpet
+
+REQUIRES GSL
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..535d955
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,15 @@
+# Interface definition for thorn Trumpet
+implements: Trumpet
+
+INHERITS: ADMBase StaticConformal grid CoordBase
+
+CCTK_INT FUNCTION GetDomainSpecification \
+  (CCTK_INT IN size, \
+   CCTK_REAL OUT ARRAY physical_min, \
+   CCTK_REAL OUT ARRAY physical_max, \
+   CCTK_REAL OUT ARRAY interior_min, \
+   CCTK_REAL OUT ARRAY interior_max, \
+   CCTK_REAL OUT ARRAY exterior_min, \
+   CCTK_REAL OUT ARRAY exterior_max, \
+   CCTK_REAL OUT ARRAY spacing)
+REQUIRES FUNCTION GetDomainSpecification
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..a009072
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,18 @@
+# Parameter definitions for thorn Trumpet
+
+SHARES: ADMBase
+
+EXTENDS KEYWORD initial_data
+{
+  "trumpet" :: "one maximal trumpet"
+}
+
+EXTENDS KEYWORD initial_lapse
+{
+ "trumpet" :: "one maximal trumpet"
+}
+
+EXTENDS KEYWORD initial_shift
+{
+ "trumpet" :: "one maximal trumpet"
+}
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..bd99c3f
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,20 @@
+# Schedule definitions for thorn TwoTrumpets
+#
+if (CCTK_Equals(initial_data, "trumpet")) {
+
+    SCHEDULE trumpet_data IN ADMBase_InitialData {
+        LANG: C
+    } "One maximal trumpet initial data"
+}
+
+if (CCTK_Equals(initial_lapse, "trumpet")) {
+    SCHEDULE trumpet_lapse IN ADMBase_InitialData {
+        LANG: C
+    } "One maximal trumpet initial lapse"
+}
+
+if (CCTK_Equals(initial_shift, "trumpet")) {
+    SCHEDULE trumpet_shift IN ADMBase_InitialData {
+        LANG: C
+    } "One maximal trumpet initial shift"
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..280a915
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,7 @@
+# Main make.code.defn file for thorn TwoPunctures
+
+# Source files in this directory
+SRCS = trumpet.c
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/trumpet.c b/src/trumpet.c
new file mode 100644
index 0000000..1598a0c
--- /dev/null
+++ b/src/trumpet.c
@@ -0,0 +1,174 @@
+#include <ctype.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gsl/gsl_spline.h>
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#define MAX(x,y) (x) > (y) ? (x) : (y)
+#define SQR(x) (x)*(x)
+#define C 3*sqrt(3)/4
+#define EPS 1E-10
+
+static inline double r_areal_to_isotropic(double r, double m)
+{
+    double par   = sqrt(4*r*r + 4*m*r + 3*m*m);
+    double term1 = (2*r + m + par)/4;
+    double term2 = (4 + 3*M_SQRT2)*(2*r - 3*m)/(8*r + 6*m + 3*M_SQRT2*par);
+    return term1*pow(term2, M_SQRT1_2);
+}
+
+/**
+ * get maximal isotropic r
+ */
+static CCTK_REAL get_max_r(CCTK_INT dim)
+{
+    CCTK_REAL phys_min[dim], phys_max[dim], ext_min[dim], ext_max[dim], int_min[dim], int_max[dim], step[dim];
+    CCTK_REAL max = 0;
+
+    GetDomainSpecification(dim, phys_min, phys_max, int_min, int_max, ext_min, ext_max, step);
+
+    for (int i = 0; i < dim; i++) {
+        max = MAX(max, abs(phys_min[i]));
+        max = MAX(max, abs(phys_max[i]));
+    }
+    return sqrt(dim)*2*max;
+}
+
+static void get_r_tables(double **pr_iso, double **pr_areal, CCTK_INT *size, CCTK_INT dim)
+{
+#define STEP 0.01
+    double *r_iso, *r_areal, max = get_max_r(dim);
+    int count = max*2/STEP, i = 0;
+
+    r_iso   = malloc(count*sizeof(*r_iso));
+    r_areal = malloc(count*sizeof(*r_areal));
+
+    for (i = 0;;i++) {
+        if (i >= count) {
+            count *= 2;
+            r_iso = realloc(r_iso, count*sizeof(*r_iso));
+            r_areal = realloc(r_areal, count*sizeof(*r_areal));
+        }
+
+        r_areal[i] = i*STEP + 1.5;
+        r_iso[i]   = r_areal_to_isotropic(r_areal[i], 1);
+        if (r_iso[i] > max)
+            break;
+    }
+    *pr_iso   = r_iso;
+    *pr_areal = r_areal;
+    *size     = i;
+}
+
+void trumpet_data(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    gsl_interp_accel *acc;
+    gsl_spline    *spline;
+
+    double *r_iso, *r_areal;
+    CCTK_INT size;
+
+    get_r_tables(&r_iso, &r_areal, &size, cctk_dim);
+
+    spline = gsl_spline_alloc(gsl_interp_cspline, size);
+    gsl_spline_init(spline, r_iso, r_areal, size);
+    acc = gsl_interp_accel_alloc();
+
+    memset(gxy, 0, sizeof(*gxy)*CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1));
+    memset(gxz, 0, sizeof(*gxy)*CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1));
+    memset(gyz, 0, sizeof(*gxy)*CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1));
+
+#pragma omp parallel for
+    for (int k = 0; k < cctk_lsh[2]; k++)
+        for (int j = 0; j < cctk_lsh[1]; j++)
+            for (int i = 0; i < cctk_lsh[0]; i++) {
+                int   index = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                CCTK_REAL r = sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]));
+                CCTK_REAL R = gsl_spline_eval(spline, r, acc);
+                if (r < EPS)
+                    r = EPS;
+                CCTK_REAL psi2 = R/r, psi4 = psi2*psi2;
+                CCTK_REAL k_fact = C/(r*r*r*r*R);
+
+                gxx[index] = gyy[index] = gzz[index] = psi4;
+
+                kxx[index] = -k_fact*(3*SQR(x[index]) - r*r);
+                kyy[index] = -k_fact*(3*SQR(y[index]) - r*r);
+                kzz[index] = -k_fact*(3*SQR(z[index]) - r*r);
+
+                kxy[index] = -k_fact*3*x[index]*y[index];
+                kxz[index] = -k_fact*3*x[index]*z[index];
+                kyz[index] = -k_fact*3*y[index]*z[index];
+            }
+}
+
+void trumpet_lapse(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    gsl_interp_accel *acc;
+    gsl_spline    *spline;
+
+    double *r_iso, *r_areal;
+    CCTK_INT size;
+
+    get_r_tables(&r_iso, &r_areal, &size, cctk_dim);
+
+    spline = gsl_spline_alloc(gsl_interp_cspline, size);
+    gsl_spline_init(spline, r_iso, r_areal, size);
+    acc = gsl_interp_accel_alloc();
+
+#pragma omp parallel for
+    for (int k = 0; k < cctk_lsh[2]; k++)
+        for (int j = 0; j < cctk_lsh[1]; j++)
+            for (int i = 0; i < cctk_lsh[0]; i++) {
+                int   index = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                CCTK_REAL r = sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS);
+                CCTK_REAL R = gsl_spline_eval(spline, r, acc);
+                if (r < EPS)
+                    r = EPS;
+
+                alp[index] = sqrt(1 - 2*1/R + C*C/(R*R*R*R));
+            }
+}
+
+void trumpet_shift(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    gsl_interp_accel *acc;
+    gsl_spline    *spline;
+
+    double *r_iso, *r_areal;
+    CCTK_INT size;
+
+    get_r_tables(&r_iso, &r_areal, &size, cctk_dim);
+
+    spline = gsl_spline_alloc(gsl_interp_cspline, size);
+    gsl_spline_init(spline, r_iso, r_areal, size);
+    acc = gsl_interp_accel_alloc();
+
+#pragma omp parallel for
+    for (int k = 0; k < cctk_lsh[2]; k++)
+        for (int j = 0; j < cctk_lsh[1]; j++)
+            for (int i = 0; i < cctk_lsh[0]; i++) {
+                int   index = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                CCTK_REAL r = sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS);
+                CCTK_REAL R = gsl_spline_eval(spline, r, acc);
+                if (r < EPS)
+                    r = EPS;
+
+                betax[index] = C*x[index]/(R*R*R);
+                betay[index] = C*y[index]/(R*R*R);
+                betaz[index] = C*z[index]/(R*R*R);
+            }
+}