From 4a623f1553f807291a7f49bdf7ca64ee8b44b47e Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 23 Aug 2018 12:42:39 +0200 Subject: simplify coordpatch initialization --- src/maximal_slicing_axi.c | 67 +++++++++++++++++++++++++---------------------- 1 file changed, 36 insertions(+), 31 deletions(-) diff --git a/src/maximal_slicing_axi.c b/src/maximal_slicing_axi.c index 9689ec0..4f4c72e 100644 --- a/src/maximal_slicing_axi.c +++ b/src/maximal_slicing_axi.c @@ -24,9 +24,7 @@ /* precomputed values for a given refined grid */ typedef struct CoordPatch { - CCTK_REAL origin[3]; - CCTK_INT delta[3]; - CCTK_INT size[3]; + int level; // basis values on the grid double *basis_val_r; @@ -54,6 +52,21 @@ typedef struct MaximalSlicingContext { double scale_factor; +static int ctz(int a) +{ + int ret = 0; + + if (!a) + return INT_MAX; + + while (!(a & 1)) { + a >>= 1; + ret++; + } + + return ret; +} + /* get an approximate "main" frequency component in a basis function */ static double calc_basis_freq(const BasisSet *b, int order) { @@ -71,10 +84,10 @@ static void coord_patch_free(CoordPatch *cp) free(cp->one); } -static CoordPatch *get_coord_patch(MaximalSlicingContext *ms, +static CoordPatch *get_coord_patch(MaximalSlicingContext *ms, int level, CCTK_REAL *x, CCTK_REAL *y, CCTK_REAL *z) { - cGH *cctkGH = ms->gh; + cGH *gh = ms->gh; int nb_coeffs_x = ms->solver->nb_coeffs[0]; int nb_coeffs_z = ms->solver->nb_coeffs[1]; @@ -85,19 +98,11 @@ static CoordPatch *get_coord_patch(MaximalSlicingContext *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; } - grid_size = cctkGH->cctk_lsh[0] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[2]; + grid_size = gh->cctk_lsh[0] * gh->cctk_lsh[1] * gh->cctk_lsh[2]; /* create a new patch */ ms->patches = realloc(ms->patches, sizeof(*ms->patches) * (ms->nb_patches + 1)); @@ -105,16 +110,14 @@ static CoordPatch *get_coord_patch(MaximalSlicingContext *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 = level; - for (i = 0; i < cp->size[1]; i++) - if (fabs(y[CCTK_GFINDEX3D(cctkGH, 0, i, 0)]) < 1e-8) { + for (i = 0; i < gh->cctk_lsh[1]; i++) + if (fabs(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"); #if 0 @@ -142,14 +145,14 @@ static CoordPatch *get_coord_patch(MaximalSlicingContext *ms, for (int i = 0; i < grid_size; i++) cp->one[i] = 1.0; #else - posix_memalign((void**)&cp->transform_matrix, 32, sizeof(*cp->transform_matrix) * nb_coeffs_x * cp->size[0] * cp->size[2]); - posix_memalign((void**)&cp->transform_matrix1, 32, sizeof(*cp->transform_matrix1) * nb_coeffs_z * cp->size[0] * cp->size[2]); + posix_memalign((void**)&cp->transform_matrix, 32, sizeof(*cp->transform_matrix) * nb_coeffs_x * gh->cctk_lsh[0] * gh->cctk_lsh[2]); + posix_memalign((void**)&cp->transform_matrix1, 32, sizeof(*cp->transform_matrix1) * nb_coeffs_z * gh->cctk_lsh[0] * gh->cctk_lsh[2]); #pragma omp parallel for - for (int j = 0; j < cp->size[2]; j++) { + for (int j = 0; j < gh->cctk_lsh[2]; j++) { CCTK_REAL zz = z[CCTK_GFINDEX3D(ms->gh, 0, 0, j)]; - for (int i = 0; i < cp->size[0]; i++) { - const int idx_grid = j * cp->size[0] + i; + for (int i = 0; i < gh->cctk_lsh[0]; i++) { + const int idx_grid = j * gh->cctk_lsh[0] + i; double xx = x[CCTK_GFINDEX3D(ms->gh, i, 0, 0)]; double rr = sqrt(SQR(xx) + SQR(zz)); @@ -165,7 +168,7 @@ static CoordPatch *get_coord_patch(MaximalSlicingContext *ms, //for (int k = 0; k < ms->nb_coeffs_z; k++) // for (int l = 0; l < ms->nb_coeffs_x; l++) { // const int idx_coeff = k * ms->nb_coeffs_x + l; - // cp->transform_matrix[idx_grid + cp->size[0] * cp->size[2] * idx_coeff] = ms->basis->eval(r, l) * ms->basis1->eval(phi, k); + // cp->transform_matrix[idx_grid + gh->cctk_lsh[0] * gh->cctk_lsh[2] * idx_coeff] = ms->basis->eval(r, l) * ms->basis1->eval(phi, k); // } for (int k = 0; k < nb_coeffs_x; k++) { double dx = calc_basis_freq(ms->solver->basis[0], k); @@ -173,7 +176,7 @@ static CoordPatch *get_coord_patch(MaximalSlicingContext *ms, double fact = exp(-36.0 * pow(rr / r0, 10)); fact = 1.0; - cp->transform_matrix[idx_grid + cp->size[0] * cp->size[2] * k] = ms->solver->basis[0]->eval(coord0, k) * fact; + cp->transform_matrix[idx_grid + gh->cctk_lsh[0] * gh->cctk_lsh[2] * k] = ms->solver->basis[0]->eval(coord0, k) * fact; } for (int k = 0; k < nb_coeffs_z; k++) { double dx = calc_basis_freq(ms->solver->basis[1], k); @@ -185,7 +188,7 @@ static CoordPatch *get_coord_patch(MaximalSlicingContext *ms, } } } - posix_memalign((void**)&cp->transform_tmp, 32, sizeof(*cp->transform_tmp) * cp->size[0] * cp->size[2] * nb_coeffs_z); + posix_memalign((void**)&cp->transform_tmp, 32, sizeof(*cp->transform_tmp) * gh->cctk_lsh[0] * gh->cctk_lsh[2] * nb_coeffs_z); #endif ms->nb_patches++; @@ -240,12 +243,13 @@ void maximal_slicing_axi_initlapse(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; + const int reflevel = ctz(cctkGH->cctk_levfac[0]); double *coeffs = NULL; int i, ret; context_init(cctkGH, &ms); - cp = get_coord_patch(ms, x, y, z); + cp = get_coord_patch(ms, reflevel, x, y, z); ms_solver_solve(ms->solver); @@ -314,6 +318,7 @@ void maximal_slicing_axi(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; + const int reflevel = ctz(cctkGH->cctk_levfac[0]); int64_t expand_start, totaltime_start; double *coeffs = NULL; @@ -325,7 +330,7 @@ void maximal_slicing_axi(CCTK_ARGUMENTS) if (!ms) context_init(cctkGH, &ms); - cp = get_coord_patch(ms, x, y, z); + cp = get_coord_patch(ms, reflevel, x, y, z); CCTK_TimerStart("MaximalSlicingAxi_Solve"); ms_solver_solve(ms->solver); -- cgit v1.2.3