From b497f900c2274ec120752fd0bcc455d3f4c385ac Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Tue, 7 Aug 2018 15:17:41 +0200 Subject: Initial commit. --- configuration.ccl | 2 + interface.ccl | 15 ++ param.ccl | 10 + schedule.ccl | 24 +++ src/make.code.defn | 7 + src/maximal_slicing_axi_mg.c | 497 +++++++++++++++++++++++++++++++++++++++++++ src/register.c | 18 ++ 7 files changed, 573 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/maximal_slicing_axi_mg.c create mode 100644 src/register.c diff --git a/configuration.ccl b/configuration.ccl new file mode 100644 index 0000000..f124283 --- /dev/null +++ b/configuration.ccl @@ -0,0 +1,2 @@ +# Configuration definition for thorn MaximalSlicingAxiMG + diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..43e512c --- /dev/null +++ b/interface.ccl @@ -0,0 +1,15 @@ +# Interface definition for thorn MaximalSlicingAxiMG +implements: MaximalSlicingAxiMG + +INHERITS: ADMBase grid CoordBase MethodOfLines ML_BSSN + +CCTK_INT FUNCTION MoLRegisterConstrained(CCTK_INT IN idx) +CCTK_INT FUNCTION MoLRegisterSaveAndRestore(CCTK_INT IN idx) +CCTK_INT FUNCTION MoLRegisterSaveAndRestoreGroup(CCTK_INT IN idx) + +CCTK_INT FUNCTION MoLNumIntegratorSubsteps () +USES FUNCTION MoLNumIntegratorSubsteps + +REQUIRES FUNCTION MoLRegisterConstrained +REQUIRES FUNCTION MoLRegisterSaveAndRestore +REQUIRES FUNCTION MoLRegisterSaveAndRestoreGroup diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..e6b39e2 --- /dev/null +++ b/param.ccl @@ -0,0 +1,10 @@ +# Parameter definitions for thorn MaximalSlicingAxiMG + +SHARES: ADMBase + +EXTENDS KEYWORD lapse_evolution_method +{ + "maximal_axi_mg" :: "Maximal slicing for an axisymmetric spacetime" +} + +RESTRICTED: diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..246aa63 --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,24 @@ +# Schedule definitions for thorn MaximalSlicingAxiMG +# +if (CCTK_Equals(lapse_evolution_method, "maximal_axi_mg")) { + + SCHEDULE maximal_slicing_axi_mg IN ML_BSSN_evolCalcGroup BEFORE ML_BSSN_RHS { + LANG: C + } "Maximal slicing in axisymmetry" + + #SCHEDULE maximal_slicing_axi IN MoL_PostStepModify { + # LANG: C + #} "Maximal slicing in axisymmetry" + + #SCHEDULE maximal_slicing_axi IN MoL_PseudoEvolution { + # LANG: C + #} "Maximal slicing in axisymmetry" + + SCHEDULE maximal_slicing_axi_mg_register AT Startup { + LANG: C + } "" + + SCHEDULE maximal_slicing_axi_mg_register_mol IN MoL_Register { + LANG: C + } "" +} diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..4bc0d9c --- /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 = maximal_slicing_axi_mg.c register.c + +# Subdirectories containing source files +SUBDIRS = diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c new file mode 100644 index 0000000..668a702 --- /dev/null +++ b/src/maximal_slicing_axi_mg.c @@ -0,0 +1,497 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Timers.h" +#include "util_Table.h" + +#define SQR(x) ((x) * (x)) +#define ARRAY_ELEMS(x) (sizeof(x) / sizeof(*x)) + +/* precomputed values for a given refined grid */ +typedef struct CoordPatch { + double origin[3]; + int delta[3]; + int size[3]; + int level; + + MG2DContext *solver; + + ptrdiff_t y_idx; + + /* number of x/z grid points by which the elliptic solver domain is offset + * from the cactus grid */ + ptrdiff_t offset_left[2]; + ptrdiff_t offset_right[2]; + + /* boundary interpolator parameters */ + double *interp_coords[3]; + + double *boundary_interp_vals; + size_t nb_boundary_interp_vals; + + int interp_values_codes[1]; + int interp_var_indices[1]; + int interp_operation_codes[1]; + int interp_operation_indices[1]; + + int coord_system; + int interp_operator; + int interp_params; + + /** + * pointers into boundary_interp_vals + * the first row are the boundary points for the solver + * [0] is the upper-z patch + * rows run from x=0 to x[lsh[0] - 1] + * [1] is the upper-x patch + * rows run from z=0 to z= + */ + double *boundary_val[2]; + ptrdiff_t boundary_val_stride[2]; + +} CoordPatch; + +typedef struct MSMGContext { + cGH *gh; + + CoordPatch *patches; + int nb_patches; +} MSMGContext; + +static int ctz(int a) +{ + int ret = 0; + + if (!a) + return INT_MAX; + + while (!(a & 1)) { + a >>= 1; + ret++; + } + + return ret; +} + +static void coord_patch_free(CoordPatch *cp) +{ + mg2d_solver_free(&cp->solver); +} + +static CoordPatch *get_coord_patch(MSMGContext *ms) +{ + cGH *gh = ms->gh; + + const double *a_x = CCTK_VarDataPtr(gh, 0, "grid::x"); + const double *a_y = CCTK_VarDataPtr(gh, 0, "grid::y"); + const double *a_z = CCTK_VarDataPtr(gh, 0, "grid::z"); + + CoordPatch *cp; + size_t domain_size[2]; + int integrator_substeps = MoLNumIntegratorSubsteps(); + int i, ret; + + for (int i = 0; i < ms->nb_patches; i++) { + cp = &ms->patches[i]; + + if (cp->origin[0] == ms->gh->cctk_origin_space[0] && + cp->origin[1] == ms->gh->cctk_origin_space[1] && + cp->origin[2] == ms->gh->cctk_origin_space[2] && + cp->size[0] == ms->gh->cctk_lsh[0] && + cp->size[1] == ms->gh->cctk_lsh[1] && + cp->size[2] == ms->gh->cctk_lsh[2] && + cp->delta[0] == ms->gh->cctk_levfac[0] && + cp->delta[1] == ms->gh->cctk_levfac[1] && + cp->delta[2] == ms->gh->cctk_levfac[2]) + return cp; + } + + /* create a new patch */ + ms->patches = realloc(ms->patches, sizeof(*ms->patches) * (ms->nb_patches + 1)); + cp = &ms->patches[ms->nb_patches]; + + memset(cp, 0, sizeof(*cp)); + + memcpy(cp->origin, ms->gh->cctk_origin_space, sizeof(cp->origin)); + memcpy(cp->size, ms->gh->cctk_lsh, sizeof(cp->size)); + memcpy(cp->delta, ms->gh->cctk_levfac, sizeof(cp->delta)); + cp->level = ctz(ms->gh->cctk_levfac[0]); + + for (i = 0; i < cp->size[1]; i++) + if (fabs(a_y[CCTK_GFINDEX3D(gh, 0, i, 0)]) < 1e-8) { + cp->y_idx = i; + break; + } + if (i == cp->size[1]) + CCTK_WARN(0, "The grid does not include y==0"); + + for (int i = 0; i < 2; i++) { + ptrdiff_t offset_left = ms->gh->cctk_nghostzones[2 * i]; + ptrdiff_t offset_right = offset_left; + + /* account for grid tapering on refined levels */ + if (cp->level > 0) + offset_right *= integrator_substeps * 2; + + domain_size[i] = ms->gh->cctk_lsh[2 * i] - offset_left - offset_right; + /* the outer boundary layer overlaps with the first ghost zone */ + if (cp->level > 0) + domain_size[i]++; + + cp->offset_left[i] = offset_left; + cp->offset_right[i] = offset_right; + } + + if (domain_size[0] != domain_size[1]) + CCTK_WARN(0, "The grid is non-square, only square grids are supported"); + + cp->solver = mg2d_solver_alloc(domain_size[0]); + if (!cp->solver) + CCTK_WARN(0, "Error allocating the solver"); + + cp->solver->step[0] = a_x[1] - a_x[0]; + cp->solver->step[1] = cp->solver->step[0]; + + cp->solver->fd_stencil = 1; + + cp->solver->boundaries[MG2D_BOUNDARY_0L]->type = MG2D_BC_TYPE_FIXDIFF; + cp->solver->boundaries[MG2D_BOUNDARY_1L]->type = MG2D_BC_TYPE_FIXDIFF; + cp->solver->boundaries[MG2D_BOUNDARY_0U]->type = MG2D_BC_TYPE_FIXVAL; + cp->solver->boundaries[MG2D_BOUNDARY_1U]->type = MG2D_BC_TYPE_FIXVAL; + + cp->solver->maxiter = 32; + cp->solver->tol = 1e-10; + cp->solver->nb_cycles = 2; + + /* initialize boundary values to zero, + * non-zero values on the outer boundaries of refined levels are filled in elsewhere */ + for (int i = 0; i < 4; i++) + memset(cp->solver->boundaries[i]->val, 0, sizeof(double) * cp->solver->boundaries[i]->val_len); + + /* initialize the interpolation state */ + if (cp->level > 0) { + cp->boundary_val_stride[0] = gh->cctk_lsh[0] - cp->offset_left[0]; + cp->boundary_val_stride[1] = cp->solver->domain_size; + + cp->nb_boundary_interp_vals = cp->boundary_val_stride[0] * cp->offset_right[1] + + cp->boundary_val_stride[1] * cp->offset_right[0]; + + ret = posix_memalign((void**)&cp->boundary_interp_vals, 32, + cp->nb_boundary_interp_vals * sizeof(*cp->boundary_interp_vals)); + if (ret != 0) + CCTK_WARN(0, "Error allocating arrays for boundary values"); + cp->boundary_val[0] = cp->boundary_interp_vals; + cp->boundary_val[1] = cp->boundary_interp_vals + cp->boundary_val_stride[0] * cp->offset_right[1]; + + for (int i = 0; i < 3; i++) { + ret = posix_memalign((void**)&cp->interp_coords[i], 32, + cp->nb_boundary_interp_vals * sizeof(*cp->interp_coords[i])); + if (ret != 0) + CCTK_WARN(0, "Error allocating interp coords"); + } + + for (ptrdiff_t idx_row = 0; idx_row < cp->offset_right[1]; idx_row++) + for (ptrdiff_t i = 0; i < cp->boundary_val_stride[0]; i++) { + ptrdiff_t idx_dst = idx_row * cp->boundary_val_stride[0] + i; + ptrdiff_t idx_src_x = CCTK_GFINDEX3D(gh, cp->offset_left[0] + i, 0, 0); + ptrdiff_t idx_src_z = CCTK_GFINDEX3D(gh, 0, 0, gh->cctk_lsh[2] - cp->offset_right[1] + idx_row); + + cp->interp_coords[0][idx_dst] = a_x[idx_src_x]; + cp->interp_coords[1][idx_dst] = 0.0; + cp->interp_coords[2][idx_dst] = a_z[idx_src_z]; + } + + for (ptrdiff_t idx_row = 0; idx_row < cp->offset_right[0]; idx_row++) + for (ptrdiff_t i = 0; i < cp->boundary_val_stride[1]; i++) { + ptrdiff_t idx_dst = (cp->boundary_val[1] - cp->boundary_val[0]) + + idx_row * cp->boundary_val_stride[1] + i; + ptrdiff_t idx_src_x = CCTK_GFINDEX3D(gh, gh->cctk_lsh[0] - cp->offset_right[0] + idx_row, + 0, 0); + ptrdiff_t idx_src_z = CCTK_GFINDEX3D(gh, 0, 0, cp->offset_left[1] + i); + + cp->interp_coords[0][idx_dst] = a_x[idx_src_x]; + cp->interp_coords[1][idx_dst] = 0.0; + cp->interp_coords[2][idx_dst] = a_z[idx_src_z]; + } + +#if 0 + for (int i = 0; i < gh->cctk_lsh[0]; i++) { + cp->interp_coords[0][i] = a_x[CCTK_GFINDEX3D(gh, i, 0, 0)]; + cp->interp_coords[1][i] = 0.0; + cp->interp_coords[2][i] = a_z[CCTK_GFINDEX3D(gh, 0, 0, gh->cctk_lsh[2] - gh->cctk_nghostzones[2] - 1)]; + } + for (int i = 0; i < gh->cctk_lsh[2]; i++) { + int idx = gh->cctk_lsh[0] + i; + cp->interp_coords[0][idx] = a_x[CCTK_GFINDEX3D(gh, gh->cctk_lsh[0] - gh->cctk_nghostzones[0] - 1, 0, 0)]; + cp->interp_coords[1][idx] = 0.0; + cp->interp_coords[2][idx] = a_z[CCTK_GFINDEX3D(gh, 0, 0, i)]; + } +#endif + + cp->interp_values_codes[0] = CCTK_VARIABLE_REAL; + cp->interp_var_indices[0] = CCTK_VarIndex("ML_BSSN::alpha"); + if (cp->interp_var_indices[0] < 0) + CCTK_WARN(0, "Error getting the index of lapse variable"); + cp->interp_operation_codes[0] = 0; + cp->interp_operation_indices[0] = 0; + + cp->coord_system = CCTK_CoordSystemHandle("cart3d"); + if (cp->coord_system < 0) + CCTK_WARN(0, "Error getting the coordinate system"); + + cp->interp_operator = CCTK_InterpHandle("Lagrange polynomial interpolation (tensor product)"); + if (cp->interp_operator < 0) + CCTK_WARN(0, "Error getting the interpolation operator"); + + cp->interp_params = Util_TableCreateFromString("order=4 want_global_mode=1"); + if (cp->interp_params < 0) + CCTK_WARN(0, "Error creating interpolation parameters table"); + + ret = Util_TableSetIntArray(cp->interp_params, ARRAY_ELEMS(cp->interp_operation_codes), + cp->interp_operation_codes, "operation_codes"); + if (ret < 0) + CCTK_WARN(0, "Error setting operation codes"); + + ret = Util_TableSetIntArray(cp->interp_params, ARRAY_ELEMS(cp->interp_operation_indices), + cp->interp_operation_indices, "operand_indices"); + if (ret < 0) + CCTK_WARN(0, "Error setting operand indices"); + } + + ms->nb_patches++; + return cp; +} + +static int context_init(cGH *gh, MSMGContext **ctx) +{ + MSMGContext *ms; + int ret; + + ms = calloc(1, sizeof(*ms)); + if (!ms) + return -ENOMEM; + + ms->gh = gh; + + *ctx = ms; + + return 0; +} + +static void context_free(MSMGContext **pms) +{ + MSMGContext *ms = *pms; + + if (!ms) + return; + + for (int i = 0; i < ms->nb_patches; i++) + coord_patch_free(&ms->patches[i]); + free(ms->patches); + + free(ms); + *pms = NULL; +} + +static void fill_eq_coeffs(cGH *gh, CoordPatch *cp) +{ + int ret; + + const ptrdiff_t stride_z = CCTK_GFINDEX3D(gh, 0, 0, 1); + + double *a_x = CCTK_VarDataPtr(gh, 0, "grid::x"); + double *a_z = CCTK_VarDataPtr(gh, 0, "grid::z"); + + const double dx = a_x[1] - a_x[0]; + const double dz = a_z[stride_z] - a_z[0]; + + double *a_gtxx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt11"); + double *a_gtyy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt22"); + double *a_gtzz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt33"); + double *a_gtxy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt12"); + double *a_gtxz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt13"); + double *a_gtyz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt23"); + double *a_Atxx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At11"); + double *a_Atyy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At22"); + double *a_Atzz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At33"); + double *a_Atxy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At12"); + double *a_Atxz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At13"); + double *a_Atyz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At23"); + double *a_phi = CCTK_VarDataPtr(gh, 0, "ML_BSSN::phi"); + double *a_Xtx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt1"); + double *a_Xtz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt3"); + + for (int idx_z = 0; idx_z < cp->solver->domain_size; idx_z++) + for (int idx_x = 0; idx_x < cp->solver->domain_size; idx_x++) { + const int idx_src = CCTK_GFINDEX3D(gh, idx_x + gh->cctk_nghostzones[0], cp->y_idx, idx_z + gh->cctk_nghostzones[2]); + const int idx_dc = idx_z * cp->solver->diff_coeffs_stride + idx_x; + const int idx_rhs = idx_z * cp->solver->rhs_stride + idx_x; + + const double x = a_x[idx_src]; + + const double gtxx = a_gtxx[idx_src]; + const double gtyy = a_gtyy[idx_src]; + const double gtzz = a_gtzz[idx_src]; + const double gtxy = a_gtxy[idx_src]; + const double gtxz = a_gtxz[idx_src]; + const double gtyz = a_gtyz[idx_src]; + const double Atxx = a_Atxx[idx_src]; + const double Atyy = a_Atyy[idx_src]; + const double Atzz = a_Atzz[idx_src]; + const double Atxy = a_Atxy[idx_src]; + const double Atxz = a_Atxz[idx_src]; + const double Atyz = a_Atyz[idx_src]; + const double phi = a_phi[idx_src]; + const double Xtx = a_Xtx[idx_src]; + const double Xtz = a_Xtz[idx_src]; + + const double phi_dx = (a_phi[idx_src + 1] - a_phi[idx_src - 1]) / (2.0 * dx); + const double phi_dz = (a_phi[idx_src + stride_z] - a_phi[idx_src - stride_z]) / (2.0 * dz); + + const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz); + + const double At[3][3] = { + { Atxx, Atxy, Atxz }, + { Atxy, Atyy, Atyz }, + { Atxz, Atyz, Atzz }}; + + + double Xx, Xz, k2; + + double Am[3][3], gtu[3][3]; + + 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]; + + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) { + double val = 0.0; + for (int l = 0; l < 3; l++) + val += gtu[j][l] * At[l][k]; + Am[j][k] = val; + } + + // K_{ij} K^{ij} + k2 = 0.0; + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + k2 += Am[j][k] * Am[k][j]; + + Xx = SQR(phi) * (Xtx + (phi_dx * gtu[0][0] + phi_dz * gtu[0][2]) / phi); + Xz = SQR(phi) * (Xtz + (phi_dx * gtu[0][2] + phi_dz * gtu[2][2]) / phi); + + cp->solver->diff_coeffs[MG2D_DIFF_COEFF_20][idx_dc] = SQR(phi) * (gtu[0][0] + ((idx_x == 0) ? gtu[1][1] : 0.0)); + cp->solver->diff_coeffs[MG2D_DIFF_COEFF_02][idx_dc] = SQR(phi) * gtu[2][2]; + cp->solver->diff_coeffs[MG2D_DIFF_COEFF_11][idx_dc] = SQR(phi) * gtu[0][2] * 2; + cp->solver->diff_coeffs[MG2D_DIFF_COEFF_10][idx_dc] = -Xx + (idx_x ? SQR(phi) * gtu[1][1] / x : 0.0); + cp->solver->diff_coeffs[MG2D_DIFF_COEFF_01][idx_dc] = -Xz; + cp->solver->diff_coeffs[MG2D_DIFF_COEFF_00][idx_dc] = -k2; + cp->solver->rhs[idx_rhs] = k2; + } + + if (cp->level > 0) { + ret = CCTK_InterpGridArrays(gh, 3, cp->interp_operator, cp->interp_params, + cp->coord_system, cp->nb_boundary_interp_vals, CCTK_VARIABLE_REAL, + (const void * const *)cp->interp_coords, ARRAY_ELEMS(cp->interp_var_indices), + cp->interp_var_indices, ARRAY_ELEMS(cp->interp_values_codes), cp->interp_values_codes, + (void * const *)&cp->boundary_interp_vals); + if (ret < 0) + CCTK_WARN(0, "Error interpolating"); + + for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++) + cp->solver->boundaries[MG2D_BOUNDARY_0U]->val[i] = cp->boundary_val[0][i] - 1.0; + for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++) + cp->solver->boundaries[MG2D_BOUNDARY_1U]->val[i] = cp->boundary_val[1][i] - 1.0; + } +} + +void maximal_slicing_axi_mg(CCTK_ARGUMENTS) +{ + static MSMGContext *ms; + + CoordPatch *cp; + + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + double *coeffs = NULL; + int i, ret; + + /* on the first run, init the solver */ + if (!ms) + context_init(cctkGH, &ms); + + fprintf(stderr, "ms mg solve: level %d time %g\n", ctz(ms->gh->cctk_levfac[0]), ms->gh->cctk_time); + + cp = get_coord_patch(ms); + + fill_eq_coeffs(ms->gh, cp); + + CCTK_TimerStart("MaximalSlicingAxiMG_Solve"); + ret = mg2d_solve(cp->solver); + if (ret < 0) + CCTK_WARN(0, "Error solving the maximal slicing equation"); + CCTK_TimerStop("MaximalSlicingAxiMG_Solve"); + + + CCTK_TimerStart("MaximalSlicingAxiMG_Copy"); +//#pragma omp parallel for + for (int j = 0; j < cp->solver->domain_size; j++) + for (int i = 0; i < cp->solver->domain_size; i++) { + const int idx_dst = CCTK_GFINDEX3D(ms->gh, i + cp->offset_left[0], cp->y_idx, j + cp->offset_left[1]); + const int idx_src = j * cp->solver->u_stride + i; + alpha[idx_dst] = 1.0 + cp->solver->u[idx_src]; + } + if (cp->level > 0) { + /* fill in the interpolated outer buffer points */ + for (int idx_z = 0; idx_z < cp->offset_right[1]; idx_z++) + for (int idx_x = 0; idx_x < cp->boundary_val_stride[0]; idx_x++) { + const ptrdiff_t idx_src = idx_z * cp->boundary_val_stride[0] + idx_x; + const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, cp->offset_left[0] + idx_x, cp->y_idx, ms->gh->cctk_lsh[2] - cp->offset_right[1] + idx_z); + alpha[idx_dst] = cp->boundary_val[0][idx_src]; + } + for (int idx_z = 0; idx_z < cp->boundary_val_stride[1]; idx_z++) + for (int idx_x = 0; idx_x < cp->offset_right[0]; idx_x++) { + const ptrdiff_t idx_src = idx_x * cp->boundary_val_stride[1] + idx_z; + const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, ms->gh->cctk_lsh[0] - cp->offset_right[0] + idx_x, cp->y_idx, + cp->offset_left[1] + idx_z); + alpha[idx_dst] = cp->boundary_val[1][idx_src]; + } + } + + /* fill in the axial symmetry ghostpoints by mirroring */ + for (int idx_z = cp->offset_left[1]; idx_z < ms->gh->cctk_lsh[2]; idx_z++) + for (int idx_x = 0; idx_x < cp->offset_left[0]; idx_x++) { + const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, idx_x, cp->y_idx, idx_z); + const ptrdiff_t idx_src = CCTK_GFINDEX3D(ms->gh, 2 * cp->offset_left[0] - idx_x, cp->y_idx, idx_z); + alpha[idx_dst] = alpha[idx_src]; + } + for (int idx_z = 0; idx_z < cp->offset_left[1]; idx_z++) + for (int idx_x = 0; idx_x < ms->gh->cctk_lsh[0]; idx_x++) { + const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, idx_x, cp->y_idx, idx_z); + const ptrdiff_t idx_src = CCTK_GFINDEX3D(ms->gh, idx_x, cp->y_idx, 2 * cp->offset_left[1] - idx_z); + alpha[idx_dst] = alpha[idx_src]; + } + CCTK_TimerStop("MaximalSlicingAxiMG_Copy"); +} diff --git a/src/register.c b/src/register.c new file mode 100644 index 0000000..9bc64cd --- /dev/null +++ b/src/register.c @@ -0,0 +1,18 @@ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "Slicing.h" + +int maximal_slicing_axi_mg_register(void) +{ + Einstein_RegisterSlicing("maximal_axi_mg"); + return 0; +} + +void maximal_slicing_axi_mg_register_mol(CCTK_ARGUMENTS) +{ + MoLRegisterConstrained(CCTK_VarIndex("ADMBase::alp")); + + MoLRegisterSaveAndRestoreGroup(CCTK_GroupIndex("ADMBase::metric")); + MoLRegisterSaveAndRestoreGroup(CCTK_GroupIndex("ADMBase::curv")); +} -- cgit v1.2.3