From ae9311a976086674bde25ce9dc67c939e2cf2cf1 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Mon, 9 Jan 2012 15:41:06 +0100 Subject: trumpet: macrofication --- src/trumpet.c | 47 ++++++++++++++++++++++++++++------------------- 1 file changed, 28 insertions(+), 19 deletions(-) diff --git a/src/trumpet.c b/src/trumpet.c index 1598a0c..7b65822 100644 --- a/src/trumpet.c +++ b/src/trumpet.c @@ -9,9 +9,21 @@ #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 +#define SQR(x) ((x)*(x)) + +#define MASS 1 +/* + * the constant C from e.g. 0908.1063v4 + */ +#define TRUMPET_CONST (3*sqrt(3)*SQR(MASS)/4) +/* + * small number to avoid r=0 singularities + */ +#define EPS 1E-08 +/* + * isotropic/coordinate radius + */ +#define ISO_R(index) (sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS)) static inline double r_areal_to_isotropic(double r, double m) { @@ -40,7 +52,7 @@ static CCTK_REAL get_max_r(CCTK_INT dim) static void get_r_tables(double **pr_iso, double **pr_areal, CCTK_INT *size, CCTK_INT dim) { -#define STEP 0.01 +#define STEP 0.0001 double *r_iso, *r_areal, max = get_max_r(dim); int count = max*2/STEP, i = 0; @@ -90,18 +102,16 @@ void trumpet_data(CCTK_ARGUMENTS) 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 = ISO_R(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); + CCTK_REAL k_fact = TRUMPET_CONST/(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); + kxx[index] = -k_fact*(3*SQR(x[index]) - SQR(r)); + kyy[index] = -k_fact*(3*SQR(y[index]) - SQR(r)); + kzz[index] = -k_fact*(3*SQR(z[index]) - SQR(r)); kxy[index] = -k_fact*3*x[index]*y[index]; kxz[index] = -k_fact*3*x[index]*z[index]; @@ -133,10 +143,8 @@ void trumpet_lapse(CCTK_ARGUMENTS) 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)); + alp[index] = sqrt(1 - 2*MASS/R + TRUMPET_CONST*TRUMPET_CONST/(R*R*R*R)); } } @@ -164,11 +172,12 @@ void trumpet_shift(CCTK_ARGUMENTS) 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); + betax[index] = TRUMPET_CONST*x[index]/(R*R*R); + betay[index] = TRUMPET_CONST*y[index]/(R*R*R); + betaz[index] = TRUMPET_CONST*z[index]/(R*R*R); + } +} + } } -- cgit v1.2.3