Rewrite constraint evaluation code.

Make it much more efficient and easier to add other seed functions.

6 files changed, 397 insertions, 251 deletions

diff --git a/init.c b/init.c
index e90bbb8..d93d5a3 100644
--- a/init.c
+++ b/init.c
@@ -36,6 +36,7 @@
 #include "cpu.h"
 #include "log.h"
 #include "nlsolve.h"
+#include "td_constraints.h"
 #include "teukolsky_data.h"
 #include "threadpool.h"
 
@@ -241,16 +242,47 @@ static int teukolsky_init_check_options(TDContext *td)
     return 0;
 }
 
+static int constraint_eval_alloc(TDContext *td, double amplitude,
+                                 TDConstraintEvalContext **pce)
+{
+    TDPriv *priv = td->priv;
+    TDConstraintEvalContext *ce;
+    int ret;
+
+    ret = tdi_constraint_eval_alloc(&ce, TD_FAMILY_AE_TIME_ANTISYM);
+    if (ret < 0) {
+        tdi_log(&priv->logger, 0, "Error allocating the constraints evaluator\n");
+        return ret;
+    }
+
+    ce->logger       = priv->logger;
+    ce->amplitude    = amplitude;
+    ce->nb_coords[0] = td->nb_coeffs[0];
+    ce->nb_coords[1] = td->nb_coeffs[1];
+    ce->coords[0]    = priv->solver->colloc_grid[0][0];
+    ce->coords[1]    = priv->solver->colloc_grid[0][1];
+
+    ret = tdi_constraint_eval_init(ce);
+    if (ret < 0) {
+        tdi_constraint_eval_free(&ce);
+        tdi_log(&priv->logger, 0, "Error initializing the constraints evaluator\n");
+        return ret;
+    }
+
+    *pce = ce;
+    return 0;
+}
+
 static int teukolsky_constraint_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)
 {
+    TDConstraintEvalContext *ce = opaque;
     const double *amplitude = opaque;
 
-    return tdi_constraints_eval(eq_idx, *amplitude, colloc_grid_order, colloc_grid,
-                                vars, dst);
+    return tdi_constraint_eval(ce, eq_idx, vars, dst);
 }
 
 #define AMPLITUDE_DIRECT_SOLVE 1e-3
@@ -260,19 +292,26 @@ static int teukolsky_constraint_eval(void *opaque, unsigned int eq_idx,
 int td_solve(TDContext *td, double *coeffs_init[3])
 {
     TDPriv *s = td->priv;
+    TDConstraintEvalContext *ce;
+    double a0;
     int ret;
 
     ret = teukolsky_init_check_options(td);
     if (ret < 0)
         return ret;
 
-    if (td->amplitude > AMPLITUDE_CONVERGE) {
+    ret = constraint_eval_alloc(td, 0.0, &ce);
+    if (ret < 0)
+        return ret;
+    if (fabs(td->amplitude) >= ce->a_diverge) {
         tdi_log(&s->logger, 0,
-                "Amplitude A=%16.16g is above the branch point, no solutions "
+                "Amplitude A=%16.16g is above the point A_{max}=%g, no solutions "
                 "are known to exist there. Set solution_branch=1 to get to the "
                 "second branch where mass increases with decreasing amplitude\n",
-                td->amplitude);
+                td->amplitude, ce->a_converge);
     }
+    a0 = ce->a_converge;
+    tdi_constraint_eval_free(&ce);
 
     if (coeffs_init) {
         for (int i = 0; i < 3; i++) {
@@ -283,8 +322,13 @@ int td_solve(TDContext *td, double *coeffs_init[3])
 
     if (td->solution_branch == 0 || coeffs_init) {
         // direct solve with default (flat space) or user-provided initial guess
+        ret = constraint_eval_alloc(td, td->amplitude, &ce);
+        if (ret < 0)
+            return ret;
+
         ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL,
-                                &td->amplitude, s->coeffs, 0);
+                                ce, s->coeffs, 0);
+        tdi_constraint_eval_free(&ce);
         if (ret < 0) {
             tdi_log(&s->logger, 0, "tdi_nlsolve_solve() failed: %d", ret);
             return ret;
@@ -294,18 +338,27 @@ int td_solve(TDContext *td, double *coeffs_init[3])
         // execute two lower-branch solutions and extrapolate to get to the upper branch
         const int N = NB_EQUATIONS * td->nb_coeffs[0] * td->nb_coeffs[1];
         const double delta = 1e-5;
-        const double a0 = AMPLITUDE_CONVERGE;
-        const double a1 = AMPLITUDE_CONVERGE - delta;
+        const double a1 = a0 - delta;
 
         double cur_amplitude, new_amplitude, inverse_step;
         int dir;
 
+        ret = constraint_eval_alloc(td, a0, &ce);
+        if (ret < 0)
+            return ret;
+
         ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL,
-                                &a0, s->coeffs, 0);
+                                ce, s->coeffs, 0);
+        tdi_constraint_eval_free(&ce);
+        if (ret < 0)
+            return ret;
+
+        ret = constraint_eval_alloc(td, a1, &ce);
         if (ret < 0)
             return ret;
         ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL,
-                               &a1, s->coeffs_tmp, 0);
+                                ce, s->coeffs_tmp, 0);
+        tdi_constraint_eval_free(&ce);
         if (ret < 0)
             return ret;
 
