simplify coordpatch initialization

1 files 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);