From cf9324106c743433640e213bb60bf9f6ebdd4b15 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Tue, 11 Oct 2011 12:04:46 +0200 Subject: Initial commit. --- configuration.ccl | 3 + interface.ccl | 15 +++++ param.ccl | 18 ++++++ schedule.ccl | 20 ++++++ src/make.code.defn | 7 +++ src/trumpet.c | 174 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 237 insertions(+) create mode 100644 configuration.ccl create mode 100644 interface.ccl create mode 100644 param.ccl create mode 100644 schedule.ccl create mode 100644 src/make.code.defn create mode 100644 src/trumpet.c 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 +#include +#include +#include +#include +#include +#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); + } +} -- cgit v1.2.3