Implement evaluating maximal lapse for the initial data.

4 files changed, 294 insertions, 5 deletions

diff --git a/init.c b/init.c
index 33e5587..e5ad3d8 100644
--- a/init.c
+++ b/init.c
@@ -53,6 +53,8 @@ typedef struct TDPriv {
     double *coeffs_tmp;
     ptrdiff_t coeffs_stride;
 
+    double *coeffs_lapse;
+
     int cpu_flags;
 
     double (*scalarproduct_metric)(size_t len1, size_t len2, double *mat,
@@ -500,12 +502,240 @@ void td_context_free(TDContext **ptd)
 
     free(s->coeffs);
     free(s->coeffs_tmp);
+    free(s->coeffs_lapse);
 
     free(td->priv);
     free(td);
     *ptd = NULL;
 }
 
+typedef struct MaximalSlicingEval {
+    double *psi;
+    double *dpsi_r;
+    double *dpsi_t;
+    double *krr;
+    double *kpp;
+    double *krt;
+} MaximalSlicingEval;
+
+static int maximal_slicing_eval(void *opaque, unsigned int eq_idx,
+                                const unsigned int *colloc_grid_order,
+                                const double * const *colloc_grid,
+                                const double * const (*vars)[PSSOLVE_DIFF_ORDER_NB],
+                                double *dst)
+{
+    MaximalSlicingEval *max = opaque;
+
+    for (int idx1 = 0; idx1 < colloc_grid_order[1]; idx1++) {
+        const double theta = colloc_grid[1][idx1];
+        const double st  = sin(theta);
+        const double st2 = SQR(st);
+        const double ct  = cos(theta);
+
+        for (int idx0 = 0; idx0 < colloc_grid_order[0]; idx0++) {
+            const double r  = colloc_grid[0][idx0];
+            const double r2 = SQR(r);
+
+            const int idx = idx1 * colloc_grid_order[0] + idx0;
+
+            const double alpha      = vars[0][PSSOLVE_DIFF_ORDER_00][idx] + 1.0;
+            const double dalpha_r   = vars[0][PSSOLVE_DIFF_ORDER_10][idx];
+            const double dalpha_t   = vars[0][PSSOLVE_DIFF_ORDER_01][idx];
+            const double d2alpha_rr = vars[0][PSSOLVE_DIFF_ORDER_20][idx];
+            const double d2alpha_tt = vars[0][PSSOLVE_DIFF_ORDER_02][idx];
+            const double d2alpha_rt = vars[0][PSSOLVE_DIFF_ORDER_11][idx];
+
+            const double d2alpha[3][3] = {
+                { d2alpha_rr, d2alpha_rt, 0.0 },
+                { d2alpha_rt, d2alpha_tt, 0.0 },
+                { 0.0, 0.0, 0.0 },
+            };
+
+            const double psi  = max->psi[idx];
+            const double psi2 = SQR(psi);
+            const double psi3 = psi * psi2;
+            const double psi4 = SQR(psi2);
+
+            const double dpsi_r = max->dpsi_r[idx];
+            const double dpsi_t = max->dpsi_t[idx];
+
+            const double g[3][3] = {
+                { psi4, 0.0, 0.0 },
+                { 0.0, psi4 * r2, 0.0 },
+                { 0.0, 0.0, psi4 * r2 * st2},
+            };
+            const double gu[3][3] = {
+                { 1.0 / psi4, 0.0, 0.0 },
+                { 0.0, 1.0 / (psi4 * r2), 0.0 },
+                { 0.0, 0.0, 1.0 / (psi4 * r2 * st2) },
+            };
+
+            const double dg[3][3][3] = {
+                {
+                    { 4.0 * dpsi_r * psi3, 0.0, 0.0 },
+                    { 0.0, 4.0 * dpsi_r * psi3 * r2 + 2.0 * r * psi4, 0.0 },
+                    { 0.0, 0.0, 4.0 * dpsi_r * psi3 * r2 * st2 + 2.0 * r * psi4 * st2 },
+                },
+                {
+                    { 4.0 * dpsi_t * psi3, 0.0, 0.0 },
+                    { 0.0, 4.0 * dpsi_t * psi3 * r2, 0.0 },
+                    { 0.0, 0.0, 4.0 * dpsi_t * psi3 * r2 * st2 + 2.0 * psi4 * r2 * st * ct },
+
+                },
+                {
+                    { 0.0, 0.0, 0.0 },
+                    { 0.0, 0.0, 0.0 },
+                    { 0.0, 0.0, 0.0 },
+                },
+            };
+
+            const double krr  = max->krr[idx];
+            const double kpp  = max->kpp[idx];
+            const double krt  = max->krt[idx];
+
+            const double Km[3][3] = {
+                { krr, krt, 0.0 },
+                { krt / r2, -(krr + kpp), 0.0 },
+                { 0.0, 0.0, kpp },
+            };
+
+            double G[3][3][3], X[3];
+            double laplace_alpha, k2;
+
+            for (int i = 0; i < 3; i++)
+                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[i][l] * (dg[k][j][l] + dg[j][k][l] - dg[l][j][k]);
+
+                        G[i][j][k] = 0.5 * val;
+                    }
+
+            for (int i = 0; i < 3; i++) {
+                double val = 0.0;
+                for (int j = 0; j < 3; j++)
+                    for (int k = 0; k < 3; k++)
+                        val += gu[j][k] * G[i][j][k];
+                X[i] = val;
+            }
+
+            laplace_alpha = 0.0;
+            for (int i = 0; i < 3; i++)
+                for (int j = 0; j < 3; j++)
+                    laplace_alpha += gu[i][j] * d2alpha[i][j];
+
+            laplace_alpha -= X[0] * dalpha_r + X[1] * dalpha_t;
+
+            k2 = 0.0;
+            for (int i = 0; i < 3; i++)
+                for (int j = 0; j < 3; j++)
+                    k2 += Km[i][j] * Km[j][i];
+
+            dst[idx] = laplace_alpha - alpha * k2; 
+        }
+    }
+
+    return 0;
+}
+
+static int lapse_solve_max(const TDContext *td)
+{
+    MaximalSlicingEval max = { NULL };
+    TDPriv *priv = td->priv;
+    NLSolveContext *nl;
+    double *coords_r = NULL;
+    double *coords_t = NULL;
+    double *coeffs;
+    int ret = 0;
+
+    ret = tdi_nlsolve_context_alloc(&nl, 1);
+
+    if (ret < 0) {
+        tdi_log(&priv->logger, 0, "Error allocating the non-linear solver\n");
+        return ret;
+    }
+
+    nl->logger  = priv->logger;
+    nl->tp      = priv->tp;
+    nl->maxiter = td->max_iter;
+    nl->atol    = td->atol;
+
+    nl->basis[0][0]       = priv->basis[0][0];
+    nl->basis[0][1]       = priv->basis[0][1];
+    nl->solve_order[0][0] = priv->basis_order[0][0];
+    nl->solve_order[0][1] = priv->basis_order[0][1];
+
+#if HAVE_OPENCL
+    nl->ocl_ctx   = priv->ocl_ctx;
+    nl->ocl_queue = priv->ocl_queue;
+#endif
+
+    ret = tdi_nlsolve_context_init(nl);
+    if (ret < 0) {
+        tdi_log(&priv->logger, 0, "Error initializing the non-linear solver\n");
+        goto fail;
+    }
+
+    ret  = posix_memalign((void**)&max.psi,  32, sizeof(*max.psi) * td->nb_coeffs[0] * td->nb_coeffs[1]);
+    ret |= posix_memalign((void**)&max.dpsi_r,  32, sizeof(*max.psi) * td->nb_coeffs[0] * td->nb_coeffs[1]);
+    ret |= posix_memalign((void**)&max.dpsi_t,  32, sizeof(*max.psi) * td->nb_coeffs[0] * td->nb_coeffs[1]);
+    ret |= posix_memalign((void**)&max.krr,  32, sizeof(*max.krr) * td->nb_coeffs[0] * td->nb_coeffs[1]);
+    ret |= posix_memalign((void**)&max.kpp,  32, sizeof(*max.kpp) * td->nb_coeffs[0] * td->nb_coeffs[1]);
+    ret |= posix_memalign((void**)&max.krt,  32, sizeof(*max.krt) * td->nb_coeffs[0] * td->nb_coeffs[1]);
+    ret |= posix_memalign((void**)&coords_r, 32, sizeof(*coords_r) * td->nb_coeffs[0] * td->nb_coeffs[1]);
+    ret |= posix_memalign((void**)&coords_t, 32, sizeof(*coords_t) * td->nb_coeffs[0] * td->nb_coeffs[1]);
+    ret |= posix_memalign((void**)&coeffs,   32, sizeof(*coeffs) * td->nb_coeffs[0] * td->nb_coeffs[1]);
+    if (ret) {
+        ret = -ENOMEM;
+        goto fail;
+    }
+
+    for (int j = 0; j < td->nb_coeffs[1]; j++) {
+        for (int i = 0; i < td->nb_coeffs[0]; i++) {
+            coords_r[j * td->nb_coeffs[0] + i] = nl->colloc_grid[0][0][i];
+            coords_t[j * td->nb_coeffs[0] + i] = nl->colloc_grid[0][1][j];
+        }
+    }
+
+    td_eval_psi(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t,
+                (const unsigned int [2]){ 0, 0 }, max.psi);
+    td_eval_psi(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t,
+                (const unsigned int [2]){ 1, 0 }, max.dpsi_r);
+    td_eval_psi(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t,
+                (const unsigned int [2]){ 0, 1 }, max.dpsi_t);
+    td_eval_krr(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t,
+                (const unsigned int [2]){ 0, 0 }, max.krr);
+    td_eval_kpp(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t,
+                (const unsigned int [2]){ 0, 0 }, max.kpp);
+    td_eval_krt(td, td->nb_coeffs[0] * td->nb_coeffs[1], coords_r, coords_t,
+                (const unsigned int [2]){ 0, 0 }, max.krt);
+
+    ret = tdi_nlsolve_solve(nl, maximal_slicing_eval, NULL, &max, coeffs, 0);
+    if (ret < 0) {
+        tdi_log(&priv->logger, 0, "Error solving the maximal slicing equation\n");
+        goto fail;
+    }
+
+    priv->coeffs_lapse = coeffs;
+    coeffs             = NULL;
+
+fail:
+    free(coords_r);
+    free(coords_t);
+    free(coeffs);
+    free(max.psi);
+    free(max.dpsi_r);
+    free(max.dpsi_t);
+    free(max.krr);
+    free(max.kpp);
+    free(max.krt);
+
+    tdi_nlsolve_context_free(&nl);
+    return ret;
+}
+
 static double scalarproduct_metric_c(size_t len1, size_t len2, double *mat,
                                      double *vec1, double *vec2)
 {
@@ -528,7 +758,7 @@ static int eval_var(const TDContext *td, unsigned int var_idx,
     double *basis_val[2] = { NULL };
     int ret = 0;
 
-    if (diff_order[0] > 0 || diff_order[1] > 0) {
+    if (diff_order[0] > 2 || diff_order[1] > 2) {
         tdi_log(&s->logger, 0, "Derivatives of higher order than 2 are not supported.\n");
         return -ENOSYS;
     }
@@ -546,12 +776,12 @@ static int eval_var(const TDContext *td, unsigned int var_idx,
         double theta_val = theta[i];
         double r_val     = r[i];
 
-        double val = (var_idx == 0) ? 1.0 : 0.0;
+        double val = (var_idx == 0 && diff_order[0] == 0 && diff_order[1] == 0) ? 1.0 : 0.0;
 
         for (int k = 0; k < td->nb_coeffs[0]; k++)
-            basis_val[0][k] = tdi_basis_eval(s->basis[var_idx][0], BS_EVAL_TYPE_VALUE, r_val, k);
+            basis_val[0][k] = tdi_basis_eval(s->basis[var_idx][0], diff_order[0], r_val, k);
         for (int k = 0; k < td->nb_coeffs[1]; k++)
-            basis_val[1][k] = tdi_basis_eval(s->basis[var_idx][1], BS_EVAL_TYPE_VALUE, theta_val, k);
+            basis_val[1][k] = tdi_basis_eval(s->basis[var_idx][1], diff_order[1], theta_val, k);
 
         val += scalarproduct_metric_c(td->nb_coeffs[0], td->nb_coeffs[1], td->coeffs[var_idx],
                                       basis_val[0], basis_val[1]);
@@ -621,3 +851,56 @@ int td_eval_krt(const TDContext *td,
 
     return 0;
 }
+
+int td_eval_lapse(const TDContext *td,
+                  size_t nb_coords, const double *r, const double *theta,
+                  const unsigned int diff_order[2],
+                  double *out)
+{
+    TDPriv *priv = td->priv;
+    double *basis_val[2] = { NULL };
+    int ret;
+
+    if (!priv->coeffs_lapse) {
+        ret = lapse_solve_max(td);
+        if (ret < 0)
+            return ret;
+    }
+
+    if (diff_order[0] > 2 || diff_order[1] > 2) {
+        tdi_log(&priv->logger, 0, "Derivatives of higher order than 2 are not supported.\n");
+        return -ENOSYS;
+    }
+
+    posix_memalign((void**)&basis_val[0], 32, sizeof(*basis_val[0]) * td->nb_coeffs[0]);
+    posix_memalign((void**)&basis_val[1], 32, sizeof(*basis_val[1]) * td->nb_coeffs[1]);
+    if (!basis_val[0] || !basis_val[1]) {
+        ret = -ENOMEM;
+        goto fail;
+    }
+    memset(basis_val[0], 0, sizeof(*basis_val[0]) * td->nb_coeffs[0]);
+    memset(basis_val[1], 0, sizeof(*basis_val[1]) * td->nb_coeffs[1]);
+
+    for (int i = 0; i < nb_coords; i++) {
+        double theta_val = theta[i];
+        double r_val     = r[i];
+
+        double val = (diff_order[0] == 0 && diff_order[1] == 0) ? 1.0 : 0.0;
+
+        for (int k = 0; k < td->nb_coeffs[0]; k++)
+            basis_val[0][k] = tdi_basis_eval(priv->basis[0][0], diff_order[0], r_val, k);
+        for (int k = 0; k < td->nb_coeffs[1]; k++)
+            basis_val[1][k] = tdi_basis_eval(priv->basis[0][1], diff_order[1], theta_val, k);
+
+        val += scalarproduct_metric_c(td->nb_coeffs[0], td->nb_coeffs[1], priv->coeffs_lapse,
+                                      basis_val[0], basis_val[1]);
+
+        out[i] = val;
+    }
+
+fail:
+    free(basis_val[0]);
+    free(basis_val[1]);
+
+    return 0;
+}
diff --git a/libteukolskydata.v b/libteukolskydata.v
index 50bd496..2ecfdde 100644
--- a/libteukolskydata.v
+++ b/libteukolskydata.v
@@ -1,4 +1,4 @@
-LIBTEUKOLSKYDATA_1 {
+LIBTEUKOLSKYDATA_2 {
     global: td_*;
     local: *;
 };
diff --git a/teukolsky_data.h b/teukolsky_data.h
index d95dbfe..4967f8f 100644
--- a/teukolsky_data.h
+++ b/teukolsky_data.h
@@ -191,5 +191,9 @@ int td_eval_krt(const TDContext *td,
                 size_t nb_coords, const double *r, const double *theta,
                 const unsigned int diff_order[2],
                 double *out);
+int td_eval_lapse(const TDContext *td,
+                  size_t nb_coords, const double *r, const double *theta,
+                  const unsigned int diff_order[2],
+                  double *out);
 
 #endif /* TEUKOLSKY_DATA_H */
diff --git a/teukolsky_data.py b/teukolsky_data.py
index 2458bf7..86db668 100644
--- a/teukolsky_data.py
+++ b/teukolsky_data.py
@@ -124,6 +124,8 @@ class TeukolskyData(object):
         return self._eval_var(self._libtd.td_eval_kpp, r, theta, diff_order)
     def eval_krt(self, r, theta, diff_order = None):
         return self._eval_var(self._libtd.td_eval_krt, r, theta, diff_order)
+    def eval_lapse(self, r, theta, diff_order = None):
+        return self._eval_var(self._libtd.td_eval_lapse, r, theta, diff_order)
 
 
     @property