From 040e858d1c4aaede251981f1a204ff5fd5c0c94c Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 15 Feb 2019 18:30:25 +0100 Subject: Stop filling the boundary values for coarsest level. --- src/maximal_slicing_axi_mg.c | 66 +++++++++++++++++++++++--------------------- 1 file changed, 34 insertions(+), 32 deletions(-) diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c index 56a7c94..c83c2fe 100644 --- a/src/maximal_slicing_axi_mg.c +++ b/src/maximal_slicing_axi_mg.c @@ -633,7 +633,7 @@ void msa_mg_eval(CCTK_ARGUMENTS) } for (int j = 0; j < cp->solver->fd_stencil; j++) { - MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_0U]; + 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, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1] + j); @@ -641,7 +641,7 @@ void msa_mg_eval(CCTK_ARGUMENTS) } } for (int j = 0; j < cp->solver->fd_stencil; j++) { - MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_1U]; + 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, cctkGH->cctk_lsh[0] - cp->offset_right[0] + j, cp->y_idx, cp->offset_left[1] + i); @@ -736,48 +736,50 @@ void msa_mg_solve(CCTK_ARGUMENTS) ms->time_solve_fill_coeffs += gettime() - start; start = gettime(); - if (reflevel == 0 || timestep % 2) { - /* outer-most level for this solve, use extrapolated lapse as the - * boundary condition */ - double time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel]; - double fact0 = (lapse_prev1_time[reflevel] - t) / time_interp_step; - double fact1 = (t - lapse_prev0_time[reflevel]) / time_interp_step; - - LOGDEBUG( "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0); + if (reflevel > 0) { + if (timestep % 2) { + /* outer-most level for this solve, use extrapolated lapse as the + * boundary condition */ + double time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel]; + double fact0 = (lapse_prev1_time[reflevel] - t) / time_interp_step; + double fact1 = (t - lapse_prev0_time[reflevel]) / time_interp_step; + + LOGDEBUG( "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0); + + 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 (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 (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, 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]; + dst[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, 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]; + dst[i] = lapse_mg[idx] - 1.0; } } } - /* if the condition above was false, then lapse_mg should be filled by - * prolongation from the coarser level */ - - /* fill the solver boundary conditions */ - 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 (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; /* do the elliptic solve */ -- cgit v1.2.3