From dbba307e91e7e44ff9e1a2b01992e07398898b4c Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 23 Apr 2020 16:53:50 +0200 Subject: Initial commit. --- configuration.ccl | 2 + interface.ccl | 29 +++++++++++ param.ccl | 12 +++++ schedule.ccl | 39 +++++++++++++++ src/make.code.defn | 7 +++ src/nullsurf.c | 142 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 231 insertions(+) create mode 100644 configuration.ccl create mode 100644 interface.ccl create mode 100644 param.ccl create mode 100644 schedule.ccl create mode 100644 src/make.code.defn create mode 100644 src/nullsurf.c diff --git a/configuration.ccl b/configuration.ccl new file mode 100644 index 0000000..12661fa --- /dev/null +++ b/configuration.ccl @@ -0,0 +1,2 @@ +# Configuration definition for thorn NullSurf + diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..612f623 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,29 @@ +# Interface definition for thorn NullSurf +implements: NullSurf + +INHERITS: ADMBase grid CoordBase MethodOfLines ML_BSSN Driver + +CCTK_INT FUNCTION MoLRegisterEvolved(CCTK_INT IN EvolvedIndex, CCTK_INT IN RHSIndex) + +REQUIRES FUNCTION MoLRegisterEvolved + +CCTK_POINTER FUNCTION \ + VarDataPtrI \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN map, \ + CCTK_INT IN reflevel, \ + CCTK_INT IN component, \ + CCTK_INT IN timelevel, \ + CCTK_INT IN varindex) +USES FUNCTION VarDataPtrI + +public: +REAL null_surface TYPE=GF TIMELEVELS=3 tags='tensortypealias="Scalar" tensorweight=0' +{ + F +} + +REAL null_surface_rhs TYPE=GF TIMELEVELS=3 tags='tensortypealias="Scalar" tensorweight=0 Prolongation="None"' +{ + F_rhs +} diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..0210523 --- /dev/null +++ b/param.ccl @@ -0,0 +1,12 @@ +# Parameter definitions for thorn NullSurf + +shares: MethodOfLines + +USES CCTK_INT MoL_Num_Evolved_Vars + +restricted: +CCTK_INT mol_evolved "Number of evolved variables used by this thorn" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Evolved_Vars STEERABLE=RECOVER +{ + 1:1 :: "Number of evolved variables used by this thorn" +} 1 + diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..d6c38b3 --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,39 @@ +# Schedule definitions for thorn NullSurf +# +SCHEDULE ns_initial IN ADMBase_InitialData { + LANG: C +} "NullSurf initial data" + +SCHEDULE ns_eval_rhs IN MoL_CalcRHS { + LANG: C + READS: ML_BSSN::alpha(Interior) + READS: ML_BSSN::phi(Interior) + READS: ML_BSSN::gt11(Interior) + READS: ML_BSSN::gt12(Interior) + READS: ML_BSSN::gt13(Interior) + READS: ML_BSSN::gt22(Interior) + READS: ML_BSSN::gt23(Interior) + READS: ML_BSSN::gt33(Interior) + READS: ML_BSSN::beta1(Interior) + READS: ML_BSSN::beta2(Interior) + READS: ML_BSSN::beta3(Interior) + WRITES: NullSurf::F(Interior) +} "NullSurf eval RHS" + +SCHEDULE ns_mol_register in MoL_Register { + LANG: C +} "NullSurf register MoL variables" + +SCHEDULE ns_register_symmetries in SymmetryRegister { + LANG: C +} "NullSurf register symmetry properties" + +schedule ns_select_bc in MoL_PostStep +{ + LANG: C + OPTIONS: level + SYNC: null_surface +} "select boundary conditions" + +STORAGE: null_surface[3] +STORAGE: null_surface_rhs[3] diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..8fa0281 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,7 @@ +# Main make.code.defn file for thorn MaximalSlicingAxi + +# Source files in this directory +SRCS = nullsurf.c + +# Subdirectories containing source files +SUBDIRS = 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 +#include +#include + +#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"); +} -- cgit v1.2.3