From 9bb2871198d34564bd9f9ece01d4e5619cf50b49 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 29 Nov 2018 18:00:17 +0100 Subject: Implement higher-order finite differences. --- param.ccl | 4 +++ src/maximal_slicing_axi_mg.c | 72 ++++++++++++++++++++++++++++++-------------- 2 files changed, 53 insertions(+), 23 deletions(-) diff --git a/param.ccl b/param.ccl index e6b39e2..1718290 100644 --- a/param.ccl +++ b/param.ccl @@ -8,3 +8,7 @@ EXTENDS KEYWORD lapse_evolution_method } RESTRICTED: +CCTK_INT fd_stencil "finite differencing stencil" +{ + 1: :: "" +} 1 diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c index b5255cc..ae851db 100644 --- a/src/maximal_slicing_axi_mg.c +++ b/src/maximal_slicing_axi_mg.c @@ -67,6 +67,8 @@ typedef struct CoordPatch { typedef struct MSMGContext { cGH *gh; + int fd_stencil; + CoordPatch *patches; int nb_patches; @@ -203,7 +205,7 @@ static CoordPatch *get_coord_patch(MSMGContext *ms, int level) 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->fd_stencil = ms->fd_stencil; cp->solver->boundaries[MG2D_BOUNDARY_0L]->type = MG2D_BC_TYPE_FIXDIFF; cp->solver->boundaries[MG2D_BOUNDARY_1L]->type = MG2D_BC_TYPE_FIXDIFF; @@ -216,8 +218,11 @@ static CoordPatch *get_coord_patch(MSMGContext *ms, int level) /* 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); + for (int i = 0; i < 4; i++) { + for (size_t j = 0; j < cp->solver->fd_stencil; j++) + memset(cp->solver->boundaries[i]->val + j * cp->solver->boundaries[i]->val_stride - j, + 0, sizeof(*cp->solver->boundaries[i]->val) * (domain_size[0] + 2 * j)); + } if (cp->level > 0) { cp->boundary_val_stride[0] = gh->cctk_lsh[0] - cp->offset_left[0]; @@ -276,7 +281,7 @@ static void print_stats(MSMGContext *ms) } } -static int context_init(cGH *gh, MSMGContext **ctx) +static int context_init(cGH *gh, int fd_stencil, MSMGContext **ctx) { MSMGContext *ms; int ret; @@ -286,6 +291,7 @@ static int context_init(cGH *gh, MSMGContext **ctx) return -ENOMEM; ms->gh = gh; + ms->fd_stencil = fd_stencil; *ctx = ms; @@ -537,13 +543,21 @@ void msa_mg_eval(CCTK_ARGUMENTS) lapse_mg_eval[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx]; } - for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++) { - const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1]); - cp->solver->boundaries[MG2D_BOUNDARY_0U]->val[i] = lapse_mg_eval[idx] - 1.0; + for (int j = 0; j < cp->solver->fd_stencil; j++) { + MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_0U]; + double *dst = bnd->val + j * bnd->val_stride; + for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) { + const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1] + j); + dst[i] = lapse_mg_eval[idx] - 1.0; + } } - for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++) { - const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0], cp->y_idx, cp->offset_left[1] + i); - cp->solver->boundaries[MG2D_BOUNDARY_1U]->val[i] = lapse_mg_eval[idx] - 1.0; + for (int j = 0; j < cp->solver->fd_stencil; j++) { + MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_1U]; + double *dst = bnd->val + j * bnd->val_stride; + for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) { + const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0] + j, cp->y_idx, cp->offset_left[1] + i); + dst[i] = lapse_mg_eval[idx] - 1.0; + } } ms->time_fine_boundaries += gettime() - start; @@ -642,26 +656,38 @@ void msa_mg_solve(CCTK_ARGUMENTS) LOGDEBUG( "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0); - for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++) { - const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1]); - lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx]; + for (int j = 0; j < cp->solver->fd_stencil; j++) { + for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) { + const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1] + j); + lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx]; + } } - for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++) { - const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0], cp->y_idx, cp->offset_left[1] + i); - lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx]; + for (int j = 0; j < cp->solver->fd_stencil; j++) { + for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) { + const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0] + j, cp->y_idx, cp->offset_left[1] + i); + lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx]; + } } } /* if the condition above was false, then lapse_mg should be filled by * prolongation from the coarser level */ /* fill the solver boundary conditions */ - for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++) { - const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1]); - cp->solver->boundaries[MG2D_BOUNDARY_0U]->val[i] = lapse_mg[idx] - 1.0; + for (int j = 0; j < cp->solver->fd_stencil; j++) { + MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_0U]; + double *dst = bnd->val + j * bnd->val_stride; + for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) { + const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1] + j); + dst[i] = lapse_mg[idx] - 1.0; + } } - for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++) { - const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0], cp->y_idx, cp->offset_left[1] + i); - cp->solver->boundaries[MG2D_BOUNDARY_1U]->val[i] = lapse_mg[idx] - 1.0; + for (int j = 0; j < cp->solver->fd_stencil; j++) { + MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_1U]; + double *dst = bnd->val + j * bnd->val_stride; + for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->domain_size + j; i++) { + const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0] + j, cp->y_idx, cp->offset_left[1] + i); + dst[i] = lapse_mg[idx] - 1.0; + } } ms->time_solve_boundaries += gettime() - start; @@ -790,7 +816,7 @@ void msa_mg_init(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; - context_init(cctkGH, &ms); + context_init(cctkGH, fd_stencil, &ms); } void msa_mg_inithist(CCTK_ARGUMENTS) -- cgit v1.2.3