#include #include #include #include #include #include #include #include #include #include #include #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" #include "cctk_Timers.h" #include "util_Table.h" #define SQR(x) ((x) * (x)) #define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0) #define MAX(x, y) ((x) > (y) ? (x) : (y)) #define MIN(x, y) ((x) > (y) ? (y) : (x)) #define ABS(x) ((x >= 0) ? (x) : -(x)) #define ARRAY_ELEMS(arr) (sizeof(arr) / sizeof(*arr)) #define CPINDEX(cp, i, j, k) ((k * cp->grid_size[1] + j) * cp->grid_size[0] + i) /* precomputed values for a given refined grid */ typedef struct CoordPatch { int level; ptrdiff_t grid_size[3]; MG2DContext *solver; int y_idx; ptrdiff_t offset_left[2]; int bnd_intercomp[2][2]; } CoordPatch; typedef struct QMSMGContext { cGH *gh; int log_level; int fd_stencil; int maxiter; int max_exact_size; int nb_cycles; int nb_relax_post; int nb_relax_pre; double tol_residual_base; double cfl_factor; int solve_level; CoordPatch *patches; int nb_patches; } QMSMGContext; static const struct { const char *str; enum MG2DLogLevel level; } log_levels[] = { { "fatal", MG2D_LOG_FATAL }, { "error", MG2D_LOG_ERROR }, { "warning", MG2D_LOG_WARNING }, { "info", MG2D_LOG_INFO }, { "verbose", MG2D_LOG_VERBOSE }, { "debug", MG2D_LOG_DEBUG }, }; static QMSMGContext *qms_context; static int ctz(int a) { int ret = 0; if (!a) return INT_MAX; while (!(a & 1)) { a >>= 1; ret++; } return ret; } static void log_callback(const MG2DContext *ctx, int level, const char *fmt, va_list vl) { QMSMGContext *ms = ctx->opaque; int target_level = ms->log_level; uint8_t buf[1024]; int ret; if (level > target_level) return; ret = snprintf(buf, sizeof(buf), "[%d] t=%g ", ctz(ms->gh->cctk_levfac[0]), ms->gh->cctk_time); if (ret >= sizeof(buf)) return; vsnprintf(buf + ret, sizeof(buf) - ret, fmt, vl); fputs(buf, stderr); } static CoordPatch *get_coord_patch(QMSMGContext *ms, int level) { cGH *gh = ms->gh; const char *omp_threads = getenv("OMP_NUM_THREADS"); const int integrator_substeps = MoLNumIntegratorSubsteps(); 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"); size_t solver_size[2]; CoordPatch *cp; MG2DContext *solver; int i, ret; for (int i = 0; i < ms->nb_patches; i++) { cp = &ms->patches[i]; if (cp->level == level) return cp; } /* create a new patch */ ms->patches = realloc(ms->patches, sizeof(*ms->patches) * (ms->nb_patches + 1)); cp = &ms->patches[ms->nb_patches]; memset(cp, 0, sizeof(*cp)); cp->level = level; 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 < gh->cctk_lsh[1]; i++) if (fabs(a_y[CCTK_GFINDEX3D(gh, 0, i, 0)]) < 1e-8) { cp->y_idx = i; break; } if (i == gh->cctk_lsh[1]) CCTK_WARN(0, "The grid does not include y==0"); for (int dir = 0; dir < 2; dir++) { ptrdiff_t offset_right; cp->offset_left[dir] = gh->cctk_nghostzones[2 * dir]; offset_right = gh->cctk_nghostzones[2 * dir]; if (cp->level > 0) { offset_right *= integrator_substeps * 2; offset_right -= 2; } solver_size[dir] = cp->grid_size[2 * dir] - cp->offset_left[dir] - offset_right; } if (solver_size[0] != solver_size[1]) CCTK_WARN(0, "Non-square domain"); solver = mg2d_solver_alloc(solver_size[0]); if (!solver) CCTK_WARN(0, "Error allocating the solver"); cp->solver = solver; solver->step[0] = a_x[1] - a_x[0]; solver->step[1] = a_z[CCTK_GFINDEX3D(gh, 0, 0, 1)] - a_z[0]; solver->fd_stencil = ms->fd_stencil; 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 = 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]); solver->nb_cycles = ms->nb_cycles; solver->nb_relax_pre = ms->nb_relax_pre; solver->nb_relax_post = ms->nb_relax_post; solver->cfl_factor = ms->cfl_factor; solver->max_exact_size = ms->max_exact_size; solver->opaque = ms; solver->log_callback = log_callback; if (omp_threads) solver->nb_threads = strtol(omp_threads, NULL, 0); if (solver->nb_threads <= 0) solver->nb_threads = 1; /* 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++) { int ci = mg2d_bnd_coord_idx(i); for (size_t j = 0; j < solver->fd_stencil; j++) { memset(solver->boundaries[i]->val + j * solver->boundaries[i]->val_stride - j, 0, sizeof(*solver->boundaries[i]->val) * (solver_size[!ci] + 2 * j)); } } ms->nb_patches++; return cp; } #define FD4(arr, stride, h_inv) \ ((-1.0 * arr[2 * stride] + 8.0 * arr[1 * stride] - 8.0 * arr[-1 * stride] + 1.0 * arr[-2 * stride]) * h_inv / 12.0) #define FD24(arr, stride, h_inv) \ ((-1.0 * arr[2 * stride] + 16.0 * arr[1 * stride] - 30.0 * arr[0] + 16.0 * arr[-1 * stride] - 1.0 * arr[-2 * stride]) * SQR(h_inv) / 12.0) static double diff4_dx(const double *arr, ptrdiff_t stride, double h_inv) { return FD4(arr, 1, h_inv); } static double diff4_dz(const double *arr, ptrdiff_t stride, double h_inv) { return FD4(arr, stride, h_inv); } static double diff4_dxx(const double *arr, ptrdiff_t stride, double h_inv) { return FD24(arr, 1, h_inv); } static double diff4_dzz(const double *arr, ptrdiff_t stride, double h_inv) { return FD24(arr, stride, h_inv); } static double diff4_dxz(const double *arr, ptrdiff_t stride, double dx_inv, double dz_inv) { return (-1.0 * (-1.0 * arr[ 2 * stride + 2] + 8.0 * arr[ 2 * stride + 1] - 8.0 * arr[ 2 * stride - 1] + 1.0 * arr[ 2 * stride - 2]) +8.0 * (-1.0 * arr[ 1 * stride + 2] + 8.0 * arr[ 1 * stride + 1] - 8.0 * arr[ 1 * stride - 1] + 1.0 * arr[ 1 * stride - 2]) -8.0 * (-1.0 * arr[-1 * stride + 2] + 8.0 * arr[-1 * stride + 1] - 8.0 * arr[-1 * stride - 1] + 1.0 * arr[-1 * stride - 2]) +1.0 * (-1.0 * arr[-2 * stride + 2] + 8.0 * arr[-2 * stride + 1] - 8.0 * arr[-2 * stride - 1] + 1.0 * arr[-2 * stride - 2])) * dx_inv * dz_inv / 144.0; } static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solver) { cGH *gh = ctx->gh; int ret; const ptrdiff_t stride_z = CCTK_GFINDEX3D(gh, 0, 0, 1); double *a_x = CCTK_VarDataPtr(gh, 0, "grid::x"); double *a_z = CCTK_VarDataPtr(gh, 0, "grid::z"); const double dx = a_x[1] - a_x[0]; const double dz = a_z[stride_z] - a_z[0]; const double dx_inv = 1.0 / dx; const double dz_inv = 1.0 / dz; double *a_gtxx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt11"); double *a_gtyy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt22"); double *a_gtzz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt33"); double *a_gtxy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt12"); double *a_gtxz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt13"); double *a_gtyz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt23"); double *a_Atxx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At11"); double *a_Atyy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At22"); double *a_Atzz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At33"); double *a_Atxy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At12"); double *a_Atxz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At13"); double *a_Atyz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At23"); double *a_trK = CCTK_VarDataPtr(gh, 0, "ML_BSSN::trK"); double *a_phi = CCTK_VarDataPtr(gh, 0, "ML_BSSN::phi"); double *a_Xtx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt1"); double *a_Xtz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt3"); double *a_kdotxx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot11"); double *a_kdotyy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot22"); double *a_kdotzz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot33"); double *a_kdotxy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot12"); double *a_kdotxz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot13"); double *a_kdotyz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot23"); double *a_alpha = CCTK_VarDataPtr(gh, 0, "ML_BSSN::alpha"); double *a_betax = CCTK_VarDataPtr(gh, 0, "ML_BSSN::beta1"); double *a_betaz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::beta3"); solver->diff_coeffs[MG2D_DIFF_COEFF_20]->boundaries[MG2D_BOUNDARY_0L].flags = MG2D_DC_FLAG_DISCONT; solver->diff_coeffs[MG2D_DIFF_COEFF_10]->boundaries[MG2D_BOUNDARY_0L].flags = MG2D_DC_FLAG_POLE; #pragma omp parallel for for (int idx_z = 0; idx_z < solver->local_size[1]; idx_z++) for (int idx_x = 0; idx_x < solver->local_size[0]; idx_x++) { const int idx_src = CCTK_GFINDEX3D(gh, idx_x + cp->offset_left[0], cp->y_idx, idx_z + cp->offset_left[1]); const int idx_dc = idx_z * solver->diff_coeffs[0]->stride + idx_x; const int idx_rhs = idx_z * solver->rhs_stride + idx_x; const double x = a_x[idx_src]; const int on_axis = fabs(x) < 1e-13; double K[3][3], Km[3][3], Ku[3][3], dK[3][3][3]; double gtu[3][3], g[3][3], gu[3][3]; double dg[3][3][3], dgu[3][3][3], d2g[3][3][3][3], gudot[3][3]; double G[3][3][3], dG[3][3][3][3], Gdot[3][3][3], X[3]; double Ric[3][3], Ricm[3][3]; double dgdot[3][3][3]; double k2, Kij_Dij_alpha, k_kdot, k3; double Gammadot_term, beta_term; double betal[3], dbetal[3][3], d2betal[3][3][3]; const double gtxx = a_gtxx[idx_src]; const double gtyy = a_gtyy[idx_src]; const double gtzz = a_gtzz[idx_src]; const double gtxy = a_gtxy[idx_src]; const double gtxz = a_gtxz[idx_src]; const double gtyz = a_gtyz[idx_src]; const double gt[3][3] = {{ gtxx, gtxy, gtxz }, { gtxy, gtyy, gtyz }, { gtxz, gtyz, gtzz }}; const double dx_gt11 = diff4_dx(&a_gtxx[idx_src], 1, dx_inv); const double dx_gt22 = diff4_dx(&a_gtyy[idx_src], 1, dx_inv); const double dx_gt33 = diff4_dx(&a_gtzz[idx_src], 1, dx_inv); const double dx_gt13 = diff4_dx(&a_gtxz[idx_src], 1, dx_inv); const double dz_gt11 = diff4_dz(&a_gtxx[idx_src], gh->cctk_lsh[0], dz_inv); const double dz_gt22 = diff4_dz(&a_gtyy[idx_src], gh->cctk_lsh[0], dz_inv); const double dz_gt33 = diff4_dz(&a_gtzz[idx_src], gh->cctk_lsh[0], dz_inv); const double dz_gt13 = diff4_dz(&a_gtxz[idx_src], gh->cctk_lsh[0], dz_inv); const double dgt[3][3][3] = { { { dx_gt11, 0.0, dx_gt13 }, { 0.0, dx_gt22, 0.0 }, { dx_gt13, 0.0, dx_gt33 }, }, { { 0.0, on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0 }, { on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0, on_axis ? dx_gt13 : gtxz / x }, { 0.0, on_axis ? dx_gt13 : gtxz / x, 0.0 }, }, { { dz_gt11, 0.0, dz_gt13 }, { 0.0, dz_gt22, 0.0 }, { dz_gt13, 0.0, dz_gt33 }, }, }; const double dxx_gt11 = diff4_dxx(&a_gtxx[idx_src], 1, dx_inv); const double dxx_gt22 = diff4_dxx(&a_gtyy[idx_src], 1, dx_inv); const double dxx_gt33 = diff4_dxx(&a_gtzz[idx_src], 1, dx_inv); const double dxx_gt13 = diff4_dxx(&a_gtxz[idx_src], 1, dx_inv); const double dxz_gt11 = diff4_dxz(&a_gtxx[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); const double dxz_gt22 = diff4_dxz(&a_gtyy[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); const double dxz_gt33 = diff4_dxz(&a_gtzz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); const double dxz_gt13 = diff4_dxz(&a_gtxz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); const double dzz_gt11 = diff4_dzz(&a_gtxx[idx_src], gh->cctk_lsh[0], dz_inv); const double dzz_gt22 = diff4_dzz(&a_gtyy[idx_src], gh->cctk_lsh[0], dz_inv); const double dzz_gt33 = diff4_dzz(&a_gtzz[idx_src], gh->cctk_lsh[0], dz_inv); const double dzz_gt13 = diff4_dzz(&a_gtxz[idx_src], gh->cctk_lsh[0], dz_inv); const double d2gt[3][3][3][3] = { { { { dxx_gt11, 0.0, dxx_gt13 }, { 0.0, dxx_gt22, 0.0 }, { dxx_gt13, 0.0, dxx_gt33 }, }, { { 0.0, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 }, { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, { 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 }, }, { { dxz_gt11, 0.0, dxz_gt13 }, { 0.0, dxz_gt22, 0.0 }, { dxz_gt13, 0.0, dxz_gt33 }, }, }, { { { 0.0, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 }, { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, { 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 }, }, { { on_axis ? dxx_gt22 : dx_gt11 / x - 2 * (gtxx - gtyy) / SQR(x), 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) }, { 0.0, on_axis ? dxx_gt11 : dx_gt22 / x + 2.0 * (gtxx - gtyy) / SQR(x), 0.0 }, { on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0, on_axis ? dxx_gt33 : dx_gt33 / x }, }, { { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 }, { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0, on_axis ? dxz_gt13 : dz_gt13 / x }, { 0.0, on_axis ? dxz_gt13 : dz_gt13 / x, 0.0 }, }, }, { { { dxz_gt11, 0.0, dxz_gt13 }, { 0.0, dxz_gt22, 0.0 }, { dxz_gt13, 0.0, dxz_gt33 }, }, { { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 }, { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0, on_axis ? dxz_gt13 : dz_gt13 / x }, { 0.0, on_axis ? dxz_gt13 : dz_gt13 / x, 0.0 }, }, { { dzz_gt11, 0.0, dzz_gt13 }, { 0.0, dzz_gt22, 0.0 }, { dzz_gt13, 0.0, dzz_gt33 }, }, }, }; const double Atxx = a_Atxx[idx_src]; const double Atyy = a_Atyy[idx_src]; const double Atzz = a_Atzz[idx_src]; const double Atxy = a_Atxy[idx_src]; const double Atxz = a_Atxz[idx_src]; const double Atyz = a_Atyz[idx_src]; const double At[3][3] = { { Atxx, Atxy, Atxz }, { Atxy, Atyy, Atyz }, { Atxz, Atyz, Atzz }}; const double dx_At11 = diff4_dx(&a_Atxx[idx_src], 1, dx_inv); const double dx_At22 = diff4_dx(&a_Atyy[idx_src], 1, dx_inv); const double dx_At33 = diff4_dx(&a_Atzz[idx_src], 1, dx_inv); const double dx_At13 = diff4_dx(&a_Atxz[idx_src], 1, dx_inv); const double dz_At11 = diff4_dz(&a_Atxx[idx_src], gh->cctk_lsh[0], dz_inv); const double dz_At22 = diff4_dz(&a_Atyy[idx_src], gh->cctk_lsh[0], dz_inv); const double dz_At33 = diff4_dz(&a_Atzz[idx_src], gh->cctk_lsh[0], dz_inv); const double dz_At13 = diff4_dz(&a_Atxz[idx_src], gh->cctk_lsh[0], dz_inv); const double dAt[3][3][3] = { { { dx_At11, 0.0, dx_At13 }, { 0.0, dx_At22, 0.0 }, { dx_At13, 0.0, dx_At33 }, }, { { 0.0, on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0 }, { on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0, on_axis ? dx_At13 : Atxz / x }, { 0.0, on_axis ? dx_At13 : Atxz / x, 0.0 }, }, { { dz_At11, 0.0, dz_At13 }, { 0.0, dz_At22, 0.0 }, { dz_At13, 0.0, dz_At33 }, }, }; const double phi = a_phi[idx_src]; const double dx_phi = diff4_dx(&a_phi[idx_src], 1, dx_inv); const double dz_phi = diff4_dz(&a_phi[idx_src], gh->cctk_lsh[0], dz_inv); const double dphi[3] = { dx_phi, 0.0, dz_phi }; const double dxx_phi = diff4_dxx(&a_phi[idx_src], gh->cctk_lsh[0], dx_inv); const double dzz_phi = diff4_dzz(&a_phi[idx_src], gh->cctk_lsh[0], dz_inv); const double dxz_phi = diff4_dxz(&a_phi[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); const double d2phi[3][3] = { { dxx_phi, 0.0, dxz_phi }, { 0.0, on_axis ? dxx_phi : dx_phi / x, 0.0 }, { dxz_phi, 0.0, dzz_phi }, }; const double trK = a_trK[idx_src]; const double dx_trk = diff4_dx(&a_trK[idx_src], 1, dx_inv); const double dz_trk = diff4_dz(&a_trK[idx_src], gh->cctk_lsh[0], dz_inv); const double dtrK[3] = { dx_trk, 0.0, dz_trk }; const double kdot_xx = a_kdotxx[idx_src]; const double kdot_yy = a_kdotyy[idx_src]; const double kdot_zz = a_kdotzz[idx_src]; const double kdot_xy = a_kdotxy[idx_src]; const double kdot_xz = a_kdotxz[idx_src]; const double kdot_yz = a_kdotyz[idx_src]; const double kdot[3][3] = {{ kdot_xx, kdot_xy, kdot_xz }, { kdot_xy, kdot_yy, kdot_yz }, { kdot_xz, kdot_yz, kdot_zz }}; const double Xtx = a_Xtx[idx_src]; const double Xtz = a_Xtz[idx_src]; const double alpha = a_alpha[idx_src]; const double dx_alpha = diff4_dx(&a_alpha[idx_src], 1, dx_inv); const double dz_alpha = diff4_dz(&a_alpha[idx_src], gh->cctk_lsh[0], dz_inv); const double dalpha[3] = { dx_alpha, 0.0, dz_alpha }; const double dxx_alpha = diff4_dxx(&a_alpha[idx_src], 1, dx_inv); const double dzz_alpha = diff4_dzz(&a_alpha[idx_src], gh->cctk_lsh[0], dz_inv); const double dxz_alpha = diff4_dxz(&a_alpha[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); const double d2alpha[3][3] = {{ dxx_alpha, 0, dxz_alpha }, { 0, on_axis ? dxx_alpha : dx_alpha / x, 0 }, { dxz_alpha, 0, dzz_alpha }}; const double betax = a_betax[idx_src]; const double betaz = a_betaz[idx_src]; const double beta[3] = { betax, 0.0, betaz }; const double dx_betax = diff4_dx(&a_betax[idx_src], 1, dx_inv); const double dz_betax = diff4_dz(&a_betax[idx_src], gh->cctk_lsh[0], dz_inv); const double dx_betaz = diff4_dx(&a_betaz[idx_src], 1, dx_inv); const double dz_betaz = diff4_dz(&a_betaz[idx_src], gh->cctk_lsh[0], dz_inv); const double dbeta[3][3] = {{ dx_betax, 0.0, dx_betaz }, { 0.0, on_axis ? dx_betax : betax / x, 0.0 }, { dz_betax, 0.0, dz_betaz }}; const double dxx_betax = diff4_dxx(&a_betax[idx_src], gh->cctk_lsh[0], dx_inv); const double dxz_betax = diff4_dxz(&a_betax[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); const double dzz_betax = diff4_dzz(&a_betax[idx_src], gh->cctk_lsh[0], dz_inv); const double dxx_betaz = diff4_dxx(&a_betaz[idx_src], gh->cctk_lsh[0], dx_inv); const double dxz_betaz = diff4_dxz(&a_betaz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); const double dzz_betaz = diff4_dzz(&a_betaz[idx_src], gh->cctk_lsh[0], dz_inv); const double d2beta[3][3][3] = { { { dxx_betax, 0.0, dxx_betaz }, { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 }, { dxz_betax, 0.0, dxz_betaz }, }, { { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 }, { on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0, on_axis ? dxx_betaz : dx_betaz / x }, { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 }, }, { { dxz_betax, 0.0, dxz_betaz }, { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 }, { dzz_betax, 0.0, dzz_betaz }, }, }; const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz); // \tilde{γ}^{ij} gtu[0][0] = (gtyy * gtzz - SQR(gtyz)) / det; gtu[1][1] = (gtxx * gtzz - SQR(gtxz)) / det; gtu[2][2] = (gtxx * gtyy - SQR(gtxy)) / det; gtu[0][1] = -(gtxy * gtzz - gtyz * gtxz) / det; gtu[0][2] = (gtxy * gtyz - gtyy * gtxz) / det; gtu[1][2] = -(gtxx * gtyz - gtxy * gtxz) / det; gtu[1][0] = gtu[0][1]; gtu[2][0] = gtu[0][2]; gtu[2][1] = gtu[1][2]; //　γ_{jk}/^{jk} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { gu[j][k] = SQR(phi) * gtu[j][k]; g[j][k] = gt[j][k] / SQR(phi); } // β_j for (int j = 0; j < 3; j++) { double val = 0.0; for (int k = 0; k < 3; k++) val += g[j][k] * beta[k]; betal[j] = val; } // ∂_j γ_{kl} // ∂_j K_{kl} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) { dg[j][k][l] = -2.0 * dphi[j] * gt[k][l] / (phi * SQR(phi)) + dgt[j][k][l] / SQR(phi); dK[j][k][l] = -2.0 * dphi[j] * At[k][l] / (phi * SQR(phi)) + dAt[j][k][l] / SQR(phi) + (1.0 / 3.0) * (dtrK[j] * g[k][l] + trK * dg[j][k][l]); } // ∂_{jk} g_{lm} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) for (int m = 0; m < 3; m++) { d2g[j][k][l][m] = 6.0 * gt [l][m] * dphi[j] * dphi[k] / SQR(SQR(phi)) - 2.0 * gt [l][m] * d2phi[j][k] / (phi * SQR(phi)) - 2.0 * dgt [j][l][m] * dphi[k] / (phi * SQR(phi)) - 2.0 * dgt [k][l][m] * dphi[j] / (phi * SQR(phi)) + d2gt[j][k][l][m] / SQR(phi); } // ∂_j γ^{kl} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) { double val = 0.0; for (int m = 0; m < 3; m++) for (int n = 0; n < 3; n++) val += -gu[k][m] * gu[l][n] * dg[j][m][n]; dgu[j][k][l] = val; } // Γ^j_{kl} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) { double val = 0.0; for (int m = 0; m < 3; m++) val += 0.5 * gu[j][m] * (dg[k][l][m] + dg[l][k][m] - dg[m][k][l]); G[j][k][l] = val; } // ∂_j Γ^k_{lm} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) for (int m = 0; m < 3; m++) { double val = 0.0; for (int n = 0; n < 3; n++) { val += dgu[j][k][n] * (dg [l][m][n] + dg [m][l][n] - dg [n][l][m]) + gu [k][n] * (d2g[j][l][m][n] + d2g[j][m][l][n] - d2g[j][n][l][m]); } dG[j][k][l][m] = 0.5 * val; } // Γ^j = γ^{kl} Γ_{kl}^j for (int j = 0; j < 3; j++) { double val = 0.0; for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) val += gu[k][l] * G[j][k][l]; X[j] = val; } // ∂_j β_k for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) val += g[k][l] * dbeta[j][l] + dg[j][k][l] * beta[l]; dbetal[j][k] = val; } // ∂_{jk} β_l for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) { double val = 0.0; for (int m = 0; m < 3; m++) val += d2g[j][k][l][m] * beta[m] + dg[k][l][m] * dbeta[j][m] + dg[j][l][m] * dbeta[k][m] + g[l][m] * d2beta[j][k][m]; d2betal[j][k][l] = val; } // K_{ij} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) K[j][k] = At[j][k] / SQR(phi) + (1.0 / 3.0) * g[j][k] * trK; for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) val += gu[j][l] * K[l][k]; Km[j][k] = val; } // K^{ij} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) val += gu[j][l] * Km[k][l]; Ku[j][k] = val; } // ∂_j \dot{γ_kl} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) { dgdot[j][k][l] = -2.0 * (dalpha[j] * K[k][l] + alpha * dK[j][k][l]) + d2betal[j][k][l] + d2betal[j][l][k]; for (int m = 0; m < 3; m++) dgdot[j][k][l] += -2.0 * (dG[j][m][k][l] * betal[m] + G[m][k][l] * dbetal[j][m]); } // \dot{γ^{jk}} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { gudot[j][k] = 2.0 * alpha * Ku[j][k]; for (int l = 0; l < 3; l++) { gudot[j][k] += -gu[j][l] * dbeta[l][k] - gu[k][l] * dbeta[l][j]; for (int m = 0; m < 3; m++) gudot[j][k] += -gu[j][l] * G[k][l][m] * beta[m] - gu[k][l] * G[j][l][m] * beta[m]; } } // \dot{Γ}^j_{kl} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) { double val = 0.0; for (int m = 0; m < 3; m++) val += gudot[j][m] * (dg [k][l][m] + dg [l][k][m] - dg [m][k][l]) + gu [j][m] * (dgdot[k][l][m] + dgdot[l][k][m] - dgdot[m][k][l]); Gdot[j][k][l] = 0.5 * val; } Gammadot_term = 0.0; for (int j = 0; j < 3; j++) { double val = 0.0; for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) { val += gu[k][l] * Gdot[j][k][l]; } Gammadot_term += val * dalpha[j]; } Kij_Dij_alpha = 0.0; for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) val += G[l][j][k] * dalpha[l]; Kij_Dij_alpha += Ku[j][k] * (d2alpha[j][k] - val); } k_kdot = 0.0; for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) k_kdot += kdot[j][k] * Ku[j][k]; k3 = 0.0; for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) val += Km[k][l] * Km[l][j]; k3 += val * Km[j][k]; } // K_{ij} K^{ij} k2 = 0.0; for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) k2 += Km[j][k] * Km[k][j]; // Ric_{jk} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { double val = 0.0; for (int m = 0; m < 3; m++) val += dG[m][m][j][k] - dG[k][m][j][m]; for (int m = 0; m < 3; m++) for (int l = 0; l < 3; l++) val += G[l][l][m] * G[m][j][k] - G[l][k][m] * G[m][j][l]; Ric[j][k] = val; } // Ric^j_k for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) val += gu[j][l] * Ric[l][k]; Ricm[j][k] = val; } beta_term = 0.0; for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) { double D2beta = 0.0; D2beta += d2beta[j][k][l]; for (int m = 0; m < 3; m++) { D2beta += dG[j][l][k][m] * beta[m] + G[l][k][m] * dbeta[j][m] + G[l][j][m] * dbeta[k][m] - G[m][j][k] * dbeta[m][l]; for (int n = 0; n < 3; n++) D2beta += G[l][j][m] * G[m][k][n] * beta[n] - G[m][j][k] * G[l][m][n] * beta[n]; } beta_term += -gu[j][k] * D2beta * dalpha[l]; } for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) beta_term += -beta[k] * Ricm[j][k] * dalpha[j]; solver->diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_dc] = gu[0][0] + (on_axis ? gu[1][1] : 0.0); solver->diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_dc] = gu[2][2]; solver->diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_dc] = 2.0 * gu[0][2]; solver->diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_dc] = -X[0] + (on_axis ? 0.0 : gu[1][1] / x); solver->diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_dc] = -X[2]; solver->diff_coeffs[MG2D_DIFF_COEFF_00]->data[idx_dc] = -k2; solver->rhs[idx_rhs] = -2.0 * alpha * Kij_Dij_alpha + Gammadot_term + 2.0 * alpha * (k_kdot + 2.0 * alpha * k3) + beta_term; if (on_axis) { solver->diff_coeffs[MG2D_DIFF_COEFF_20]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1]; solver->diff_coeffs[MG2D_DIFF_COEFF_10]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1]; } } } static void solution_to_grid(CoordPatch *cp, MG2DContext *solver, double *dst) { #pragma omp parallel for for (int j = 0; j < solver->local_size[1]; j++) for (int i = 0; i < solver->local_size[0]; 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 * solver->u_stride + i; dst[idx_dst] = solver->u[idx_src]; } /* fill in the axial symmetry ghostpoints by mirroring */ if (!cp->bnd_intercomp[0][0]) { #pragma omp parallel for 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]; } } if (!cp->bnd_intercomp[1][0]) { #pragma omp parallel for 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]; } } } void qms_mg_solve(CCTK_ARGUMENTS) { QMSMGContext *ms = qms_context; CoordPatch *cp; DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; const int reflevel = ctz(cctkGH->cctk_levfac[0]); const double time = ms->gh->cctk_time; int64_t grid_size; int ret; if (cctkGH->cctk_time >= switchoff_time + 1.0) return; if (reflevel < ms->solve_level) return; if (reflevel > ms->solve_level && fabs(time - W_val1_time[ms->solve_level]) > 1e-13) return; cp = get_coord_patch(ms, reflevel); grid_size = cctkGH->cctk_lsh[0] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[2]; fill_eq_coeffs(ms, cp, cp->solver); if (reflevel > ms->solve_level) { /* fill the solver boundary conditions */ 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++) { 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->local_size[0] + j; i++) { const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, ABS(i) + cp->offset_left[0], cp->y_idx, cp->offset_left[1] + cp->solver->local_size[1] - 1 + j); dst[i] = W_mg[idx]; } } } if (cp->solver->local_start[0] + cp->solver->local_size[0] == cp->solver->domain_size) { #pragma omp parallel for 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->local_size[1] + j; i++) { const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cp->offset_left[1] + cp->solver->local_size[0] - 1 + j, cp->y_idx, cp->offset_left[1] + ABS(i)); dst[i] = W_mg[idx]; } } } { /* use the solution from the coarser level as the initial guess */ CoordPatch *cp_coarse = get_coord_patch(ms, reflevel - 1); ret = mg2d_init_guess(cp->solver, cp_coarse->solver->u, cp_coarse->solver->u_stride, cp_coarse->solver->local_start, cp_coarse->solver->local_size, cp_coarse->solver->step); if (ret < 0) CCTK_WARN(0, "Error setting the initial guess"); } } fprintf(stderr, "%d qms solve: time %g\n", reflevel, time); ret = mg2d_solve(cp->solver); if (ret < 0) CCTK_WARN(0, "Error solving the quasi-maximal slicing equation"); { double *W_mg_1 = CCTK_VarDataPtr(cctkGH, 1, "QuasiMaximalSlicingMG::W_mg"); solution_to_grid(cp, cp->solver, W_mg); memcpy(W_mg_1, W_mg, grid_size * sizeof(*W_mg_1)); } /* add the solution to the list of past solutions */ { memcpy(W_val0, W_val1, sizeof(*W_val0) * grid_size); W_val0_time[reflevel] = W_val1_time[reflevel]; memcpy(W_val1, W_mg, sizeof(*W_val1) * grid_size); W_val1_time[reflevel] = ms->gh->cctk_time; } /* linearly extrapolate the past two solution to predict the next one */ { double *u0 = W_val0; double *u1 = W_val1; double time0 = W_val0_time[reflevel]; double time1 = W_val1_time[reflevel]; double time = time1 + ms->gh->cctk_delta_time / (1.0 * (1 << ms->solve_level)); double fact0, fact1; if (time0 == DBL_MAX && time1 == DBL_MAX) { fact0 = 0.0; fact1 = 0.0; } else if (time0 == DBL_MAX) { fact0 = 0.0; fact1 = 0.0; } else { fact0 = (time - time1) / (time0 - time1); fact1 = (time - time0) / (time1 - time0); } memcpy(W_pred0, W_pred1, sizeof(*W_pred0) * grid_size); W_pred0_time[reflevel] = W_pred1_time[reflevel]; for (int i = 0; i < grid_size; i++) W_pred1[i] = (u1[i] * fact1 + u0[i] * fact0); W_pred1_time[reflevel] = time; } } void qms_mg_sync(CCTK_ARGUMENTS) { } void qms_mg_eval(CCTK_ARGUMENTS) { QMSMGContext *ms; DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; const int reflevel = ctz(cctkGH->cctk_levfac[0]); const double time = cctkGH->cctk_time; const int64_t grid_size = cctkGH->cctk_lsh[0] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[2]; int i, ret; ms = qms_context; if (reflevel < ms->solve_level || time >= switchoff_time + 1.0) return; { double *u0 = W_pred0; double *u1 = W_pred1; double time0 = W_pred0_time[reflevel]; double time1 = W_pred1_time[reflevel]; double fact, fact0, fact1; if (time0 == DBL_MAX && time1 == DBL_MAX) { fact0 = 0.0; fact1 = 0.0; } else if (time0 == DBL_MAX) { fact0 = 0.0; fact1 = 1.0; } else { fact0 = (time - time1) / (time0 - time1); fact1 = (time - time0) / (time1 - time0); } if (time < switchoff_time) fact = 1.0; else fact = exp(-36.0 * pow((time - switchoff_time), 4.0)); #pragma omp parallel for for (int i = 0; i < grid_size; i++) W[i] = (u1[i] * fact1 + u0[i] * fact0) * fact; } } static int context_init(cGH *gh, int solve_level, int fd_stencil, int maxiter, int exact_size, int nb_cycles, int nb_relax_pre, int nb_relax_post, double tol_residual_base, double cfl_factor, const char *loglevel_str, QMSMGContext **ctx) { QMSMGContext *qms; int loglevel = -1; int ret; qms = calloc(1, sizeof(*qms)); if (!qms) return -ENOMEM; qms->gh = gh; qms->solve_level = solve_level; qms->fd_stencil = fd_stencil; qms->maxiter = maxiter; qms->max_exact_size = exact_size; qms->nb_cycles = nb_cycles; qms->nb_relax_pre = nb_relax_pre; qms->nb_relax_post = nb_relax_post; qms->tol_residual_base = tol_residual_base; qms->cfl_factor = cfl_factor; for (int i = 0; i < ARRAY_ELEMS(log_levels); i++) { if (!strcmp(loglevel_str, log_levels[i].str)) { loglevel = log_levels[i].level; break; } } if (loglevel < 0) { fprintf(stderr, "Invalid loglevel: %s\n", loglevel_str); return -EINVAL; } qms->log_level = loglevel; *ctx = qms; return 0; } void qms_mg_init(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; int use_tapered_grids = *(int*)CCTK_ParameterGet("use_tapered_grids", "Carpet", &use_tapered_grids); int ret; fprintf(stderr, "%d init\n", ctz(cctk_levfac[0])); if (!use_tapered_grids) CCTK_WARN(0, "MaximalSlicingAxiMG only works with use_tapered_grids=1"); if (!qms_context) { ret = context_init(cctkGH, solve_level, fd_stencil, maxiter, exact_size, nb_cycles, nb_relax_pre, nb_relax_post, tol_residual_base, cfl_factor, loglevel, &qms_context); if (ret < 0) CCTK_WARN(0, "Error initializing the solver context"); } } void qms_mg_inithist(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; size_t grid_size = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]; fprintf(stderr, "%d inithist\n", ctz(cctk_levfac[0])); for (int i = 0; i < 32; i++) { W_pred0_time[i] = DBL_MAX; W_pred1_time[i] = DBL_MAX; W_val0_time[i] = DBL_MAX; W_val1_time[i] = DBL_MAX; } memset(Kdot11, 0, sizeof(double) * grid_size); memset(Kdot22, 0, sizeof(double) * grid_size); memset(Kdot33, 0, sizeof(double) * grid_size); memset(Kdot12, 0, sizeof(double) * grid_size); memset(Kdot13, 0, sizeof(double) * grid_size); memset(Kdot23, 0, sizeof(double) * grid_size); memset(Xtdot1, 0, sizeof(double) * grid_size); memset(Xtdot2, 0, sizeof(double) * grid_size); memset(Xtdot3, 0, sizeof(double) * grid_size); memset(W_pred0, 0, sizeof(double) * grid_size); memset(W_pred1, 0, sizeof(double) * grid_size); memset(W_val0, 0, sizeof(double) * grid_size); memset(W_val1, 0, sizeof(double) * grid_size); }