@@ -313,8 +366,13 @@ int td_solve(TDContext *td, double *coeffs_init[3])
         cblas_daxpy(N, -1.0, s->coeffs_tmp, 1, s->coeffs,     1);
 
         // obtain solution for a1 in the upper branch
+        ret = constraint_eval_alloc(td, a1, &ce);
+        if (ret < 0)
+            return ret;
+
         ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL,
-                                &a1, s->coeffs, 0);
+                                ce, s->coeffs, 0);
+        tdi_constraint_eval_free(&ce);
         if (ret < 0) {
             tdi_log(&s->logger, 0, "Failed to get into the upper branch\n");
             return ret;
@@ -333,9 +391,15 @@ int td_solve(TDContext *td, double *coeffs_init[3])
             if (dir * (td->amplitude - new_amplitude) < 0)
                 new_amplitude = td->amplitude;
             tdi_log(&s->logger, 2, "Trying amplitude %g\n", new_amplitude);
+
+            ret = constraint_eval_alloc(td, new_amplitude, &ce);
+            if (ret < 0)
+                return ret;
+
             memcpy(s->coeffs_tmp, s->coeffs, sizeof(*s->coeffs) * N);
             ret = tdi_nlsolve_solve(s->solver, teukolsky_constraint_eval, NULL,
-                                    &new_amplitude, s->coeffs_tmp, 1);
+                                    ce, s->coeffs_tmp, 1);
+            tdi_constraint_eval_free(&ce);
             if (ret == -EDOM) {
                 inverse_step = 0.5 * inverse_step;
                 if (fabs(inverse_step) < 1e-2)
@@ -360,7 +424,7 @@ int td_solve(TDContext *td, double *coeffs_init[3])
 
     }
 finish:
-    tdi_solver_print_stats(s->solver);
+    tdi_nlsolve_print_stats(s->solver);
 
     return 0;
 }
diff --git a/nlsolve.c b/nlsolve.c
index 708cd08..2137767 100644
--- a/nlsolve.c
+++ b/nlsolve.c
@@ -235,7 +235,7 @@ finish:
     return ret;
 }
 
-void tdi_solver_print_stats(NLSolveContext *ctx)
+void tdi_nlsolve_print_stats(NLSolveContext *ctx)
 {
     NLSolvePriv *s = ctx->priv;
 
@@ -497,6 +497,8 @@ int tdi_nlsolve_context_init(NLSolveContext *ctx)
         return ret;
     }
 
