qfunc: disentangle the API from the global solver API

6 files changed, 159 insertions, 97 deletions

diff --git a/eval.c b/eval.c
index 6c70790..98fd0a7 100644
--- a/eval.c
+++ b/eval.c
@@ -1,5 +1,5 @@
 /*
- * Copyright 2014-2015 Anton Khirnov <anton@khirnov.net>
+ * Copyright 2014-2016 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
@@ -27,6 +27,7 @@
 
 #include "brill_data.h"
 #include "internal.h"
+#include "qfunc.h"
 
 int bd_eval_psi(const BDContext *bd,
                 const double *rho, int nb_coords_rho,
@@ -172,14 +173,14 @@ do {                                                                \
  * derivatives as needed */
 #define DO_EVAL_Q_WRAPPER(DO_EVAL, eval_drho, eval_dz, eval_d2rho, eval_d2z, eval_d2rhoz)                                                    \
 do {                                                                                                                                         \
-    const double base = psi4 * exp(2 * s->q_func->q(bd, rrho, zz));                                                                          \
+    const double base = psi4 * exp(2 * bdi_qfunc_eval(s->qfunc, rrho, zz, 0));                                                               \
     double drho, dz, d2rho, d2z, d2rho_z;                                                                                                    \
                                                                                                                                              \
-    if (eval_drho)   drho    = s->q_func->dq_rho(bd, rrho, zz)    + 2 * GRID(dpsi_rho)    / ppsi;                                            \
-    if (eval_dz)     dz      = s->q_func->dq_z(bd, rrho, zz)      + 2 * GRID(dpsi_z)      / ppsi;                                            \
-    if (eval_d2rho)  d2rho   = s->q_func->d2q_rho(bd, rrho, zz)   + 2 * GRID(d2psi_rho)   / ppsi - 2 * SQR(GRID(dpsi_rho) / ppsi);           \
-    if (eval_d2z)    d2z     = s->q_func->d2q_z(bd, rrho, zz)     + 2 * GRID(d2psi_z)     / ppsi - 2 * SQR(GRID(dpsi_z)   / ppsi);           \
-    if (eval_d2rhoz) d2rho_z = s->q_func->d2q_rho_z(bd, rrho, zz) + 2 * GRID(d2psi_rho_z) / ppsi - 2 * GRID(dpsi_rho) * GRID(dpsi_z) / psi2; \
+    if (eval_drho)   drho    = bdi_qfunc_eval(s->qfunc, rrho, zz, 1) + 2 * GRID(dpsi_rho)    / ppsi;                                            \
+    if (eval_dz)     dz      = bdi_qfunc_eval(s->qfunc, rrho, zz, 3) + 2 * GRID(dpsi_z)      / ppsi;                                            \
+    if (eval_d2rho)  d2rho   = bdi_qfunc_eval(s->qfunc, rrho, zz, 2) + 2 * GRID(d2psi_rho)   / ppsi - 2 * SQR(GRID(dpsi_rho) / ppsi);           \
+    if (eval_d2z)    d2z     = bdi_qfunc_eval(s->qfunc, rrho, zz, 6) + 2 * GRID(d2psi_z)     / ppsi - 2 * SQR(GRID(dpsi_z)   / ppsi);           \
+    if (eval_d2rhoz) d2rho_z = bdi_qfunc_eval(s->qfunc, rrho, zz, 4) + 2 * GRID(d2psi_rho_z) / ppsi - 2 * GRID(dpsi_rho) * GRID(dpsi_z) / psi2; \
                                                                                                                                              \
     DO_EVAL;                                                                                                                                 \
 } while (0)
