From 1fc5d6f33504ba348c389ac33d1036826511e906 Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Sun, 29 Apr 2018 12:38:49 +0200
Subject: Add support for setting maximal lapse in initial data only.

---
 param.ccl                 |   5 +++
 schedule.ccl              |   5 +++
 src/maximal_slicing_axi.c | 100 ++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 110 insertions(+)

diff --git a/param.ccl b/param.ccl
index 5309214..ab3bc38 100644
--- a/param.ccl
+++ b/param.ccl
@@ -7,6 +7,11 @@ EXTENDS KEYWORD lapse_evolution_method
   "maximal_axi" :: "Maximal slicing for an axisymmetric spacetime"
 }
 
+EXTENDS KEYWORD initial_lapse
+{
+  "maximal_axi" :: "Maximal slicing for an axisymmetric spacetime"
+}
+
 RESTRICTED:
 CCTK_REAL amplitude "Wave amplitude A."
 {
diff --git a/schedule.ccl b/schedule.ccl
index eafabd1..5641401 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -27,3 +27,8 @@ if (CCTK_Equals(lapse_evolution_method, "maximal_axi")) {
     }
 }
 
+if (CCTK_Equals(initial_lapse, "maximal_axi")) {
+    SCHEDULE maximal_slicing_axi_initlapse IN CCTK_POSTPOSTINITIAL {
+        LANG: C
+    } "Maximal slicing in axisymmetry"
+}
diff --git a/src/maximal_slicing_axi.c b/src/maximal_slicing_axi.c
index f1758cb..9689ec0 100644
--- a/src/maximal_slicing_axi.c
+++ b/src/maximal_slicing_axi.c
@@ -60,6 +60,17 @@ static double calc_basis_freq(const BasisSet *b, int order)
     return b->colloc_point(order, 1);
 }
 
+static void coord_patch_free(CoordPatch *cp)
+{
+    free(cp->basis_val_r);
+    free(cp->basis_val_z);
+    free(cp->transform_z);
+    free(cp->transform_matrix);
+    free(cp->transform_matrix1);
+    free(cp->transform_tmp);
+    free(cp->one);
+}
+
 static CoordPatch *get_coord_patch(MaximalSlicingContext *ms,
                                    CCTK_REAL *x, CCTK_REAL *y, CCTK_REAL *z)
 {
@@ -205,6 +216,95 @@ static int context_init(cGH *cctkGH, MaximalSlicingContext **ctx)
     return 0;
 }
 
+static void context_free(MaximalSlicingContext **pms)
+{
+    MaximalSlicingContext *ms = *pms;
+
+    if (!ms)
+        return;
+
+    for (int i = 0; i < ms->nb_patches; i++)
+        coord_patch_free(&ms->patches[i]);
+    free(ms->patches);
+
+    free(ms);
+    *pms = NULL;
+}
+
+void maximal_slicing_axi_initlapse(CCTK_ARGUMENTS)
+{
+    MaximalSlicingContext *ms;
+
+    CoordPatch *cp;
+
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    double *coeffs = NULL;
+    int i, ret;
+
+    context_init(cctkGH, &ms);
+
+    cp = get_coord_patch(ms, x, y, z);
+
+    ms_solver_solve(ms->solver);
+
+    coeffs = ms->solver->coeffs;
+
+#if 0
+#pragma omp parallel for
+    for (int k = 0; k < cctk_lsh[2]; k++) {
+        for (int i = 0; i < cctk_lsh[0]; i++) {
+            int idx = CCTK_GFINDEX3D(cctkGH, i, cp->y_idx, k);
+            double xx = x[idx];
+            double zz = z[idx];
+            double r = sqrt(SQR(xx) + SQR(zz));
+            double phi = atan2(zz, xx);
+
+            double val = 1.0;
+
+            for (int l = 0; l < ms->nb_coeffs_z; l++) {
+                double tmp = 0.0;
+                for (int m = 0; m < ms->nb_coeffs_x; m++) {
+                    const int idx_coeff = l * ms->nb_coeffs_x + m;
+                    tmp += ms->coeffs[idx_coeff] * ms->basis->eval(r, m);
+                }
+                val += tmp * ms->basis1->eval(phi, l);
+            }
+
+            alp[idx] = val;
+        }
+    }
+#else
+    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
+                cctk_lsh[0] * cctk_lsh[2], ms->solver->nb_coeffs[1], ms->solver->nb_coeffs[0],
+                1.0, cp->transform_matrix, cctk_lsh[0] * cctk_lsh[2],
+                coeffs, ms->solver->nb_coeffs[0], 0.0, cp->transform_tmp, cctk_lsh[0] * cctk_lsh[2]);
+#pragma omp parallel for
+    for (int j = 0; j < cctk_lsh[2]; j++)
+        for (int i = 0; i < cctk_lsh[0]; i++) {
+            const int idx_grid = j * cctk_lsh[0] + i;
+            const double val = cblas_ddot(ms->solver->nb_coeffs[1], cp->transform_matrix1 + idx_grid * ms->solver->nb_coeffs[1], 1,
+                                          cp->transform_tmp + idx_grid, cctk_lsh[0] * cctk_lsh[2]);
+            alpha[CCTK_GFINDEX3D(cctkGH, i, cp->y_idx, j)] = 1.0 + val;
+        }
+
+    memcpy(alp, alpha, sizeof(*alp) * cctk_lsh[0] * cctk_lsh[2]);
+#endif
+    //memcpy(alp, cp->one, cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2] * sizeof(*alp));
+    //memset(alp, 0, cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2] * sizeof(*alp));
+    //cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
+    //            ms->nb_coeffs_x, cctk_lsh[2], ms->nb_coeffs_z, 1.0,
+    //            coeffs, ms->nb_coeffs_x, cp->basis_val_z, ms->nb_coeffs_z,
+    //            0.0, cp->transform_z, ms->nb_coeffs_x);
+    //cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
+    //            cctk_lsh[1] * cctk_lsh[0], cctk_lsh[2], ms->nb_coeffs_x, 1.0,
+    //            cp->basis_val_r, cctk_lsh[0] * cctk_lsh[1], cp->transform_z, ms->nb_coeffs_x,
+    //            1.0, alp, cctk_lsh[0] * cctk_lsh[1]);
+
+    context_free(&ms);
+}
+
 void maximal_slicing_axi(CCTK_ARGUMENTS)
 {
     static MaximalSlicingContext *ms;
-- 
cgit v1.2.3