From 11272e93ce2e5a39e7cd197e111fa3aa86f7f15e Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 13 Mar 2015 16:16:08 +0100 Subject: Split the basis sets definitions into a separate file. --- src/basis.c | 189 +++++++++++++++++++++++++++++++++++++++++ src/make.code.defn | 2 +- src/maximal_slicing_axi.c | 210 +--------------------------------------------- src/maximal_slicing_axi.h | 27 ++++++ 4 files changed, 220 insertions(+), 208 deletions(-) create mode 100644 src/basis.c create mode 100644 src/maximal_slicing_axi.h diff --git a/src/basis.c b/src/basis.c new file mode 100644 index 0000000..5ea043e --- /dev/null +++ b/src/basis.c @@ -0,0 +1,189 @@ + +#include + +#include "maximal_slicing_axi.h" + +static double sb_even_eval(double coord, int idx) +{ + double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord)); + + idx *= 2; // even only + + return sin((idx + 1) * val); +} + +static double sb_even_eval_diff1(double coord, int idx) +{ + double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord)); + + idx *= 2; // even only + + return - SCALE_FACTOR * (idx + 1) * SGN(coord) * cos((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord)); +} + +static double sb_even_eval_diff2(double coord, int idx) +{ + double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord)); + + idx *= 2; // even only + + return SCALE_FACTOR * (idx + 1) * SGN(coord) * (2 * fabs(coord) * cos((idx + 1) * val) - SCALE_FACTOR * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord)); +} + +static double sb_even_colloc_point(int order, int idx) +{ + double t; + + idx = order - idx - 1; + //order *= 2; + + //t = (idx + 2) * M_PI / (order + 4); + t = (idx + 2) * M_PI / (2 * order + 2); + return SCALE_FACTOR / tan(t); +} + +const BasisSet msa_sb_even_basis = { + .eval = sb_even_eval, + .eval_diff1 = sb_even_eval_diff1, + .eval_diff2 = sb_even_eval_diff2, + .colloc_point = sb_even_colloc_point, +}; + +static double tb_even_eval(double coord, int idx) +{ + double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord)); + + idx++; + idx *= 2; // even only + + return cos(idx * val) - 1.0; +} + +static double tb_even_eval_diff1(double coord, int idx) +{ + double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord)); + + idx++; + idx *= 2; // even only + + return SCALE_FACTOR * idx * SGN(coord) * sin(idx * val) / (SQR(SCALE_FACTOR) + SQR(coord)); +} + +static double tb_even_eval_diff2(double coord, int idx) +{ + double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord)); + + idx++; + idx *= 2; // even only + + return -SCALE_FACTOR * idx * SGN(coord) * (2 * fabs(coord) * sin(idx * val) + SCALE_FACTOR * idx * cos(idx * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord)); +} + +static double tb_even_colloc_point(int order, int idx) +{ + double t; + + idx = order - idx - 1; + //order *= 2; + + //t = (idx + 2) * M_PI / (order + 4); + t = (idx + 3) * M_PI / (2 * (order + 2)); + return SCALE_FACTOR / tan(t); +} + +const BasisSet msa_tb_even_basis = { + .eval = tb_even_eval, + .eval_diff1 = tb_even_eval_diff1, + .eval_diff2 = tb_even_eval_diff2, + .colloc_point = tb_even_colloc_point, +}; + +static double full_eval(double coord, int idx) +{ + double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord)); + int flag = idx & 1; + + idx /= 2; + + if (flag) return sin((2 * idx + 1) * 4 * val); + else return cos((2 * idx + 2) * 4 * val) - 1; +} + +static double full_eval_diff1(double coord, int idx) +{ + double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord)); + int flag = idx & 1; + + idx /= 2; + + if (flag) { + idx = 2 * idx + 1; + return -4 * idx * SCALE_FACTOR * cos(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord)); + } else { + idx = 2 * (idx + 1); + return 4 * idx * SCALE_FACTOR * sin(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord)); + } +} + +static double full_eval_diff2(double coord, int idx) +{ + double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord)); + int flag = idx & 1; + + idx /= 2; + + if (flag) { + idx = 2 * idx + 1; + return (16 * SQR(idx * SCALE_FACTOR) * cos(idx * 4 * val) - 4 * idx * SCALE_FACTOR * sin(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord)); + } else { + idx = 2 * (idx + 1); + return (-16 * SQR(idx * SCALE_FACTOR) * sin(idx * 4 * val) - 4 * idx * SCALE_FACTOR * cos(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord)); + } +} + +static double full_colloc_point(int order, int idx) +{ + double t; + + idx = order - idx - 1; + + t = (idx + 0.5) * M_PI / (2 * order); + + return SCALE_FACTOR / tan(t); + +} + +const BasisSet msa_full_basis = { + .eval = full_eval, + .eval_diff1 = full_eval_diff1, + .eval_diff2 = full_eval_diff2, + .colloc_point = full_colloc_point, +}; + +static double cheb_eval(double coord, int idx) +{ + return cos(2 * idx * acos(coord / SCALE_FACTOR)); +} + +static double cheb_eval_diff1(double coord, int idx) +{ + return 2 * idx * sin(2 * idx * acos(coord / SCALE_FACTOR)) / sqrt(SQR(SCALE_FACTOR) - SQR(coord)); +} + +static double cheb_eval_diff2(double coord, int idx) +{ + double t = acos(coord / SCALE_FACTOR); + return 2 * idx * (cos(2 * idx * t) * 2 * idx / (SQR(SCALE_FACTOR) - SQR(coord)) + sin(2 * idx * t) * coord / pow(SQR(SCALE_FACTOR) - SQR(coord), 1.5)); +} + +static double cheb_colloc_point(int order, int idx) +{ + return SCALE_FACTOR * cos((idx + 0.01) * M_PI / (4 * order + 4)); +} + +const BasisSet msa_cheb_basis = { + .eval = cheb_eval, + .eval_diff1 = cheb_eval_diff1, + .eval_diff2 = cheb_eval_diff2, + .colloc_point = cheb_colloc_point, +}; diff --git a/src/make.code.defn b/src/make.code.defn index 27625ce..a5091a4 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -1,7 +1,7 @@ # Main make.code.defn file for thorn MaximalSlicingAxi # Source files in this directory -SRCS = maximal_slicing_axi.c +SRCS = basis.c maximal_slicing_axi.c # Subdirectories containing source files SUBDIRS = diff --git a/src/maximal_slicing_axi.c b/src/maximal_slicing_axi.c index 7d3e6f4..ee39b3f 100644 --- a/src/maximal_slicing_axi.c +++ b/src/maximal_slicing_axi.c @@ -21,15 +21,12 @@ #include "util_Table.h" #include "Slicing.h" -double ms_scalarproduct_metric_avx(int offset_r, int offset_z, const double *coeffs, - const double *basis_r, const double *basis_z); +#include "maximal_slicing_axi.h" #define ACC_TEST 0 #define MAX(x, y) ((x) > (y) ? (x) : (y)) #define MIN(x, y) ((x) > (y) ? (y) : (x)) -#define SQR(x) ((x) * (x)) -#define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0) #define ARRAY_ELEMS(arr) (sizeof(arr) / sizeof(*arr)) /* @@ -45,215 +42,14 @@ int64_t gettime(void) return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; } -/* a set of basis functions */ -typedef struct BasisSet { - /* evaluate the idx-th basis function at the specified point*/ - long double (*eval) (long double coord, int idx); - /* evaluate the first derivative of the idx-th basis function at the specified point*/ - long double (*eval_diff1)(long double coord, int idx); - /* evaluate the second derivative of the idx-th basis function at the specified point*/ - long double (*eval_diff2)(long double coord, int idx); - /** - * Get the idx-th collocation point for the specified order. - * idx runs from 0 to order - 1 (inclusive) - */ - long double (*colloc_point)(int order, int idx); -} BasisSet; - /* * The basis of even (n = 2 * idx) SB functions (Boyd 2000, Ch 17.9) * SB(x, n) = sin((n + 1) arccot(|x| / L)) * They are symmetric wrt origin and decay as 1/x in infinity. */ -static CCTK_REAL scale_factor; - -#define SCALE_FACTOR scale_factor - -static long double sb_even_eval(long double coord, int idx) -{ - long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord)); - - idx *= 2; // even only - - return sinl((idx + 1) * val); -} - -static long double sb_even_eval_diff1(long double coord, int idx) -{ - long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord)); - - idx *= 2; // even only - - return - SCALE_FACTOR * (idx + 1) * SGN(coord) * cosl((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord)); -} - -static long double sb_even_eval_diff2(long double coord, int idx) -{ - long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord)); - - idx *= 2; // even only - - return SCALE_FACTOR * (idx + 1) * SGN(coord) * (2 * fabsl(coord) * cosl((idx + 1) * val) - SCALE_FACTOR * (idx + 1) * sinl((idx + 1) * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord)); -} - -static long double sb_even_colloc_point(int order, int idx) -{ - long double t; - - idx = order - idx - 1; - //order *= 2; - - //t = (idx + 2) * M_PI / (order + 4); - t = (idx + 2) * M_PI / (2 * order + 2); - return SCALE_FACTOR / tanl(t); -} +double scale_factor; -static const BasisSet sb_even_basis = { - .eval = sb_even_eval, - .eval_diff1 = sb_even_eval_diff1, - .eval_diff2 = sb_even_eval_diff2, - .colloc_point = sb_even_colloc_point, -}; - -static long double tb_even_eval(long double coord, int idx) -{ - long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord)); - - idx++; - idx *= 2; // even only - - return cosl(idx * val) - 1.0; -} - -static long double tb_even_eval_diff1(long double coord, int idx) -{ - long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord)); - - idx++; - idx *= 2; // even only - - return SCALE_FACTOR * idx * SGN(coord) * sinl(idx * val) / (SQR(SCALE_FACTOR) + SQR(coord)); -} - -static long double tb_even_eval_diff2(long double coord, int idx) -{ - long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord)); - - idx++; - idx *= 2; // even only - - return -SCALE_FACTOR * idx * SGN(coord) * (2 * fabsl(coord) * sinl(idx * val) + SCALE_FACTOR * idx * cosl(idx * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord)); -} - -static long double tb_even_colloc_point(int order, int idx) -{ - long double t; - - idx = order - idx - 1; - //order *= 2; - - //t = (idx + 2) * M_PI / (order + 4); - t = (idx + 3) * M_PI / (2 * (order + 2)); - return SCALE_FACTOR / tanl(t); -} - -static const BasisSet tb_even_basis = { - .eval = tb_even_eval, - .eval_diff1 = tb_even_eval_diff1, - .eval_diff2 = tb_even_eval_diff2, - .colloc_point = tb_even_colloc_point, -}; - -static long double full_eval(long double coord, int idx) -{ - long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord)); - int flag = idx & 1; - - idx /= 2; - - if (flag) return sinl((2 * idx + 1) * 4 * val); - else return cosl((2 * idx + 2) * 4 * val) - 1; -} - -static long double full_eval_diff1(long double coord, int idx) -{ - long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord)); - int flag = idx & 1; - - idx /= 2; - - if (flag) { - idx = 2 * idx + 1; - return -4 * idx * SCALE_FACTOR * cosl(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord)); - } else { - idx = 2 * (idx + 1); - return 4 * idx * SCALE_FACTOR * sinl(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord)); - } -} - -static long double full_eval_diff2(long double coord, int idx) -{ - long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord)); - int flag = idx & 1; - - idx /= 2; - - if (flag) { - idx = 2 * idx + 1; - return (16 * SQR(idx * SCALE_FACTOR) * cosl(idx * 4 * val) - 4 * idx * SCALE_FACTOR * sinl(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord)); - } else { - idx = 2 * (idx + 1); - return (-16 * SQR(idx * SCALE_FACTOR) * sinl(idx * 4 * val) - 4 * idx * SCALE_FACTOR * cosl(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord)); - } -} - -static long double full_colloc_point(int order, int idx) -{ - long double t; - - idx = order - idx - 1; - - t = (idx + 0.5) * M_PI / (2 * order); - - return SCALE_FACTOR / tanl(t); - -} - -static const BasisSet full_basis = { - .eval = full_eval, - .eval_diff1 = full_eval_diff1, - .eval_diff2 = full_eval_diff2, - .colloc_point = full_colloc_point, -}; - -static long double cheb_eval(long double coord, int idx) -{ - return cosl(2 * idx * acosl(coord / SCALE_FACTOR)); -} - -static long double cheb_eval_diff1(long double coord, int idx) -{ - return 2 * idx * sinl(2 * idx * acosl(coord / SCALE_FACTOR)) / sqrtl(SQR(SCALE_FACTOR) - SQR(coord)); -} - -static long double cheb_eval_diff2(long double coord, int idx) -{ - long double t = acosl(coord / SCALE_FACTOR); - return 2 * idx * (cosl(2 * idx * t) * 2 * idx / (SQR(SCALE_FACTOR) - SQR(coord)) + sinl(2 * idx * t) * coord / pow(SQR(SCALE_FACTOR) - SQR(coord), 1.5)); -} - -static long double cheb_colloc_point(int order, int idx) -{ - return SCALE_FACTOR * cosl((idx + 0.01) * M_PI / (4 * order + 4)); -} - -static const BasisSet cheb_basis = { - .eval = cheb_eval, - .eval_diff1 = cheb_eval_diff1, - .eval_diff2 = cheb_eval_diff2, - .colloc_point = cheb_colloc_point, -}; /* indices (in our code, not cactus structs) of the grid functions which we'll need to * interpolate on the pseudospectral grid */ @@ -1163,7 +959,7 @@ static MaximalSlicingContext *init_ms(cGH *cctkGH, ms->gh = cctkGH; - ms->basis = &sb_even_basis; + ms->basis = &msa_sb_even_basis; //ms->basis = &full_basis; //ms->basis = &cheb_basis; diff --git a/src/maximal_slicing_axi.h b/src/maximal_slicing_axi.h new file mode 100644 index 0000000..b15420c --- /dev/null +++ b/src/maximal_slicing_axi.h @@ -0,0 +1,27 @@ + +#define SQR(x) ((x) * (x)) +#define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0) + +#define SCALE_FACTOR scale_factor + +/* a set of basis functions */ +typedef struct BasisSet { + /* evaluate the idx-th basis function at the specified point*/ + double (*eval) (double coord, int idx); + /* evaluate the first derivative of the idx-th basis function at the specified point*/ + double (*eval_diff1)(double coord, int idx); + /* evaluate the second derivative of the idx-th basis function at the specified point*/ + double (*eval_diff2)(double coord, int idx); + /** + * Get the idx-th collocation point for the specified order. + * idx runs from 0 to order - 1 (inclusive) + */ + double (*colloc_point)(int order, int idx); +} BasisSet; + +extern const BasisSet msa_cheb_basis; +extern const BasisSet msa_full_basis; +extern const BasisSet msa_tb_even_basis; +extern const BasisSet msa_sb_even_basis; + +extern double scale_factor; -- cgit v1.2.3