summaryrefslogtreecommitdiff
path: root/src/maximal_slicing_axi.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/maximal_slicing_axi.c')
-rw-r--r--src/maximal_slicing_axi.c210
1 files changed, 3 insertions, 207 deletions
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;
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-17 07:58:22 +0000