From 8ee77841e5d32d3243213730e6d12a880200d840 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Tue, 25 Sep 2018 10:02:18 +0200 Subject: Rewrite using a new method. Based on 'Adaptive mesh refinement for coupled elliptic-hyperbolic problems' by Pretorius and Choptuik, 2006. --- interface.ccl | 8 + schedule.ccl | 27 ++- src/maximal_slicing_axi_mg.c | 434 +++++++++++++++++++++++++++++-------------- 3 files changed, 326 insertions(+), 143 deletions(-) diff --git a/interface.ccl b/interface.ccl index 43e512c..124823e 100644 --- a/interface.ccl +++ b/interface.ccl @@ -13,3 +13,11 @@ USES FUNCTION MoLNumIntegratorSubsteps REQUIRES FUNCTION MoLRegisterConstrained REQUIRES FUNCTION MoLRegisterSaveAndRestore REQUIRES FUNCTION MoLRegisterSaveAndRestoreGroup + +public: +REAL lapse_mg TYPE=GF TIMELEVELS=2 +REAL lapse_mg_eval TYPE=GF tags='Prolongation="None"' +REAL lapse_prev0 TYPE=GF tags='Prolongation="None"' +REAL lapse_prev1 TYPE=GF tags='Prolongation="None"' +REAL lapse_prev0_time TYPE=ARRAY DIM=1 SIZE=32 +REAL lapse_prev1_time TYPE=ARRAY DIM=1 SIZE=32 diff --git a/schedule.ccl b/schedule.ccl index 246aa63..0b992b9 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -1,10 +1,26 @@ # Schedule definitions for thorn MaximalSlicingAxiMG # if (CCTK_Equals(lapse_evolution_method, "maximal_axi_mg")) { + SCHEDULE msa_mg_inithist in CCTK_INITIAL { + LANG: C + } "" - SCHEDULE maximal_slicing_axi_mg IN ML_BSSN_evolCalcGroup BEFORE ML_BSSN_RHS { + SCHEDULE msa_mg_init IN CCTK_BASEGRID { LANG: C - } "Maximal slicing in axisymmetry" + } "" + + SCHEDULE msa_mg_eval IN MoL_CalcRHS BEFORE ML_BSSN_evolCalcGroup { + LANG: C + } "" + + SCHEDULE msa_mg_sync IN CCTK_POSTSTEP BEFORE msa_mg_solve { + SYNC: lapse_mg + LANG: C + } "" + + SCHEDULE msa_mg_solve IN CCTK_POSTSTEP { + LANG: C + } "" #SCHEDULE maximal_slicing_axi IN MoL_PostStepModify { # LANG: C @@ -21,4 +37,11 @@ if (CCTK_Equals(lapse_evolution_method, "maximal_axi_mg")) { SCHEDULE maximal_slicing_axi_mg_register_mol IN MoL_Register { LANG: C } "" + + STORAGE: lapse_mg[2] + STORAGE: lapse_mg_eval + STORAGE: lapse_prev0 + STORAGE: lapse_prev1 + STORAGE: lapse_prev0_time + STORAGE: lapse_prev1_time } diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c index 974ef86..0838d6a 100644 --- a/src/maximal_slicing_axi_mg.c +++ b/src/maximal_slicing_axi_mg.c @@ -21,12 +21,14 @@ #define SQR(x) ((x) * (x)) #define ARRAY_ELEMS(x) (sizeof(x) / sizeof(*x)) +#define CPINDEX(cp, i, j, k) ((k * cp->grid_size[1] + j) * cp->grid_size[0] + i) + +typedef struct MaximalSlicingContext MaximalSlicingContext; + /* precomputed values for a given refined grid */ typedef struct CoordPatch { - double origin[3]; - int delta[3]; - int size[3]; int level; + ptrdiff_t grid_size[3]; MG2DContext *solver; @@ -37,21 +39,6 @@ typedef struct CoordPatch { 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 @@ -60,9 +47,7 @@ typedef struct CoordPatch { * [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 { @@ -70,8 +55,17 @@ typedef struct MSMGContext { CoordPatch *patches; int nb_patches; + + MaximalSlicingContext *ps_solver; } MSMGContext; +int msa_context_init(cGH *gh, MaximalSlicingContext **ctx, int basis_order0, int basis_order1, double scale_factor, double filter_power); +void msa_context_free(MaximalSlicingContext **ctx); +void msa_lapse_solve(MaximalSlicingContext *ctx); +void msa_lapse_eval(MaximalSlicingContext *ctx, double *dst); + +static MSMGContext *ms; + static int ctz(int a) { int ret = 0; @@ -92,10 +86,11 @@ static void coord_patch_free(CoordPatch *cp) mg2d_solver_free(&cp->solver); } -static CoordPatch *get_coord_patch(MSMGContext *ms) +static CoordPatch *get_coord_patch(MSMGContext *ms, int level) { cGH *gh = ms->gh; + const size_t grid_size = gh->cctk_lsh[2] * gh->cctk_lsh[1] * gh->cctk_lsh[0]; 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"); @@ -108,15 +103,7 @@ static CoordPatch *get_coord_patch(MSMGContext *ms) 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]) + if (cp->level == level) return cp; } @@ -126,17 +113,17 @@ static CoordPatch *get_coord_patch(MSMGContext *ms) 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]); + cp->grid_size[0] = gh->cctk_lsh[0]; + cp->grid_size[1] = gh->cctk_lsh[1]; + cp->grid_size[2] = gh->cctk_lsh[2]; - for (i = 0; i < cp->size[1]; i++) + for (i = 0; i < gh->cctk_lsh[1]; i++) if (fabs(a_y[CCTK_GFINDEX3D(gh, 0, i, 0)]) < 1e-8) { cp->y_idx = i; break; } - if (i == cp->size[1]) + if (i == gh->cctk_lsh[1]) CCTK_WARN(0, "The grid does not include y==0"); for (int i = 0; i < 2; i++) { @@ -146,6 +133,11 @@ static CoordPatch *get_coord_patch(MSMGContext *ms) /* account for grid tapering on refined levels */ if (cp->level > 0) offset_right *= integrator_substeps * 2; + ///* the outer boundary layer overlaps with the first ghost zone */ + //if (cp->level > 0) + // offset_right++; + if (cp->level > 0) + offset_right--; domain_size[i] = ms->gh->cctk_lsh[2 * i] - offset_left - offset_right; /* the outer boundary layer overlaps with the first ghost zone */ @@ -182,98 +174,9 @@ static CoordPatch *get_coord_patch(MSMGContext *ms) 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_TableSetInt(cp->interp_params, cp->level, "max_reflevel"); - if (ret < 0) - CCTK_WARN(0, "Error setting maximum reflevel for interpolation"); - - 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++; @@ -363,8 +266,13 @@ static void fill_eq_coeffs(cGH *gh, CoordPatch *cp) const double Xtx = a_Xtx[idx_src]; const double Xtz = a_Xtz[idx_src]; +#if 1 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); +#else + const double phi_dx = (-1.0 * a_phi[idx_src + 2] + 8.0 * a_phi[idx_src + 1] - 8.0 * a_phi[idx_src - 1] + a_phi[idx_src - 1]) / (12.0 * dx); + const double phi_dz = (-1.0 * a_phi[idx_src + 2 * stride_z] + 8.0 * a_phi[idx_src + 1 * stride_z] - 8.0 * a_phi[idx_src - 1 * stride_z] + a_phi[idx_src - 1 * stride_z]) / (12.0 * dz); +#endif const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz); @@ -413,39 +321,259 @@ static void fill_eq_coeffs(cGH *gh, CoordPatch *cp) 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); +static void solution_to_grid(CoordPatch *cp, double *dst) +{ + for (int j = 0; j < cp->solver->domain_size; j++) + for (int i = 0; i < cp->solver->domain_size; i++) { + const ptrdiff_t idx_dst = CPINDEX(cp, i + cp->offset_left[0], cp->y_idx, j + cp->offset_left[1]); + const int idx_src = j * cp->solver->u_stride + i; + dst[idx_dst] = 1.0 + cp->solver->u[idx_src]; + } + + /* fill in the axial symmetry ghostpoints by mirroring */ + for (int idx_z = cp->offset_left[1]; idx_z < cp->grid_size[2]; idx_z++) + for (int idx_x = 0; idx_x < cp->offset_left[0]; idx_x++) { + const ptrdiff_t idx_dst = CPINDEX(cp, idx_x, cp->y_idx, idx_z); + const ptrdiff_t idx_src = CPINDEX(cp, 2 * cp->offset_left[0] - idx_x, cp->y_idx, idx_z); + dst[idx_dst] = dst[idx_src]; + } + for (int idx_z = 0; idx_z < cp->offset_left[1]; idx_z++) + for (int idx_x = 0; idx_x < cp->grid_size[0]; idx_x++) { + const ptrdiff_t idx_dst = CPINDEX(cp, idx_x, cp->y_idx, idx_z); + const ptrdiff_t idx_src = CPINDEX(cp, idx_x, cp->y_idx, 2 * cp->offset_left[1] - idx_z); + dst[idx_dst] = dst[idx_src]; + } +} + +static void solve_ps(MSMGContext *ctx, double *dst) +{ + if (!ctx->ps_solver) + msa_context_init(ctx->gh, &ctx->ps_solver, 60, 60, 8.0, 64.0); + + msa_lapse_solve(ctx->ps_solver); + msa_lapse_eval(ctx->ps_solver, dst); +} + +static int history_initialized(double *lapse_prev0_time, double *lapse_prev1_time, int nb_levels) +{ + for (int i = 0; i < nb_levels; i++) { + if (lapse_prev0_time[i] == DBL_MAX || lapse_prev1_time[i] == DBL_MAX) + return 0; + } + return 1; +} + +void msa_mg_eval(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + int nb_levels; + int nb_levels_type, ret; + + const size_t grid_size = cctkGH->cctk_lsh[2] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[0]; + const double t = cctkGH->cctk_time; + const int reflevel = ctz(cctkGH->cctk_levfac[0]); + double time_interp_step, fact0, fact1; + + fprintf(stderr, "evaluating lapse at rl=%d, t=%g\n", reflevel, t); + + nb_levels = *(int*)CCTK_ParameterGet("num_levels_1", "CarpetRegrid2", &nb_levels_type); + + if (!history_initialized(lapse_prev0_time, lapse_prev1_time, nb_levels)) { + solve_ps(ms, lapse_mg_eval); + memcpy(alpha, lapse_mg_eval, sizeof(*alpha) * grid_size); + fprintf(stderr, "pseudospectral solve\n"); + return; + } + + if (lapse_prev1_time[reflevel] <= lapse_prev0_time[reflevel] || + fabs(lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel]) < 1e-14) + CCTK_WARN(0, "Invalid past value times"); + + time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel]; + fact0 = (lapse_prev1_time[reflevel] - t) / time_interp_step; + fact1 = (t - lapse_prev0_time[reflevel]) / time_interp_step; + + if (reflevel < nb_levels - 1) { + /* on coarse levels use extrapolated past solutions */ + fprintf(stderr, "extrapolating 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 < grid_size; i++) + lapse_mg_eval[i] = fact0 * lapse_prev0[i] + fact1 * lapse_prev1[i]; + } else { + /* on the finest level, use the extrapolation only for the boundary + * values and solve for the interior */ + CoordPatch *cp = get_coord_patch(ms, reflevel); + + fill_eq_coeffs(ms->gh, cp); + + fprintf(stderr, "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0); + + for (ptrdiff_t idx_z = 0; idx_z < cp->offset_right[1]; idx_z++) + for (ptrdiff_t idx_x = 0; idx_x < cp->boundary_val_stride[0]; idx_x++) { + const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, idx_x + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1] + idx_z); + lapse_mg_eval[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx]; + } + 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 = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0] + idx_x, cp->y_idx, + cp->offset_left[1] + idx_z); + 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 (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; + } + + fprintf(stderr, "mg solve\n"); + + ret = mg2d_solve(cp->solver); if (ret < 0) - CCTK_WARN(0, "Error interpolating"); + CCTK_WARN(0, "Error solving the maximal slicing equation"); - 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; + solution_to_grid(cp, lapse_mg_eval); } + memcpy(alpha, lapse_mg_eval, sizeof(*alpha) * grid_size); } -void maximal_slicing_axi_mg(CCTK_ARGUMENTS) +void msa_mg_solve(CCTK_ARGUMENTS) { - static MSMGContext *ms; + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; CoordPatch *cp; + int nb_levels, reflevel_top, ts_tmp; + int nb_levels_type, ret; + + const size_t grid_size = cctkGH->cctk_lsh[2] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[0]; + const double t = cctkGH->cctk_time; + const int timestep = lrint(t * cctkGH->cctk_levfac[0] / cctkGH->cctk_delta_time); + const int reflevel = ctz(cctkGH->cctk_levfac[0]); + + double *lapse_mg1; + + fprintf(stderr, "solve lapse at rl=%d, t=%g; step %d\n", reflevel, t, timestep); + + nb_levels = *(int*)CCTK_ParameterGet("num_levels_1", "CarpetRegrid2", &nb_levels_type); + + reflevel_top = reflevel; + ts_tmp = timestep; + while (reflevel_top > 0 && !(ts_tmp % 2)) { + ts_tmp /= 2; + reflevel_top--; + } + + if (!history_initialized(lapse_prev0_time, lapse_prev1_time, nb_levels)) { + if (reflevel == 0 || !(timestep % 2)) { + if (lapse_prev1_time[reflevel] != DBL_MAX && + reflevel_top == 0) { + fprintf(stderr, "copy %g\n", lapse_prev1_time[reflevel]); + memcpy(lapse_prev0, lapse_prev1, grid_size * sizeof(*lapse_prev1)); + lapse_prev0_time[reflevel] = lapse_prev1_time[reflevel]; + } + + fprintf(stderr, "solve ps\n"); + solve_ps(ms, alpha); + memcpy(lapse_prev1, alpha, grid_size * sizeof(*lapse_prev0)); + lapse_prev1_time[reflevel] = t; + } + return; + } + + fprintf(stderr, "mg solve cur %d top %d\n", reflevel, reflevel_top); + + /* fill in the equation coefficients */ + cp = get_coord_patch(ms, reflevel); + fill_eq_coeffs(cctkGH, cp); + + 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; + + fprintf(stderr, "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 (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]; + } + } + /* 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 (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; + } + + /* do the elliptic solve */ + ret = mg2d_solve(cp->solver); + if (ret < 0) + CCTK_WARN(0, "Error solving the maximal slicing equation"); + + /* export the result */ + solution_to_grid(cp, lapse_mg); + + lapse_mg1 = CCTK_VarDataPtr(cctkGH, 1, "MaximalSlicingAxiMG::lapse_mg"); + memcpy(lapse_mg1, lapse_mg, grid_size * sizeof(*lapse_mg1)); + + /* update lapse history for extrapolation */ + if (reflevel == 0 || !(timestep % 2)) { + const double vel_fact = 1.0 / (1 << (reflevel - reflevel_top)); + + for (size_t j = 0; j < grid_size; j++) { + const double sol_new = lapse_mg[j]; + const double delta = sol_new - lapse_mg_eval[j]; + lapse_prev0[j] = lapse_prev1[j] + delta - delta * vel_fact; + lapse_prev1[j] = sol_new; + alpha[j] = sol_new; + } + lapse_prev0_time[reflevel] = lapse_prev1_time[reflevel]; + lapse_prev1_time[reflevel] = cctkGH->cctk_time; + } +} + +void msa_mg_sync(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + const int reflevel = ctz(cctkGH->cctk_levfac[0]); + fprintf(stderr, "\nsync %g %d\n\n", cctkGH->cctk_time, reflevel); + + return; +} + +#if 0 +static void maximal_slicing_axi_mg_old(CCTK_ARGUMENTS) +{ + 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); @@ -499,3 +627,27 @@ void maximal_slicing_axi_mg(CCTK_ARGUMENTS) } CCTK_TimerStop("MaximalSlicingAxiMG_Copy"); } +#endif + +void msa_mg_init(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + context_init(cctkGH, &ms); +} + +void msa_mg_inithist(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + const size_t grid_size = cctkGH->cctk_lsh[2] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[0]; + + for (size_t i = 0; i < grid_size; i++) + lapse_mg[i] = 1.0; + + for (int i = 0; i < 32; i++) { + lapse_prev0_time[i] = DBL_MAX; + lapse_prev1_time[i] = DBL_MAX; + } +} -- cgit v1.2.3