+    ctx->colloc_grid = s->ps_ctx->colloc_grid;
+
     /* init the per-equation state */
     for (int i = 0; i < s->nb_equations; i++) {
         ret = eq_ctx_init(ctx, i);
diff --git a/nlsolve.h b/nlsolve.h
index f9b69fe..426cb24 100644
--- a/nlsolve.h
+++ b/nlsolve.h
@@ -72,6 +72,17 @@ typedef struct NLSolveContext {
      */
     unsigned int (*solve_order)[2];
 
+    /**
+     * Locations of the collocation points. The equation coefficients in NLEqCallback and
+     * NLEqJacobianCallback should be evaluated at those grid positions.
+     *
+     * colloc_grid[i][j] is an array of length solve_order[i][j] and contains
+     * the collocation points for the i-th equation in the j-th direction.
+     *
+     * Set by the solver after tdi_nlsolve_context_init().
+     */
+    double *(*colloc_grid)[2];
+
 #if HAVE_OPENCL
     cl_context       ocl_ctx;
     cl_command_queue ocl_queue;
diff --git a/td_constraints.c b/td_constraints.c
index 6795891..fe6bc18 100644
--- a/td_constraints.c
+++ b/td_constraints.c
@@ -15,49 +15,275 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <errno.h>
 #include <math.h>
 
 #include "common.h"
 #include "pssolve.h"
 #include "td_constraints.h"
 
-#include "eval_constraints.c"
-#include "eval_k_rtheta.c"
-
-#define EQNAME   constraint_eq_ham
-#define EQUATION TD_CONSTRAINT_EQ_HAM
-#include "td_constraints_template.c"
-#undef EQNAME
-#undef EQUATION
-
-#define EQNAME   constraint_eq_mom_0
-#define EQUATION TD_CONSTRAINT_EQ_MOM_0
-#include "td_constraints_template.c"
-#undef EQNAME
-#undef EQUATION
-
-#define EQNAME   constraint_eq_mom_1
-#define EQUATION TD_CONSTRAINT_EQ_MOM_1
-#include "td_constraints_template.c"
-#undef EQNAME
-#undef EQUATION
-
-void
-(*constraint_funcs[TD_CONSTRAINT_EQ_NB])(const double a,
-                                         unsigned int nb_coords[2], double *coords[2],
+struct TDConstraintEvalPriv {
+    double *k_rtheta;
+    double *dk_rtheta[2];
+};
+
+typedef struct TDFamilyDef {
+    double (*eval_krt)   (TDConstraintEvalContext *ctx, double r, double theta);
+    double (*eval_krt_dr)(TDConstraintEvalContext *ctx, double r, double theta);
+    double (*eval_krt_dt)(TDConstraintEvalContext *ctx, double r, double theta);
+    double a_converge;
+    double a_diverge;
+} TDFamilyDef;
+
+static double k_rtheta_eval_ae_time_antisym(TDConstraintEvalContext *ctx,
+                                            double r, double theta)
+{
+    const double r2 = SQR(r);
+    const double r3 = r * r2;
+    return -60.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI);
+}
+
+static double k_rtheta_eval_dr_ae_time_antisym(TDConstraintEvalContext *ctx,
+                                               double r, double theta)
+{
+    const double r2 = SQR(r);
+    const double r4 = SQR(r2);
+    return -60.0 * sqrt(2) * (3.0 - 12.0 * r2 + 4.0 * r4) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI);
+}
+
+static double k_rtheta_eval_dt_ae_time_antisym(TDConstraintEvalContext *ctx,
+                                               double r, double theta)
+{
+    const double r2 = SQR(r);
+    const double r3 = r * r2;
+    return -120.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * cos(2.0 * theta) * ctx->amplitude / sqrt(M_PI);
+}
+
+static const TDFamilyDef ae_time_antisym = {
+    .eval_krt    = k_rtheta_eval_ae_time_antisym,
+    .eval_krt_dr = k_rtheta_eval_dr_ae_time_antisym,
+    .eval_krt_dt = k_rtheta_eval_dt_ae_time_antisym,
+    .a_converge  = 0.0092078837420,
+    .a_diverge   = 0.0092078837427,
+};
+
+#if 0
+static const TDFamilyDef simple_time_antisym = {
+    .eval_krt    = k_rtheta_eval_simple_time_antisym,
+    .eval_krt_dr = k_rtheta_eval_dr_simple_time_antisym,
+    .eval_krt_dt = k_rtheta_eval_dt_simple_time_antisym,
+    .a_converge  = 1.018918628476058,
+    .a_diverge   = 1.018918628476968,
+};
+#endif
+
+static const TDFamilyDef *td_families[] = {
+    [TD_FAMILY_AE_TIME_ANTISYM] = &ae_time_antisym,
+};
+
+static void
+constraint_eval_ham_confflat(TDConstraintEvalContext *ctx,
+                             double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+                             double *out)
+{
+    TDConstraintEvalPriv *priv = ctx->priv;
+
+    for (int idx1 = 0; idx1 < ctx->nb_coords[1]; idx1++) {
+        const double theta = ctx->coords[1][idx1];
+        const double tt = tan(theta);
+
+        for (int idx0 = 0; idx0 < ctx->nb_coords[0]; idx0++) {
+            const double r = ctx->coords[0][idx0];
+
+            const int idx = idx1 * ctx->nb_coords[0] + idx0;
+
+            const double K_rt = priv->k_rtheta[idx];
+
+            const double K_rr     = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_00][idx];
+            const double K_pp     = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_00][idx];
+            const double psi      = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_00][idx] + 1.0;
+            const double dpsi_r   = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_10][idx];
+            const double dpsi_t   = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_01][idx];
+            const double d2psi_rr = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_20][idx];
+            const double d2psi_tt = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_02][idx];
+
+            const double Km[3][3] = {{ K_rr,                    K_rt,    0 },
+                                     { K_rt / SQR(r), -(K_rr + K_pp),    0 },
+                                     { 0,                          0, K_pp }};
+
+            const double r2   = SQR(r);
+            const double psi2 = SQR(psi);
+            const double psi4 = SQR(psi2);
+            const double Rscal = -8.0 * (d2psi_rr + 2.0 * dpsi_r / r + d2psi_tt / r2 + dpsi_t / (r2 * tt)) / (psi4 * psi);
+
+            double K2 = 0.0;
+            for (int i = 0; i < 3; i++)
+                for (int j = 0; j < 3; j++)
+                    K2 += Km[i][j] * Km[j][i];
+
+            out[idx] = Rscal - K2;
+        }
+    }
+}
+
+static void
+constraint_eval_mom_r_confflat(TDConstraintEvalContext *ctx,
+                               double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+                               double *out)
+{
+    TDConstraintEvalPriv *priv = ctx->priv;
+
+    for (int idx1 = 0; idx1 < ctx->nb_coords[1]; idx1++) {
+        const double theta = ctx->coords[1][idx1];
+        const double tt = tan(theta);
+
+        for (int idx0 = 0; idx0 < ctx->nb_coords[0]; idx0++) {
+            const double r = ctx->coords[0][idx0];
+
+            const int idx = idx1 * ctx->nb_coords[0] + idx0;
+
+            const double K_rt     = priv->k_rtheta[idx];
+            const double dK_rt_t  = priv->dk_rtheta[1][idx];
+
+            const double K_rr     = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_00][idx];
+            const double dK_rr_r  = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_10][idx];
+            const double psi      = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_00][idx] + 1.0;
+            const double dpsi_r   = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_10][idx];
+            const double dpsi_t   = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_01][idx];
+
+            const double r2   = SQR(r);
+
+            out[idx] = 6.0 * (K_rr * dpsi_r + K_rt * dpsi_t / r2) / psi +
+                       dK_rr_r + 3.0 * K_rr / r + K_rt / (r2 * tt) + dK_rt_t / r2;
+        }
+    }
+}
+
+static void
+constraint_eval_mom_t_confflat(TDConstraintEvalContext *ctx,
+                               double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+                               double *out)
+{
+    TDConstraintEvalPriv *priv = ctx->priv;
+
+    for (int idx1 = 0; idx1 < ctx->nb_coords[1]; idx1++) {
+        const double theta = ctx->coords[1][idx1];
+        const double tt = tan(theta);
+
+        for (int idx0 = 0; idx0 < ctx->nb_coords[0]; idx0++) {
+            const double r = ctx->coords[0][idx0];
+
+            const int idx = idx1 * ctx->nb_coords[0] + idx0;
+
+            const double K_rt     = priv->k_rtheta[idx];
+            const double dK_rt_r  = priv->dk_rtheta[0][idx];
+
+            const double K_rr     = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_00][idx];
+            const double dK_rr_t  = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_01][idx];
+            const double K_pp     = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_00][idx];
+            const double dK_pp_t  = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_01][idx];
+            const double psi      = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_00][idx] + 1.0;
+            const double dpsi_r   = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_10][idx];
+            const double dpsi_t   = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_01][idx];
+
+            out[idx] = -(K_rr + 2.0 * K_pp) / tt - 6.0 * (K_rr + K_pp) * dpsi_t / psi +
+                       6.0 * K_rt * dpsi_r / psi - dK_pp_t - dK_rr_t + dK_rt_r + 2.0 * K_rt / r;
+        }
+    }
+}
+
+int tdi_constraint_eval_init(TDConstraintEvalContext *ctx)
+{
+    TDConstraintEvalPriv *priv = ctx->priv;
+    const TDFamilyDef *fd = td_families[ctx->family];
+    int ret;
+
+    ret = posix_memalign((void**)&priv->k_rtheta, 32,
+                         sizeof(*priv->k_rtheta) * ctx->nb_coords[0] * ctx->nb_coords[1]);
+    if (ret)
+        return -ENOMEM;
+
+    for (int i = 0; i < ARRAY_ELEMS(priv->dk_rtheta); i++) {
+        ret = posix_memalign((void**)&priv->dk_rtheta[i], 32,
+                             sizeof(*priv->dk_rtheta[i]) * ctx->nb_coords[0] * ctx->nb_coords[1]);
+        if (ret)
+            return -ENOMEM;
+    }
+
+    ctx->a_converge = fd->a_converge;
+    ctx->a_diverge  = fd->a_diverge;
+
+    for (int idx1 = 0; idx1 < ctx->nb_coords[1]; idx1++) {
+        const double theta = ctx->coords[1][idx1];
+
+        for (int idx0 = 0; idx0 < ctx->nb_coords[0]; idx0++) {
+            const double r = ctx->coords[0][idx0];
+
+            const int idx = idx1 * ctx->nb_coords[0] + idx0;
+
+            priv->k_rtheta[idx]     = fd->eval_krt   (ctx, r, theta);
+            priv->dk_rtheta[0][idx] = fd->eval_krt_dr(ctx, r, theta);
+            priv->dk_rtheta[1][idx] = fd->eval_krt_dt(ctx, r, theta);
+        }
+    }
+
+    return 0;
+}
+
+int tdi_constraint_eval_alloc(TDConstraintEvalContext **pctx, enum TDFamily family)
+{
+    TDConstraintEvalContext *ctx;
+    int ret;
+
+    if (family < 0 || family >= ARRAY_ELEMS(td_families))
+        return -EINVAL;
+
+    ctx = calloc(1, sizeof(*ctx));
+    if (!ctx)
+        return -ENOMEM;
+
+    ctx->priv = calloc(1, sizeof(*ctx->priv));
+    if (!ctx->priv) {
+        tdi_constraint_eval_free(&ctx);
+        return -ENOMEM;
+    }
+
+    ctx->family = family;
+
+    *pctx = ctx;
+    return 0;
+}
+
+void tdi_constraint_eval_free(TDConstraintEvalContext **pctx)
+{
+    TDConstraintEvalContext *ctx = *pctx;
+
+    if (!ctx)
+        return;
+
+    if (ctx->priv) {
+        free(ctx->priv->k_rtheta);
+        free(ctx->priv->dk_rtheta[0]);
+        free(ctx->priv->dk_rtheta[1]);
+    }
+    free(ctx->priv);
+
+    free(ctx);
+    *pctx = NULL;
+}
+
+static const void
+(*constraint_funcs[TD_CONSTRAINT_EQ_NB])(TDConstraintEvalContext *ctx,
                                          double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
                                          double *out) = {
-    [TD_CONSTRAINT_EQ_HAM]   = constraint_eq_ham,
-    [TD_CONSTRAINT_EQ_MOM_0] = constraint_eq_mom_0,
-    [TD_CONSTRAINT_EQ_MOM_1] = constraint_eq_mom_1,
+    [TD_CONSTRAINT_EQ_HAM]   = constraint_eval_ham_confflat,
+    [TD_CONSTRAINT_EQ_MOM_0] = constraint_eval_mom_r_confflat,
+    [TD_CONSTRAINT_EQ_MOM_1] = constraint_eval_mom_t_confflat,
 };
 