@@ -286,34 +287,19 @@ int bd_eval_q(const BDContext *bd, const double *rho, int nb_coords_rho,
               double *out, unsigned int out_stride)
 {
     BDPriv *s = bd->priv;
-    double (*eval)(const struct BDContext *bd, double rho, double z);
+    int diff_idx;
 
     if (diff_order[0] > 2 || diff_order[1] > 2) {
         bdi_log(bd, 0, "Derivatives of higher order than 2 are not supported.\n");
         return -ENOSYS;
     }
 
-    switch (diff_order[0]) {
-    case 0:
-        switch (diff_order[1]) {
-        case 0: eval = s->q_func->q;     break;
-        case 1: eval = s->q_func->dq_z;  break;
-        case 2: eval = s->q_func->d2q_z; break;
-        }
-        break;
-    case 1:
-        switch (diff_order[1]) {
-        case 0: eval = s->q_func->dq_rho;    break;
-        case 1: eval = s->q_func->d2q_rho_z; break;
-        }
-        break;
-    case 2:     eval = s->q_func->d2q_rho;   break;
-    }
+    diff_idx = diff_order[1] * 3 + diff_order[0];
 
     for (int i = 0; i < nb_coords_z; i++) {
         double *dst = out + i * out_stride;
         for (int j = 0; j < nb_coords_rho; j++)
-            *dst++ = eval(bd, rho[j], z[i]);
+            *dst++ = bdi_qfunc_eval(s->qfunc, rho[j], z[i], diff_idx);
     }
 
     return 0;
diff --git a/init.c b/init.c
index 292c555..5f9e3e6 100644
--- a/init.c
+++ b/init.c
@@ -23,27 +23,17 @@
 
 #include "brill_data.h"
 #include "internal.h"
+#include "qfunc.h"
 
 static int brill_init_check_options(BDContext *bd)
 {
     BDPriv *s = bd->priv;
+    int ret;
 
-    switch (bd->q_func_type) {
-    case BD_Q_FUNC_GUNDLACH:
-        s->q_func   = &bdi_q_func_gundlach;
-        break;
-    case BD_Q_FUNC_EPPLEY:
-        if (bd->eppley_n < 4) {
-            bdi_log(bd, 0, "Invalid n: %d < 4\n", bd->eppley_n);
-            return -EINVAL;
-        }
-
-        s->q_func   = &bdi_q_func_eppley;
-        break;
-    default:
-        bdi_log(bd, 0, "Unknown q function type: %d\n", bd->q_func_type);
-        return -EINVAL;
-    }
+    ret = bdi_qfunc_init(bd, &s->qfunc, bd->q_func_type,
+                         bd->amplitude, bd->eppley_n);
+    if (!s->qfunc)
+        return ret;
 
     s->basis[0] = &bdi_sb_even_basis;
     s->basis[1] = &bdi_sb_even_basis;
@@ -118,6 +108,8 @@ void bd_context_free(BDContext **pbd)
 
     s = bd->priv;
 
+    bdi_qfunc_free(&s->qfunc);
+
     free(s->coeffs);
 
     free(bd->priv);
diff --git a/internal.h b/internal.h
index 19e6d7f..ae9932a 100644
--- a/internal.h
+++ b/internal.h
@@ -22,6 +22,8 @@
 
 #include "brill_data.h"
 
+#include "qfunc.h"
+
 #define MAX(x, y) ((x) > (y) ? (x) : (y))
 #define SQR(x) ((x) * (x))
 
@@ -50,18 +52,9 @@ typedef struct BasisSet {
     double (*colloc_point)(int order, int idx, double sf);
 } BasisSet;
 
-typedef struct QFunc {
-    double         (*q)(const BDContext *bd, double rho, double z);
-    double    (*dq_rho)(const BDContext *bd, double rho, double z);
-    double      (*dq_z)(const BDContext *bd, double rho, double z);
-    double   (*d2q_rho)(const BDContext *bd, double rho, double z);
-    double     (*d2q_z)(const BDContext *bd, double rho, double z);
-    double (*d2q_rho_z)(const BDContext *bd, double rho, double z);
-} QFunc;
-
 typedef struct BDPriv {
     const BasisSet *basis[2];
-    const QFunc    *q_func;
+    QFuncContext    *qfunc;
 
     int nb_colloc_points[2];
     int nb_coeffs[2];
@@ -75,8 +68,6 @@ typedef struct BDPriv {
 
 extern const BasisSet bdi_sb_even_basis;
 
-extern const QFunc bdi_q_func_gundlach;
-extern const QFunc bdi_q_func_eppley;
 
 int bdi_solve(BDContext *bd);
 
diff --git a/qfunc.c b/qfunc.c
index e7827df..82557ba 100644
--- a/qfunc.c
+++ b/qfunc.c
@@ -20,63 +20,74 @@
  * definitions of the q functions in the exponential in the Brill data
  */
 
+#include <errno.h>
 #include <math.h>
+#include <stdlib.h>
 
 #include "brill_data.h"
 #include "internal.h"
 
-static double q_gundlach(const BDContext *bd, double rho, double z)
+typedef double (*QFunc)(const struct QFuncContext *s, double rho, double z);
+
+struct QFuncContext {
+    const QFunc *qfunc;
+
+    double amplitude;
+    int n;
+};
+
+static double q_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return bd->amplitude * SQR(rho) * exp(- (SQR(rho) + SQR(z)));
+    return s->amplitude * SQR(rho) * exp(-(SQR(rho) + SQR(z)));
 }
 
-static double dq_rho_gundlach(const BDContext *bd, double rho, double z)
+static double dq_rho_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return bd->amplitude * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
+    return s->amplitude * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
 }
 
-static double d2q_rho_gundlach(const BDContext *bd, double rho, double z)
+static double d2q_rho_gundlach(const QFuncContext *s, double rho, double z)
 {
     double rho2 = SQR(rho);
-    return bd->amplitude * 2 * exp(-rho2 - SQR(z)) * ((1 - rho2) * (1 - 2 * rho2) - 2 * rho2);
+    return s->amplitude * 2 * exp(-rho2 - SQR(z)) * ((1 - rho2) * (1 - 2 * rho2) - 2 * rho2);
 }
 
-static double d2q_rho_z_gundlach(const BDContext *bd, double rho, double z)
+static double d2q_rho_z_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return -bd->amplitude * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
+    return -s->amplitude * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
 }
 
-static double dq_z_gundlach(const BDContext *bd, double rho, double z)
+static double dq_z_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return - bd->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z));
+    return -s->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z));
 }
 
