diff options
Diffstat (limited to 'src/nullsurf.c')
| -rw-r--r-- | src/nullsurf.c | 142 |
1 files changed, 142 insertions, 0 deletions
diff --git a/src/nullsurf.c b/src/nullsurf.c new file mode 100644 index 0000000..2b88437 --- /dev/null +++ b/src/nullsurf.c @@ -0,0 +1,142 @@ + +#include <math.h> +#include <stddef.h> +#include <stdint.h> + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "Symmetry.h" + +#define SQR(x) ((x) * (x)) + +void ns_mol_register(CCTK_ARGUMENTS) +{ + MoLRegisterEvolved(CCTK_VarIndex("NullSurf::F"), CCTK_VarIndex("NullSurf::F_rhs")); +} + +void ns_initial(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + for (int j = 0; j < cctkGH->cctk_lsh[2]; j++) { + for (int i = 0; i < cctkGH->cctk_lsh[0]; i++) { + const int idx = CCTK_GFINDEX3D(cctkGH, i, 0, j); + const double r2 = SQR(x[idx]) + SQR(z[idx]); + + F[idx] = r2; + } + } +} + +#define FD4(arr, stride, h_inv) \ + ((-1.0 * arr[2 * stride] + 8.0 * arr[1 * stride] - 8.0 * arr[-1 * stride] + 1.0 * arr[-2 * stride]) * h_inv / 12.0) +#define FD24(arr, stride, h_inv) \ + ((-1.0 * arr[2 * stride] + 16.0 * arr[1 * stride] - 30.0 * arr[0] + 16.0 * arr[-1 * stride] - 1.0 * arr[-2 * stride]) * SQR(h_inv) / 12.0) + +static double diff4_dx(const double *arr, ptrdiff_t stride, double h_inv) +{ + return FD4(arr, 1, h_inv); +} +static double diff4_dz(const double *arr, ptrdiff_t stride, double h_inv) +{ + return FD4(arr, stride, h_inv); +} + +void ns_eval_rhs(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + const int ghosts_x = cctkGH->cctk_nghostzones[0]; + const int ghosts_z = cctkGH->cctk_nghostzones[2]; + + const ptrdiff_t stride_z = CCTK_GFINDEX3D(cctkGH, 0, 0, 1); + + const double dx = x[1] - x[0]; + const double dz = z[stride_z] - z[0]; + const double dx_inv = 1.0 / dx; + const double dz_inv = 1.0 / dz; + + for (int idx_z = ghosts_z; idx_z < cctkGH->cctk_lsh[2] - ghosts_z; idx_z++) { + for (int idx_x = ghosts_x; idx_x < cctkGH->cctk_lsh[0] - ghosts_x; idx_x++) { + const int idx = CCTK_GFINDEX3D(cctkGH, idx_x, 0, idx_z); + + const double dF[3] = {diff4_dx(&F[idx], 1, dx_inv), 0.0, diff4_dz(&F[idx], stride_z, dz_inv)}; + + const double alpha_val = alpha[idx]; + const double beta_val[3] = { beta1[idx], 0.0, beta3[idx] }; + const double phi_val = phi[idx]; + + const double gtxx = gt11[idx]; + const double gtxy = gt12[idx]; + const double gtxz = gt13[idx]; + const double gtyy = gt22[idx]; + const double gtyz = gt23[idx]; + const double gtzz = gt33[idx]; + + const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz); + + double gtu[3][3], g4u[4][4]; + double shift_term; + + // \tilde{γ}^{ij} + gtu[0][0] = (gtyy * gtzz - SQR(gtyz)) / det; + gtu[1][1] = (gtxx * gtzz - SQR(gtxz)) / det; + gtu[2][2] = (gtxx * gtyy - SQR(gtxy)) / det; + gtu[0][1] = -(gtxy * gtzz - gtyz * gtxz) / det; + gtu[0][2] = (gtxy * gtyz - gtyy * gtxz) / det; + gtu[1][2] = -(gtxx * gtyz - gtxy * gtxz) / det; + gtu[1][0] = gtu[0][1]; + gtu[2][0] = gtu[0][2]; + gtu[2][1] = gtu[1][2]; + + g4u[0][0] = - 1.0 / SQR(alpha_val); + for (int i = 0; i < 3; i++) { + const double val = beta_val[i] / SQR(alpha_val); + g4u[1 + i][0] = val; + g4u[0][1 + i] = val; + } + + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { + g4u[1 + i][1 + j] = SQR(phi_val) * gtu[i][j] - beta_val[i] * beta_val[j] / SQR(alpha_val); + } + + shift_term = 0.0; + for (int i = 0; i < 3; i++) + shift_term += g4u[0][1 + i] * dF[i]; + + double dF2 = 0.0; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + dF2 += g4u[1 + i][1 + j] * dF[i] * dF[j]; + + F_rhs[idx] = (-shift_term + sqrt(SQR(shift_term) - g4u[0][0] * dF2)) / g4u[0][0]; + } + } +} + +void ns_register_symmetries(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + int sym[3] = { 1, 1, 1 }; + int ret; + + ret = SetCartSymVN(cctkGH, sym, "NullSurf::F"); + if (ret != 0) + CCTK_WARN(0, "Error registering symmetries"); +} + +void ns_select_bc(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + int ret = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, "NullSurf::F", "none"); + if (ret != 0) + CCTK_WARN(0, "Error registering boundaries"); +} |