From 788f68c3bbf7c74b717d8894be4610d39d0d97f2 Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Mon, 23 Apr 2018 08:45:07 +0200
Subject: Add the simple time-antisym initial data family.

---
 init.c             | 12 +++++++-----
 libteukolskydata.v |  2 +-
 td_constraints.c   | 25 +++++++++++++++++++++++--
 td_constraints.h   | 15 +--------------
 teukolsky_data.h   | 24 ++++++++++++++++++++++++
 5 files changed, 56 insertions(+), 22 deletions(-)

diff --git a/init.c b/init.c
index a7d54a8..ff995b6 100644
--- a/init.c
+++ b/init.c
@@ -249,14 +249,14 @@ static int teukolsky_init_check_options(TDContext *td)
     return 0;
 }
 
-static int constraint_eval_alloc(TDContext *td, double amplitude,
+static int constraint_eval_alloc(const 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);
+    ret = tdi_constraint_eval_alloc(&ce, td->family);
     if (ret < 0) {
         tdi_log(&priv->logger, 0, "Error allocating the constraints evaluator\n");
         return ret;
@@ -405,8 +405,8 @@ int td_solve(TDContext *td, double *coeffs_init[3])
             tdi_constraint_eval_free(&ce);
             if (ret == -EDOM) {
                 inverse_step = 0.5 * inverse_step;
-                if (fabs(inverse_step) < 1e-2)
-                    return ret;
+                //if (fabs(inverse_step) < 1e-2)
+                //    return ret;
                 continue;
             } else if (ret < 0)
                 return ret;
@@ -446,7 +446,9 @@ TDContext *td_context_alloc(void)
 
     td->nb_threads = 1;
 
-    td->amplitude   = 1.0;
+    td->family          = TD_FAMILY_SIMPLE_TIME_ANTISYM;
+    td->solution_branch = 0;
+    td->amplitude       = 1.0;
 
     td->nb_coeffs[0] = 32;
     td->nb_coeffs[1] = 16;
diff --git a/libteukolskydata.v b/libteukolskydata.v
index 18295a2..50bd496 100644
--- a/libteukolskydata.v
+++ b/libteukolskydata.v
@@ -1,4 +1,4 @@
-LIBTEUKOLSKYDATA_0 {
+LIBTEUKOLSKYDATA_1 {
     global: td_*;
     local: *;
 };
diff --git a/td_constraints.c b/td_constraints.c
index 663b83a..dd62150 100644
--- a/td_constraints.c
+++ b/td_constraints.c
@@ -67,7 +67,27 @@ static const TDFamilyDef ae_time_antisym = {
     .a_diverge   = 0.0092078837427,
 };
 
-#if 0
+static double k_rtheta_eval_simple_time_antisym(TDConstraintEvalContext *ctx,
+                                                double r, double theta)
+{
+    const double r2 = SQR(r);
+    return ctx->amplitude * r * exp(-r2) * sin(2.0 * theta);
+}
+
+static double k_rtheta_eval_dr_simple_time_antisym(TDConstraintEvalContext *ctx,
+                                                   double r, double theta)
+{
+    const double r2 = SQR(r);
+    return ctx->amplitude * (1.0 - 2.0 * r2) * exp(-r2) * sin(2.0 * theta);
+}
+
+static double k_rtheta_eval_dt_simple_time_antisym(TDConstraintEvalContext *ctx,
+                                                   double r, double theta)
+{
+    const double r2 = SQR(r);
+    return 2.0 * ctx->amplitude * r * exp(-r2) * cos(2.0 * theta);
+}
+
 static const TDFamilyDef simple_time_antisym = {
     .eval_krt    = k_rtheta_eval_simple_time_antisym,
     .eval_krt_dr = k_rtheta_eval_dr_simple_time_antisym,
@@ -75,10 +95,10 @@ static const TDFamilyDef 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,
+    [TD_FAMILY_SIMPLE_TIME_ANTISYM] = &simple_time_antisym,
 };
 
 static void
@@ -292,4 +312,5 @@ int tdi_constraint_eval(TDConstraintEvalContext *ctx, enum TDConstraintEq eq,
                         double *out)
 {
     constraint_funcs[eq](ctx, vars, out);
+    return 0;
 }
diff --git a/td_constraints.h b/td_constraints.h
index 6746c3e..259cbfb 100644
--- a/td_constraints.h
+++ b/td_constraints.h
@@ -19,20 +19,7 @@
 #define TEUKOLSKY_DATA_CONSTRAINTS_H
 
 #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,
-};
+#include "teukolsky_data.h"
 
 enum TDConstraintEq {
     TD_CONSTRAINT_EQ_HAM = 0,
diff --git a/teukolsky_data.h b/teukolsky_data.h
index e0877cd..d95dbfe 100644
--- a/teukolsky_data.h
+++ b/teukolsky_data.h
@@ -36,6 +36,28 @@
  * Finally, free the solver context with td_context_free().
  */
 
+/**
+ * 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,
+    /**
+     * Simpler time-antisymmetric initial data.
+     * Conformally flat spatial metric.
+     *  r       x       /    x  2  \
+     * K  = a  ---  exp| - (---)    | sin(2θ)
+     *  θ       L       \    L     /
+     */
+    TD_FAMILY_SIMPLE_TIME_ANTISYM,
+};
+
 typedef struct TDContext {
     /**
      *  Solver internals, not to be accessed by the caller
@@ -105,6 +127,8 @@ typedef struct TDContext {
     double *coeffs[3];
 
     unsigned int solution_branch;
+
+    enum TDFamily family;
 } TDContext;
 
 /**
-- 
cgit v1.2.3