From c7c5610d529844104886d8e4a2fc44510a82bc7b Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 27 Jun 2019 16:41:48 +0200 Subject: Use extrapolated solution as initial guess for outermost solve level. --- src/maximal_slicing_axi_mg.c | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c index d9d1508..23051df 100644 --- a/src/maximal_slicing_axi_mg.c +++ b/src/maximal_slicing_axi_mg.c @@ -1339,7 +1339,7 @@ void msa_mg_solve(CCTK_ARGUMENTS) if (reflevel > 0) { if (timestep % 2) { /* outer-most level for this solve, use extrapolated lapse as the - * boundary condition */ + * boundary condition and initial guess */ double time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel]; double fact0, fact1; @@ -1357,24 +1357,18 @@ 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); - if (cp->solver->local_start[1] + cp->solver->local_size[1] == cp->solver->domain_size) { #pragma omp parallel for - for (int j = 0; j < cp->solver->fd_stencil; j++) { - for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->local_size[0] + j; i++) { - const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cp->offset_left[1] + cp->solver->local_size[1] - 1 + j); - lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx]; - } - } - } - if (cp->solver->local_start[0] + cp->solver->local_size[0] == cp->solver->domain_size) { + for (size_t i = 0; i < grid_size; i++) + lapse_mg[i] = fact0 * lapse_prev0[i] + fact1 * lapse_prev1[i]; + #pragma omp parallel for - for (int j = 0; j < cp->solver->fd_stencil; j++) { - for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->local_size[1] + j; i++) { - const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cp->offset_left[0] + cp->solver->local_size[0] - 1 + j, 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->local_size[1]; j++) + for (int i = 0; i < cp->solver->local_size[0]; i++) { + const ptrdiff_t idx_src = CPINDEX(cp, i + cp->offset_left[0], cp->y_idx, j + cp->offset_left[1]); + const int idx_dst = j * cp->solver->u_stride + i; + cp->solver->u[idx_dst] = lapse_mg[idx_src] - 1.0; } + } else { /* use the solution from the coarser level as the initial guess */ CoordPatch *cp1 = get_coord_patch(ms, reflevel - 1); -- cgit v1.2.3