From d1beabbdbeedab43f2ae372acdab26a2f32d1443 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Mon, 9 Jan 2012 15:41:26 +0100 Subject: Add boosting along the x-axis. --- src/mdefs.h | 41 ++++++++++++++ src/trumpet.c | 177 +++++++++++++++++++++++++++++++++++++++++++++++++--------- 2 files changed, 193 insertions(+), 25 deletions(-) create mode 100644 src/mdefs.h (limited to 'src') diff --git a/src/mdefs.h b/src/mdefs.h new file mode 100644 index 0000000..dc5d0c1 --- /dev/null +++ b/src/mdefs.h @@ -0,0 +1,41 @@ +/************************************************************************* + + Mathematica source file + + Copyright 1986 through 1999 by Wolfram Research Inc. + + +*************************************************************************/ + +/* C language definitions for use with Mathematica output */ + + +#define Power(x, y) (pow((double)(x), (double)(y))) +#define Sqrt(x) (sqrt((double)(x))) + +#define Abs(x) (fabs((double)(x))) + +#define Exp(x) (exp((double)(x))) +#define Log(x) (log((double)(x))) + +#define Sin(x) (sin((double)(x))) +#define Cos(x) (cos((double)(x))) +#define Tan(x) (tan((double)(x))) + +#define ArcSin(x) (asin((double)(x))) +#define ArcCos(x) (acos((double)(x))) +#define ArcTan(x) (atan((double)(x))) + +#define Sinh(x) (sinh((double)(x))) +#define Cosh(x) (cosh((double)(x))) +#define Tanh(x) (tanh((double)(x))) + + +#define E 2.71828182845904523536029 +#define Pi 3.14159265358979323846264 +#define Degree 0.01745329251994329576924 + + +/** Could add definitions for Random(), SeedRandom(), etc. **/ + + diff --git a/src/trumpet.c b/src/trumpet.c index 7b65822..0933c04 100644 --- a/src/trumpet.c +++ b/src/trumpet.c @@ -8,6 +8,8 @@ #include "cctk_Arguments.h" #include "cctk_Parameters.h" +#include "mdefs.h" + #define MAX(x,y) (x) > (y) ? (x) : (y) #define SQR(x) ((x)*(x)) @@ -23,11 +25,13 @@ /* * isotropic/coordinate radius */ -#define ISO_R(index) (sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS)) +#define ISO_R(x, y, z, index, gamma) (sqrt(SQR(gamma*x[index]) + SQR(y[index]) + SQR(z[index]) + EPS)) + +#define TRUMPET_ALPHA(R) (sqrt(1 - 2*MASS/R + SQR(TRUMPET_CONST)/SQR(SQR(R)))) static inline double r_areal_to_isotropic(double r, double m) { - double par = sqrt(4*r*r + 4*m*r + 3*m*m); + double par = sqrt(4*SQR(r) + 4*m*r + 3*SQR(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); @@ -76,11 +80,101 @@ static void get_r_tables(double **pr_iso, double **pr_areal, CCTK_INT *size, CCT *size = i; } +/* those equations come from trumpet.nb */ +static long double sqrt_factor(long double R, long double r, long double x, long double alpha, long double beta) +{ + long double R2 = SQR(R); + return sqrtl((R2*(R + alpha*beta*r) - beta*TRUMPET_CONST*x)*(R2*(R - alpha*beta*r) - beta*TRUMPET_CONST*x)); +} + +static inline CCTK_REAL K_11(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, + CCTK_REAL R, CCTK_REAL r, + CCTK_REAL beta, CCTK_REAL M, + CCTK_REAL alpha) +{ + long double b2 = SQR(beta); + long double r2 = SQR(r); + long double R2 = SQR(R); + + long double n1 = TRUMPET_CONST*(1.0 - 3*SQR(x/r))*(beta*x*TRUMPET_CONST/R - R2)/r2; + long double n2 = b2*beta*M*(-Power(TRUMPET_CONST,2)/(R2*R2) + Power(alpha,2))*x; + long double n3 = b2*TRUMPET_CONST*M*(1.0 + 2*SQR(x/r))/R; + long double n4 = beta*x*(SQR(R/r)*(-2*M + (-1 + alpha)*alpha*R)); + + long double den = (b2 - 1)*sqrt_factor(R, r, x, alpha, beta); + + return (n1 + n2 + n3 + n4)/den; +} + +static inline CCTK_REAL K_22(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, + CCTK_REAL R, CCTK_REAL r, + CCTK_REAL beta, CCTK_REAL M, + CCTK_REAL alpha) +{ + long double r2 = SQR(r); + long double R2 = SQR(R); + + long double n1 = TRUMPET_CONST*(R2 - beta*TRUMPET_CONST*(x/R))*(1.0 - 3*SQR(y/r))/r2; + long double n2 = alpha*(alpha - 1)*beta*SQR(R/r)*R*x; + + long double sq = sqrt_factor(R, r, x, alpha, beta); + + return (n1 + n2)/sq; +} + +static inline CCTK_REAL K_33(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, + CCTK_REAL R, CCTK_REAL r, + CCTK_REAL beta, CCTK_REAL M, + CCTK_REAL alpha) +{ + return K_22(x, z, y, R, r, beta, M, alpha); +} + +static inline CCTK_REAL K_12(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, + CCTK_REAL R, CCTK_REAL r, + CCTK_REAL beta, CCTK_REAL M, + CCTK_REAL alpha) +{ + long double b2 = SQR(beta); + long double r2 = SQR(r); + long double R2 = SQR(R); + + long double n1 = TRUMPET_CONST*(b2*M/R + 3*SQR(R/r))*x - 3*beta*SQR(TRUMPET_CONST)*SQR(x/r)/R; + long double n2 = beta*SQR(R/r)*(-M + (-1 + alpha)*alpha*R); + long double den = -sqrt_factor(R, r, x, alpha, beta)*sqrtl(1 - b2); // really minus? + + return y*(n1/r2 + n2)/den; +} + +static inline CCTK_REAL K_13(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, + CCTK_REAL R, CCTK_REAL r, + CCTK_REAL beta, CCTK_REAL M, + CCTK_REAL alpha) +{ + return K_12(x, z, y, R, r, beta, M, alpha); +} + +static inline CCTK_REAL K_23(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, + CCTK_REAL R, CCTK_REAL r, + CCTK_REAL beta, CCTK_REAL M, + CCTK_REAL alpha) +{ + long double r2 = SQR(r); + long double R2 = SQR(R); + + long double num = -3*TRUMPET_CONST*(R2 - beta*TRUMPET_CONST*x/R)*(y/r)*(z/r); + long double den = r2*sqrt_factor(R, r, x, alpha, beta); + + return num/den; +} + void trumpet_data(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; + double gamma = 1.0/sqrt(1.0 - SQR(boost_velocity)); + gsl_interp_accel *acc; gsl_spline *spline; @@ -93,8 +187,6 @@ void trumpet_data(CCTK_ARGUMENTS) 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 @@ -102,21 +194,32 @@ 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 = ISO_R(index); + CCTK_REAL xx = gamma*x[index], yy = y[index], zz = z[index]; + CCTK_REAL r = ISO_R(x, y, z, index, gamma); CCTK_REAL R = gsl_spline_eval(spline, r, acc); - CCTK_REAL psi2 = R/r, psi4 = psi2*psi2; - 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]) - 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]; - kyz[index] = -k_fact*3*y[index]*z[index]; + CCTK_REAL alpha = TRUMPET_ALPHA(R); + + long double r2 = SQR(r); + long double R2 = SQR(R); + long double b2 = SQR(boost_velocity); + + kxx[index] = K_11(xx, yy, zz, R, r, boost_velocity, MASS, alpha); + kyy[index] = K_22(xx, yy, zz, R, r, boost_velocity, MASS, alpha); + kzz[index] = K_33(xx, yy, zz, R, r, boost_velocity, MASS, alpha); + kxy[index] = K_12(xx, yy, zz, R, r, boost_velocity, MASS, alpha); + kxz[index] = K_13(xx, yy, zz, R, r, boost_velocity, MASS, alpha); + kyz[index] = K_23(xx, yy, zz, R, r, boost_velocity, MASS, alpha); + + gxx[index] = (-SQR(R/r) + b2*(-SQR(TRUMPET_CONST/R2) + Power(alpha,2)) + 2*boost_velocity*TRUMPET_CONST*xx/(R*r2)) / (b2 - 1); + gyy[index] = SQR(R/r); + gzz[index] = SQR(R/r); + gxy[index] = -((boost_velocity*TRUMPET_CONST*yy)/(Sqrt(1 - b2)*r2*R)); + gxz[index] = -((boost_velocity*TRUMPET_CONST*zz)/(Sqrt(1 - b2)*r2*R)); } + free(r_iso); + free(r_areal); + gsl_interp_accel_free(acc); + gsl_spline_free(spline); } void trumpet_lapse(CCTK_ARGUMENTS) @@ -124,6 +227,8 @@ void trumpet_lapse(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; + double gamma = 1.0/sqrt(1.0 - SQR(boost_velocity)); + gsl_interp_accel *acc; gsl_spline *spline; @@ -141,11 +246,21 @@ void trumpet_lapse(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]) + EPS); + CCTK_REAL xx = gamma*x[index], yy = y[index], zz = z[index]; + CCTK_REAL r = ISO_R(x, y, z, index, gamma); CCTK_REAL R = gsl_spline_eval(spline, r, acc); + CCTK_REAL alpha = TRUMPET_ALPHA(R); - alp[index] = sqrt(1 - 2*MASS/R + TRUMPET_CONST*TRUMPET_CONST/(R*R*R*R)); + long double r2 = SQR(r); + long double R2 = SQR(R); + long double b2 = SQR(boost_velocity); + + alp[index] = alpha*R2*R*sqrt((b2 - 1) / (Power(alpha,2)*b2*r2*R2*R2 - Power(R2*R - boost_velocity*TRUMPET_CONST*xx,2))); } + free(r_iso); + free(r_areal); + gsl_interp_accel_free(acc); + gsl_spline_free(spline); } void trumpet_shift(CCTK_ARGUMENTS) @@ -153,6 +268,8 @@ void trumpet_shift(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; + double gamma = 1.0/sqrt(1.0 - SQR(boost_velocity)); + gsl_interp_accel *acc; gsl_spline *spline; @@ -170,14 +287,24 @@ void trumpet_shift(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]) + EPS); + CCTK_REAL xx = gamma*x[index], yy = y[index], zz = z[index]; + CCTK_REAL r = ISO_R(x, y, z, index, gamma); CCTK_REAL R = gsl_spline_eval(spline, r, acc); + CCTK_REAL alpha = TRUMPET_ALPHA(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); - } -} + long double r2 = SQR(r); + long double R2 = SQR(R); + long double b2 = SQR(boost_velocity); + long double a2 = SQR(alpha); + + betax[index] = -((a2*boost_velocity*r2*R2*R2 + TRUMPET_CONST*Power(R,3)*xx + b2*TRUMPET_CONST*Power(R,3)*xx - boost_velocity*(Power(R,6) + Power(TRUMPET_CONST,2)*Power(xx,2)))/ + (a2*b2*r2*R2*R2 - Power(Power(R,3) - boost_velocity*TRUMPET_CONST*xx,2))); + betay[index] = (Sqrt(1 - b2)*TRUMPET_CONST*(-Power(R,3) + boost_velocity*TRUMPET_CONST*xx)*yy)/(a2*b2*r2*R2*R2 - Power(Power(R,3) - boost_velocity*TRUMPET_CONST*xx,2)); + betaz[index] = (Sqrt(1 - b2)*TRUMPET_CONST*(-Power(R,3) + boost_velocity*TRUMPET_CONST*xx)*zz)/(a2*b2*r2*R2*R2 - Power(Power(R,3) - boost_velocity*TRUMPET_CONST*xx,2)); } + free(r_iso); + free(r_areal); + gsl_interp_accel_free(acc); + gsl_spline_free(spline); } -- cgit v1.2.3