From bc7d230e517e25303680c4cbe367a15769d16e7c Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Tue, 10 Sep 2019 10:40:58 +0200 Subject: Smoothly zero the lapse source at the outer boundary. --- src/qms.c | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/qms.c b/src/qms.c index 1005290..d782d0c 100644 --- a/src/qms.c +++ b/src/qms.c @@ -362,8 +362,8 @@ static CoordPatch *get_coord_patch(QMSMGContext *ms, int level) solver->boundaries[MG2D_BOUNDARY_0L]->type = MG2D_BC_TYPE_REFLECT; solver->boundaries[MG2D_BOUNDARY_1L]->type = MG2D_BC_TYPE_REFLECT; - solver->boundaries[MG2D_BOUNDARY_0U]->type = level ? MG2D_BC_TYPE_FIXVAL : MG2D_BC_TYPE_FALLOFF; - solver->boundaries[MG2D_BOUNDARY_1U]->type = level ? MG2D_BC_TYPE_FIXVAL : MG2D_BC_TYPE_FALLOFF; + solver->boundaries[MG2D_BOUNDARY_0U]->type = MG2D_BC_TYPE_FIXVAL; + solver->boundaries[MG2D_BOUNDARY_1U]->type = MG2D_BC_TYPE_FIXVAL; solver->maxiter = ms->maxiter; solver->tol = ms->tol_residual_base / (solver->step[0] * solver->step[1]); @@ -1420,6 +1420,25 @@ void qms_mg_solve(CCTK_ARGUMENTS) W_pred1[idx] = (W_val1[idx] * fact1 + W_val0[idx] * fact0) * (1.0 - bnd_fact) + coarse_val * bnd_fact; } mirror(cp, W_pred1); + } else if (reflevel == 0) { + const int interior_size = cctk_lsh[0] - cp->offset_left[0] - cctk_nghostzones[0]; + + const double bnd_loc = x[CCTK_GFINDEX3D(ms->gh, cp->offset_left[0] + interior_size, 0, 0)]; + const double transition_start = bnd_loc * 0.7; + +#pragma omp parallel for + for (int idx = 0; idx < grid_size; idx++) { + const double x_val = x[idx]; + const double z_val = z[idx]; + const double coarse_val = W_pred1[idx]; + + const double r = sqrt(SQR(x_val) + SQR(z_val)); + const double dist = (r - transition_start) / (bnd_loc - transition_start); + + const double bnd_fact = 1.0 - smooth_step(dist); + + W_pred1[idx] = (W_val1[idx] * fact1 + W_val0[idx] * fact0) * bnd_fact; + } } else { #pragma omp parallel for for (int i = 0; i < grid_size; i++) -- cgit v1.2.3