diff options
| -rw-r--r-- | param.ccl | 8 | ||||
| -rw-r--r-- | src/TwoPunctures.c | 34 |
2 files changed, 38 insertions, 4 deletions
@@ -30,6 +30,14 @@ BOOLEAN keep_u_around "Keep the variable u around after solving" +KEYWORD grid_setup_method "How to fill the 3D grid from the spectral grid" +{ + "Taylor expansion" :: "use a Taylor expansion about the nearest collocation point (fast, but might be inaccurate)" + "interpolation" :: "interpolate using all spectral coefficients (slow)" +} "Taylor expansion" + + + INT npoints_A "Number of coefficients in the compactified radial direction" { 4:* :: "" diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c index 519aabb..005b98f 100644 --- a/src/TwoPunctures.c +++ b/src/TwoPunctures.c @@ -29,6 +29,9 @@ TwoPunctures (CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; + enum GRID_SETUP_METHOD { GSM_Taylor_expansion, GSM_interpolation }; + enum GRID_SETUP_METHOD gsm; + int nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi; int i, j, k, ntotal = n1 * n2 * n3 * nvar; @@ -47,6 +50,19 @@ TwoPunctures (CCTK_ARGUMENTS) F_of_v (nvar, n1, n2, n3, v, F, u); } + if (CCTK_EQUALS(grid_setup_method, "Taylor expansion")) + { + gsm = GSM_Taylor_expansion; + } + else if (CCTK_EQUALS(grid_setup_method, "interpolation")) + { + gsm = GSM_interpolation; + } + else + { + CCTK_WARN (0, "internal error"); + } + CCTK_INFO ("Interpolating result"); if (CCTK_EQUALS(metric_type, "static conformal")) { if (CCTK_EQUALS(conformal_storage, "factor")) { @@ -74,10 +90,20 @@ TwoPunctures (CCTK_ARGUMENTS) const double r_minus = sqrt(pow2(x[ind] + par_b) + pow2(y[ind]) + pow2(z[ind])); - const double U = PunctTaylorExpandAtArbitPosition - (0, nvar, n1, n2, n3, v, x[ind], y[ind], z[ind]); -/* const double U = PunctIntPolAtArbitPosition */ -/* (0, nvar, n1, n2, n3, v, x[ind], y[ind], z[ind]); */ + double U; + switch (gsm) + { + case GSM_Taylor_expansion: + U = PunctTaylorExpandAtArbitPosition + (0, nvar, n1, n2, n3, v, x[ind], y[ind], z[ind]); + break; + case GSM_interpolation: + U = PunctIntPolAtArbitPosition + (0, nvar, n1, n2, n3, v, x[ind], y[ind], z[ind]); + break; + default: + assert (0); + } const double psi1 = 1 + 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus + U; |