-int tdi_constraints_eval(enum TDConstraintEq eq,
-                         const double a,
-                         unsigned int nb_coords[2], double *coords[2],
-                         double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
-                         double *out)
+int tdi_constraint_eval(TDConstraintEvalContext *ctx, enum TDConstraintEq eq,
+                        const double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+                        double *out)
 {
-    constraint_funcs[eq](a, nb_coords, coords, vars, out);
-    return 0;
+    constraint_funcs[eq](ctx, vars, out);
 }
diff --git a/td_constraints.h b/td_constraints.h
index 44f5257..f0a892c 100644
--- a/td_constraints.h
+++ b/td_constraints.h
@@ -20,6 +20,20 @@
 
 #include "pssolve.h"
 
+/**
+ * Identifies the specific initial data family to use.
+ */
+enum TDFamily {
+    /**
+     * The time-antisymmetric initial data used in Abrahams&Evans PhysRevD v49,n8 (1994).
+     * Conformally flat spatial metric.
+     *  r                / x        x  3 \        /    x  2  \
+     * K  = -60√(2/π) a | ---  -  (---)   | * exp| - (---)    | sin(2θ)
+     *  θ                \ L        L    /        \    L     /
+     */
+    TD_FAMILY_AE_TIME_ANTISYM = 0,
+};
+
 enum TDConstraintEq {
     TD_CONSTRAINT_EQ_HAM = 0,
     TD_CONSTRAINT_EQ_MOM_0,
@@ -34,10 +48,37 @@ enum TDConstraintVar {
     TD_CONSTRAINT_VAR_NB,
 };
 
-int tdi_constraints_eval(enum TDConstraintEq eq,
-                         const double a,
-                         unsigned int nb_coords[2], double *coords[2],
-                         double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
-                         double *out);
+
+typedef struct TDConstraintEvalPriv TDConstraintEvalPriv;
+
+typedef struct TDConstraintEvalContext {
+    /**
+     * Private data, not to be touched by the caller.
+     */
+    TDConstraintEvalPriv *priv;
+
+    /**
+     * The logging context.
+     * Set by the caller before tdi_constraint_eval_init().
+     */
+    TDLogger logger;
+
+    unsigned int nb_coords[2];
+    double *coords[2];
+
+    enum TDFamily family;
+
+    double amplitude;
+    double a_converge;
+    double a_diverge;
+} TDConstraintEvalContext;
+
+int tdi_constraint_eval_alloc(TDConstraintEvalContext **ctx, enum TDFamily family);
+int tdi_constraint_eval_init (TDConstraintEvalContext  *ctx);
+void tdi_constraint_eval_free(TDConstraintEvalContext **ctx);
+
+int tdi_constraint_eval(TDConstraintEvalContext *ctx, enum TDConstraintEq eq,
+                        const double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+                        double *out);
 
 #endif // TEUKOLSKY_DATA_CONSTRAINTS_H
diff --git a/td_constraints_template.c b/td_constraints_template.c
deleted file mode 100644
index df09290..0000000
--- a/td_constraints_template.c
+++ /dev/null
@@ -1,198 +0,0 @@
-/*
- * Teukolsky data -- template for the constraint equations definitions
- * Copyright (C) 2017 Anton Khirnov <anton@khirnov.net>
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#define L 1.0
-#define R_0 0.0
-
-static void EQNAME(double a,
-                   unsigned int nb_coords[2], double *coords[2],
-                   double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
-                   double *out)
-{
-    for (int idx1 = 0; idx1 < nb_coords[1]; idx1++) {
-        const double coord1 = coords[1][idx1];
-
-        for (int idx0 = 0; idx0 < nb_coords[0]; idx0++) {
-            const double coord0 = coords[0][idx0];
-
-            const double r     = coord0;
-            const double theta = coord1;
-
-            const int idx = idx1 * nb_coords[0] + idx0;
-
-            const double psi      = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_00][idx] + 1.0;
-            const double K_rr     = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_00][idx];
-            const double K_phiphi = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_00][idx];
-
-            const double dpsi_10 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_10][idx];
-            const double dpsi_01 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_01][idx];
-            const double dpsi_11 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_11][idx];
-            const double dpsi_20 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_20][idx];
-            const double dpsi_02 = vars[TD_CONSTRAINT_VAR_CONFFAC][PSSOLVE_DIFF_ORDER_02][idx];
-
-            const double dK_rr_10      = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_10][idx];
-            const double dK_rr_01      = vars[TD_CONSTRAINT_VAR_CURV_RR][PSSOLVE_DIFF_ORDER_01][idx];
-            const double dK_phiphi_10  = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_10][idx];
-            const double dK_phiphi_01  = vars[TD_CONSTRAINT_VAR_CURV_PHIPHI][PSSOLVE_DIFF_ORDER_01][idx];
-
-            //const double K_rtheta_val = K_rtheta(a, r, theta);
-            //const double dK_rtheta_r_val = K_rtheta_dr(a, r, theta);
-            //const double dK_rtheta_theta_val = K_rtheta_dtheta(a, r, theta);
-            const double K_rtheta_val        = eval_K_rtheta(a, L, R_0, coord0, coord1);
-            const double dK_rtheta_r_val     = eval_dK_rtheta_r(a, L, R_0, coord0, coord1);
-            const double dK_rtheta_theta_val = eval_dK_rtheta_theta(a, L, R_0, coord0, coord1);
-
-            const double Km[3][3] = {{ K_rr,                               K_rtheta_val,       0 },
-                                     { K_rtheta_val / SQR(r), -(K_rr + K_phiphi), 0 },
-                                     { 0,                                  0,                  K_phiphi }};
-
-#if 0
-            double gu[3][3];
-            double dgu[3][3][3], G[3][3][3], dG[3][3][3][3];
-            double Ric[3][3];
-
-            double g[3][3], dg[3][3][3], d2g[3][3][3][3];
-
-            eval_metric                (a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, g);
-            for (int i = 0; i < 3; i++)
-                eval_dmetric[i]        (a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, dg[i]);
-            for (int i = 0; i < 3; i++)
-                for (int j = 0; j < 3; j++)
-                    eval_d2metric[i][j](a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, d2g[i][j]);
-
-            const double det = g[0][0] * g[1][1] * g[2][2] + 2 * g[0][1] * g[1][2] * g[0][2] - g[2][2] * SQR(g[0][1]) - SQR(g[0][2]) * g[1][1] - g[0][0] * SQR(g[1][2]);
-
-            // γ^{ij}
-            gu[0][0] =  (g[1][1] * g[2][2] - SQR(g[1][2])) / det;
-            gu[1][1] =  (g[0][0] * g[2][2] - SQR(g[0][2])) / det;
-            gu[2][2] =  (g[0][0] * g[1][1] - SQR(g[0][1])) / det;
-            gu[0][1] = -(g[0][1] * g[2][2] - g[1][2] * g[0][2]) / det;
-            gu[0][2] =  (g[0][1] * g[1][2] - g[1][1] * g[0][2]) / det;
-            gu[1][2] = -(g[0][0] * g[1][2] - g[0][1] * g[0][2]) / det;
-            gu[1][0] = gu[0][1];
-            gu[2][0] = gu[0][2];
-            gu[2][1] = gu[1][2];
-
-            // ∂_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;
-                    }
-
-            // 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;
-                }
-
-            double Rscal = 0.0;
-            for (int i = 0; i < 3; i++)
-                for (int j = 0; j < 3; j++)
-                    Rscal += gu[i][j] * Ric[i][j];
-
-            const double dKm[3][3][3] = {
-                {
-                    { dK_rr_10, dK_rtheta_r_val, 0 },
-                    { gu[1][1] / gu[0][0] * dK_rtheta_r_val + dgu[0][1][1] / gu[0][0] * K_rtheta_val - dgu[0][0][0] / SQR(gu[0][0]) * gu[1][1] * K_rtheta_val,
-                        -(dK_rr_10 + dK_phiphi_10), 0 },
-                    { 0, 0, dK_phiphi_10 },
-                },
-                {
-                    { dK_rr_01, dK_rtheta_theta_val, 0 },
-                    { gu[1][1] / gu[0][0] * dK_rtheta_theta_val + dgu[1][1][1] / gu[0][0] * K_rtheta_val - dgu[1][0][0] / SQR(gu[0][0]) * gu[1][1] * K_rtheta_val,
-                        -(dK_rr_01 + dK_phiphi_01), 0 },
-                    { 0, 0, dK_phiphi_01 },
-                },
-                { 0 },
-            };
-#endif
-
-            if (EQUATION == TD_CONSTRAINT_EQ_HAM) {
-                double Rscal = eval_R_scalar(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20,
-                                             K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10);
-                double K2 = 0.0;
-                for (int i = 0; i < 3; i++)
-                    for (int j = 0; j < 3; j++)
-                        K2 += Km[i][j] * Km[j][i];
-
-                out[idx] = Rscal - K2;
-            } else if (EQUATION == TD_CONSTRAINT_EQ_MOM_0) {
-                out[idx] = eval_Dcurv_m_r(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20,
-                                          K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10);
-            } else {
-                out[idx] = eval_Dcurv_m_t(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20,
-                                          K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10);
-            }
-#if 0
-            } else {
-                int mom_comp = (EQUATION == TD_CONSTRAINT_EQ_MOM_0) ? 0 : 1;
-
-                double dKm_sum = 0.0;
-                for (int i = 0; i < 3; i++)
-                    dKm_sum += dKm[i][i][mom_comp];
-
-                double Gterm1 = 0.0;
-                for (int i = 0; i < 3; i++)
-                    for (int j = 0; j < 3; j++)
-                        Gterm1 += G[i][i][j] * Km[j][mom_comp];
-
-                double Gterm2 = 0.0;
-                for (int i = 0; i < 3; i++)
-                    for (int j = 0; j < 3; j++)
-                        Gterm2 += G[j][i][mom_comp] * Km[i][j];
-
-                out[idx] = dKm_sum + Gterm1 - Gterm2;
-            }
-#endif
-        }
-    }
-}