-static double d2q_z_gundlach(const BDContext *bd, double rho, double z)
+static double d2q_z_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return bd->amplitude * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1);
+    return s->amplitude * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1);
 }
 
 // q function from from PHYSICAL REVIEW D 88, 103009 (2013)
 // with σ and ρ_0 hardcoded to 1 for now
-const QFunc bdi_q_func_gundlach = {
-    .q         = q_gundlach,
-    .dq_rho    = dq_rho_gundlach,
-    .dq_z      = dq_z_gundlach,
-    .d2q_rho   = d2q_rho_gundlach,
-    .d2q_z     = d2q_z_gundlach,
-    .d2q_rho_z = d2q_rho_z_gundlach,
+static const QFunc q_func_gundlach[9] = {
+    [0] = q_gundlach,
+    [1] = dq_rho_gundlach,
+    [3] = dq_z_gundlach,
+    [2] = d2q_rho_gundlach,
+    [6] = d2q_z_gundlach,
+    [4] = d2q_rho_z_gundlach,
 };
 
-static double q_eppley(const BDContext *bd, double rho, double z)
+static double q_eppley(const QFuncContext *s, double rho, double z)
 {
     double r2 = SQR(rho) + SQR(z);
-    return bd->amplitude * SQR(rho) / (1.0 + pow(r2, bd->eppley_n / 2.0));
+    return s->amplitude * SQR(rho) / (1.0 + pow(r2, s->n / 2.0));
 }
 
-static double dq_rho_eppley(const BDContext *bd, double rho, double z)
+static double dq_rho_eppley(const QFuncContext *s, double rho, double z)
 {
-    double A = bd->amplitude;
-    double n = bd->eppley_n;
+    double A = s->amplitude;
+    double n = s->n;
     double rho2 = SQR(rho);
     double z2   = SQR(z);
 
@@ -87,10 +98,10 @@ static double dq_rho_eppley(const BDContext *bd, double rho, double z)
     return A * rho * (2 * (1 + rnm2 * r2) - n * rho2 * rnm2) / SQR(1 + rnm2 * r2);
 }
 
-static double dq_z_eppley(const BDContext *bd, double rho, double z)
+static double dq_z_eppley(const QFuncContext *s, double rho, double z)
 {
-    double A = bd->amplitude;
-    double n = bd->eppley_n;
+    double A = s->amplitude;
+    double n = s->n;
     double rho2 = SQR(rho);
     double z2   = SQR(z);
 
@@ -101,10 +112,10 @@ static double dq_z_eppley(const BDContext *bd, double rho, double z)
     return - A * n * rho2 * z * rnm2 / SQR(1 + rnm2 * r2);
 }
 
-static double d2q_rho_eppley(const BDContext *bd, double rho, double z)
+static double d2q_rho_eppley(const QFuncContext *s, double rho, double z)
 {
-    double A = bd->amplitude;
-    double n = bd->eppley_n;
+    double A = s->amplitude;
+    double n = s->n;
     double rho2 = SQR(rho);
     double z2   = SQR(z);
 
@@ -117,10 +128,10 @@ static double d2q_rho_eppley(const BDContext *bd, double rho, double z)
     return A * (SQR(n) * SQR(rho2) * rnm4 * (rn - 1) - n * rho2 * rnm4 * (3 * rho2 + 5 * z2) * rnp1 + 2 * SQR(rnp1)) / (rnp1 * SQR(rnp1));
 }
 
-static double d2q_z_eppley(const BDContext *bd, double rho, double z)
+static double d2q_z_eppley(const QFuncContext *s, double rho, double z)
 {
-    double A = bd->amplitude;
-    double n = bd->eppley_n;
+    double A = s->amplitude;
+    double n = s->n;
     double rho2 = SQR(rho);
     double z2   = SQR(z);
 
@@ -133,10 +144,10 @@ static double d2q_z_eppley(const BDContext *bd, double rho, double z)
     return A * n * rho2 * rnm4 * (z2 - rho2 - n * z2 + rn * ((1 + n) * z2 - rho2)) / (rnp1 * SQR(rnp1));
 }
 
-static double d2q_rho_z_eppley(const BDContext *bd, double rho, double z)
+static double d2q_rho_z_eppley(const QFuncContext *s, double rho, double z)
 {
-    double A = bd->amplitude;
-    double n = bd->eppley_n;
+    double A = s->amplitude;
+    double n = s->n;
     double rho2 = SQR(rho);
     double z2   = SQR(z);
 
@@ -151,11 +162,60 @@ static double d2q_rho_z_eppley(const BDContext *bd, double rho, double z)
 
 // q function from from PHYSICAL REVIEW D 16, 1609 (1977)
 // with λ hardcoded to 1 for now
-const QFunc bdi_q_func_eppley = {
-    .q         = q_eppley,
-    .dq_rho    = dq_rho_eppley,
-    .dq_z      = dq_z_eppley,
-    .d2q_rho   = d2q_rho_eppley,
-    .d2q_z     = d2q_z_eppley,
-    .d2q_rho_z = d2q_rho_z_eppley,
+static const QFunc q_func_eppley[9] = {
+    [0] = q_eppley,
+    [1] = dq_rho_eppley,
+    [3] = dq_z_eppley,
+    [2] = d2q_rho_eppley,
+    [6] = d2q_z_eppley,
+    [4] = d2q_rho_z_eppley,
 };
+
+double bdi_qfunc_eval(QFuncContext *s, double rho, double z, int diff_idx)
+{
+    return s->qfunc[diff_idx](s, rho, z);
+}
+
+void bdi_qfunc_free(QFuncContext **ps)
+{
+    QFuncContext *s = *ps;
+    free(s);
+    *ps = NULL;
+}
+
+int bdi_qfunc_init(const BDContext *bd, QFuncContext **out, enum BDQFuncType type,
+                   double amplitude, int n)
+{
+    QFuncContext *s = calloc(1, sizeof(*s));
+    int ret = 0;
+
+    if (!s)
+        return -ENOMEM;
+
+    switch (type) {
+    case BD_Q_FUNC_GUNDLACH:
+        s->qfunc = q_func_gundlach;
+        break;
+    case BD_Q_FUNC_EPPLEY:
+        if (n < 4) {
+            bdi_log(bd, 0, "Invalid n: %d < 4\n", n);
+            ret = -EINVAL;
+            goto fail;
+        }
+        s->qfunc = q_func_eppley;
+        break;
+    default:
+        bdi_log(bd, 0, "Unknown q function type: %d\n", type);
+        ret = -EINVAL;
+        goto fail;
+    }
+
+    s->amplitude = amplitude;
+    s->n         = n;
+
+    *out = s;
+    return 0;
+fail:
+    free(s);
+    return ret;
+}
diff --git a/qfunc.h b/qfunc.h
new file mode 100644
index 0000000..be3f60e
--- /dev/null
+++ b/qfunc.h
@@ -0,0 +1,31 @@
+/*
+ * Copyright 2016 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/>.
+ */
+
+#ifndef BRILL_DATA_QFUNC_H
+#define BRILL_DATA_QFUNC_H
+
+#include "brill_data.h"
+
+typedef struct QFuncContext QFuncContext;
+
+int bdi_qfunc_init(const BDContext *bd, QFuncContext **out, enum BDQFuncType type,
+                   double amplitude, int n);
+void bdi_qfunc_free(QFuncContext **ps);
+
+double bdi_qfunc_eval(QFuncContext *s, double rho, double z, int diff_idx);
+
+#endif /* BRILL_DATA_QFUNC_H */
diff --git a/solve.c b/solve.c
index 64dcb55..5803d87 100644
--- a/solve.c
+++ b/solve.c
@@ -32,6 +32,7 @@
 
 #include "brill_data.h"
 #include "internal.h"
+#include "qfunc.h"
 
 static int brill_construct_matrix(const BDContext *bd, double *mat,
                                   double *rhs)
@@ -76,7 +77,8 @@ static int brill_construct_matrix(const BDContext *bd, double *mat,
 
         for (idx_grid_rho = 0; idx_grid_rho < s->nb_colloc_points[0]; idx_grid_rho++) {
             double x_val = s->basis[0]->colloc_point(s->nb_coeffs[0], idx_grid_rho, sf[0]);
-            double d2q   = s->q_func->d2q_rho(bd, x_val, z_val) + s->q_func->d2q_z(bd, x_val, z_val);
+            double d2q = bdi_qfunc_eval(s->qfunc, x_val, z_val, 2) +
+                         bdi_qfunc_eval(s->qfunc, x_val, z_val, 6);
             int idx_grid = idx_grid_z * s->nb_colloc_points[0] + idx_grid_rho;
 
             for (idx_coeff_z = 0; idx_coeff_z < s->nb_coeffs[1]; idx_coeff_z++)