Rewrite to use the external libbrilldata library.

5 files changed, 62 insertions, 442 deletions

diff --git a/configuration.ccl b/configuration.ccl
index af74fec..7fa37b8 100644
--- a/configuration.ccl
+++ b/configuration.ccl
@@ -1,3 +1 @@
 # Configuration definition for thorn BrillData
-
-REQUIRES GSL
diff --git a/interface.ccl b/interface.ccl
index bdd570c..6c1213a 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,5 +2,3 @@
 implements: BrillData
 
 INHERITS: ADMBase grid CoordBase
-
-CCTK_REAL brill_coeffs TYPE=array DIM=2 SIZE=basis_order_z,basis_order_r
diff --git a/param.ccl b/param.ccl
index 1f80535..2834034 100644
--- a/param.ccl
+++ b/param.ccl
@@ -18,36 +18,17 @@ CCTK_INT eppley_n ""
     1: :: ""
 } 5
 
-CCTK_INT basis_order_r "Number of the basis functions in the radial direction"
+CCTK_INT basis_order_0 "Number of the basis functions in the radial direction"
 {
     1: :: ""
-} 40
+} 80
 
-CCTK_INT basis_order_z "Number of the basis functions in the z direction"
+CCTK_INT basis_order_1 "Number of the basis functions in the angular direction"
 {
     1: :: ""
-} 40
+} 20
 
 CCTK_REAL scale_factor "Scaling factor L for the SB basis"
 {
     0: :: ""
 } 3.0
-
-CCTK_INT overdet "Solve the system overdetermined by this many equations"
-{
-    0: :: ""
-} 0
-
-BOOLEAN export_coeffs "Export the coefficients of the spectral expansion in brill_coeffs"
-{
-} "no"
-
-CCTK_INT colloc_grid_offset_r "Difference between the order of the basis and the collocation grid"
-{
-    0: :: ""
-} 3
-
-CCTK_INT colloc_grid_offset_z "Difference between the order of the basis and the collocation grid"
-{
-    0: :: ""
-} 3
diff --git a/schedule.ccl b/schedule.ccl
index 9f41c25..9ad28f4 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -5,8 +5,4 @@ if (CCTK_Equals(initial_data, "brilldata")) {
     SCHEDULE brill_data IN ADMBase_InitialData {
         LANG: C
     } "Brill wave initial data"
-
-    if (export_coeffs) {
-        STORAGE: brill_coeffs
-    }
 }
diff --git a/src/brill.c b/src/brill.c
index a2cfbd5..2bf40a1 100644
--- a/src/brill.c
+++ b/src/brill.c
@@ -6,8 +6,8 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
-#include <gsl/gsl_blas.h>
-#include <gsl/gsl_linalg.h>
+
+#include <brill_data.h>
 
 #include <lapacke.h>
 
@@ -15,396 +15,23 @@
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
-#define ACC_TEST 1
-
-#define MAX(x, y) ((x) > (y) ? (x) : (y))
-#define MIN(x, y) ((x) > (y) ? (y) : (x))
 #define SQR(x) ((x) * (x))
-#define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
 
 /*
  * small number to avoid r=0 singularities
  */
 #define EPS 1E-08
 
-/* a set of basis functions */
-typedef struct BasisSet {
-    /* evaluate the idx-th basis function at the specified point*/
-    long double (*eval)      (long double coord, int idx);
-    /* evaluate the first derivative of the idx-th basis function at the specified point*/
-    long double (*eval_diff1)(long double coord, int idx);
-    /* evaluate the second derivative of the idx-th basis function at the specified point*/
-    long double (*eval_diff2)(long double coord, int idx);
-    /**
-     * Get the idx-th collocation point for the specified order.
-     * idx runs from 0 to order - 1 (inclusive)
-     */
-    long double (*colloc_point)(int order, int idx);
-} BasisSet;
-
-/*
- * The basis of even (n = 2 * idx) SB functions (Boyd 2000, Ch 17.9)
- * SB(x, n) = sin((n + 1) arccot(|x| / L))
- * They are symmetric wrt origin and decay as 1/x in infinity.
- */
-
-static CCTK_REAL scale_factor;
-
-#define SCALE_FACTOR scale_factor
-
-static long double sb_even_eval(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-
-    idx *= 2;   // even only
-
-    return sinl((idx + 1) * val);
-}
-
-static long double sb_even_eval_diff1(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-
-    idx *= 2;   // even only
-
-    return - SCALE_FACTOR * (idx + 1) * SGN(coord) * cosl((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord));
-}
-
-static long double sb_even_eval_diff2(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-
-    idx *= 2;   // even only
-
-    return SCALE_FACTOR * (idx + 1) * SGN(coord) * (2 * fabsl(coord) * cosl((idx + 1) * val) - SCALE_FACTOR * (idx + 1) * sinl((idx + 1) * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
-}
-
-static long double sb_even_colloc_point(int order, int idx)
-{
-    long double t;
-
-    idx = order - idx - 1;
-    order *= 2;
-
-    t = (idx + 2) * M_PI / (order + 4);
-    return SCALE_FACTOR / tanl(t);
-}
-
-static const BasisSet sb_even_basis = {
-    .eval         = sb_even_eval,
-    .eval_diff1   = sb_even_eval_diff1,
-    .eval_diff2   = sb_even_eval_diff2,
-    .colloc_point = sb_even_colloc_point,
-};
-
-typedef struct BrillDataContext {
-    /* options */
-    CCTK_REAL amplitude;
-    int       eppley_n;
-
-    double (*q_func)(struct BrillDataContext *bd, double rho, double z);
-    double (*laplace_q_func)(struct BrillDataContext *bd, double rho, double z);
-
-    /* state */
-    const BasisSet *basis;
-
-    int nb_coeffs_x;
-    int nb_coeffs_z;
-    int nb_coeffs;
-
-    int nb_colloc_points_x;
-    int nb_colloc_points_z;
-    int nb_colloc_points;
-
-    int colloc_grid_order_x;
-    int colloc_grid_order_z;
-
-    gsl_vector *psi_coeffs;
-} BrillDataContext;
-
-static int brill_construct_matrix(BrillDataContext *bd, gsl_matrix *mat,
-                                  gsl_vector *rhs)
-{
-    CCTK_REAL *basis_val, *basis_dval, *basis_d2val;
-    int idx_coeff_x, idx_coeff_z, idx_grid_x, idx_grid_z;
-
-    /* precompute the basis values we will need */
-    // FIXME: assumes same number of points in x and z directions
-    basis_val   = malloc(sizeof(*basis_val)   * bd->nb_colloc_points_x * bd->nb_coeffs_x);
-    basis_dval  = malloc(sizeof(*basis_dval)  * bd->nb_colloc_points_x * bd->nb_coeffs_x);
-    basis_d2val = malloc(sizeof(*basis_d2val) * bd->nb_colloc_points_x * bd->nb_coeffs_x);
-    for (int i = 0; i < bd->nb_colloc_points_x; i++) {
-        CCTK_REAL coord = bd->basis->colloc_point(bd->colloc_grid_order_x, i);
-        for (int j = 0; j < bd->nb_coeffs_x; j++) {
-            basis_val  [i * bd->nb_coeffs_x + j] = bd->basis->eval(coord, j);
-            basis_dval [i * bd->nb_coeffs_x + j] = bd->basis->eval_diff1(coord, j);
-            basis_d2val[i * bd->nb_coeffs_x + j] = bd->basis->eval_diff2(coord, j);
-        }
-    }
-
-    for (idx_grid_z = 0; idx_grid_z < bd->nb_colloc_points_z; idx_grid_z++) {
-        CCTK_REAL z_physical = bd->basis->colloc_point(bd->colloc_grid_order_z, idx_grid_z);
-
-        for (idx_grid_x = 0; idx_grid_x < bd->nb_colloc_points_x; idx_grid_x++) {
-            CCTK_REAL x_physical = bd->basis->colloc_point(bd->colloc_grid_order_x, idx_grid_x);
-            CCTK_REAL d2q        = bd->laplace_q_func(bd, x_physical, z_physical);
-            int idx_grid         = idx_grid_z * bd->nb_colloc_points_z + idx_grid_x;
-
-            for (idx_coeff_z = 0; idx_coeff_z < bd->nb_coeffs_z; idx_coeff_z++)
-                for (idx_coeff_x = 0; idx_coeff_x < bd->nb_coeffs_x; idx_coeff_x++) {
-                    int idx_coeff = idx_coeff_z * bd->nb_coeffs_z + idx_coeff_x;
-
-                    mat->data[idx_grid + mat->tda * idx_coeff] = basis_d2val[idx_grid_x * bd->nb_coeffs_x + idx_coeff_x] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] * x_physical +
-                                                                 basis_dval [idx_grid_x * bd->nb_coeffs_x + idx_coeff_x] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] +
-                                                                 basis_d2val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] * basis_val[idx_grid_x * bd->nb_coeffs_x + idx_coeff_x] * x_physical +
-                                                                 basis_val  [idx_grid_x * bd->nb_coeffs_x + idx_coeff_x] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] * 0.25 * d2q * x_physical;
-                }
-            rhs->data[idx_grid] = -0.25 * d2q * x_physical;
-        }
-    }
-
-    free(basis_val);
-    free(basis_dval);
-    free(basis_d2val);
-
-    return 0;
-}
-
-static int solve_linear(gsl_matrix *mat, gsl_vector **px, gsl_vector **prhs)
-{
-    int *ipiv;
-    gsl_matrix *mat_f;
-    double cond, ferr, berr, rpivot;
-
-    gsl_vector   *x = *px;
-    gsl_vector *rhs = *prhs;
-    char      equed = 'N';
-    int         ret = 0;
-
-    ipiv  = malloc(mat->size1 * sizeof(*ipiv));
-    mat_f = gsl_matrix_alloc(mat->size1, mat->size2);
-    if (!ipiv || !mat_f) {
-        ret = -ENOMEM;
-        goto fail;
-    }
-
-    LAPACKE_dgesvx(LAPACK_COL_MAJOR, 'N', 'N', mat->size1, 1,
-                   mat->data, mat->tda, mat_f->data, mat_f->tda, ipiv, &equed,
-                   NULL, NULL, rhs->data, rhs->size, x->data, x->size,
-                   &cond, &ferr, &berr, &rpivot);
-
-    CCTK_VInfo(CCTK_THORNSTRING, "LU factorization solution to a %zdx%zd matrix: "
-               "condition number %16.16g; forward error %16.16g backward error %16.16g",
-               mat->size1, mat->size2, cond, ferr, berr);
-fail:
-    gsl_matrix_free(mat_f);
-    free(ipiv);
-
-    return ret;
-}
-
-static int solve_svd(gsl_matrix *mat, gsl_vector **px, gsl_vector **prhs)
-{
-    double *sv;
-    int rank;
-
-    gsl_vector   *x = *px;
-    gsl_vector *rhs = *prhs;
-
-    sv = malloc(sizeof(*sv) * x->size);
-    if (!sv)
-        return -ENOMEM;
-
-    LAPACKE_dgelsd(LAPACK_COL_MAJOR, rhs->size, x->size, 1, mat->data, mat->tda,
-                   rhs->data, rhs->size, sv, -1, &rank);
-
-    CCTK_VInfo(CCTK_THORNSTRING, "Least squares SVD solution to a %zdx%zd matrix: "
-               "rank %d, condition number %16.16g", mat->size1, mat->size2,
-               rank, sv[x->size - 1] / sv[0]);
-
-    gsl_vector_free(x);
-
-    *px         = rhs;
-    (*px)->size = mat->size1;
-    *prhs       = NULL;
-
-    free(sv);
-
-    return 0;
-}
-
-/*
- * Solve the equation
- * Δψ + ¼ ψ (∂²q/∂r² + ∂²q/∂z²) = 0
- * for the coefficients of spectral approximation of ψ:
- * ψ(r, z) = 1 + ΣaᵢⱼTᵢ(r)Tⱼ(z)
- * where i =  { 0, ... , bd->nb_coeffs_x };
- *       j =  { 0, ... , bd->nb_coeffs_z };
- * q(r, z) and Tᵢ(x) are defined by bd->q_func, bd->laplace_q_func and
- * bd->basis.
- *
- * The cofficients are exported in bd->psi_coeffs
- */
-static int brill_solve(BrillDataContext *bd)
-{
-    gsl_matrix *mat = NULL;
-    gsl_vector *coeffs = NULL, *rhs = NULL;
-
-#if ACC_TEST
-    gsl_vector *rhs_bkp = NULL;
-    gsl_matrix *mat_bkp = NULL;
-#endif
-
-    int ret = 0;
-
-    /* allocate and fill the pseudospectral matrix */
-    mat    = gsl_matrix_alloc (bd->nb_coeffs, bd->nb_colloc_points); // inverted order for lapack
-    coeffs = gsl_vector_alloc (bd->nb_coeffs);
-    rhs    = gsl_vector_calloc(bd->nb_colloc_points);
-    if (!mat || !coeffs || !rhs) {
-        ret = -ENOMEM;
-        goto fail;
-    }
-
-    /* fill the matrix */
-    ret = brill_construct_matrix(bd, mat, rhs);
-    if (ret < 0)
-        goto fail;
-
-#if ACC_TEST
-    /* make backups of the matrix and RHS, since they might be destroyed later */
-    mat_bkp = gsl_matrix_alloc(mat->size2, mat->size1);
-    rhs_bkp = gsl_vector_alloc(rhs->size);
-    if (!mat_bkp || !rhs_bkp) {
-        ret = -ENOMEM;
-        goto fail;
-    }
-
-    gsl_vector_memcpy(rhs_bkp, rhs);
-    gsl_matrix_transpose_memcpy(mat_bkp, mat);
-#endif
-
-    /* solve for the coeffs */
-    if (bd->nb_colloc_points > bd->nb_coeffs)
-        ret = solve_svd(mat, &coeffs, &rhs);
-    else
-        ret = solve_linear(mat, &coeffs, &rhs);
-
-    /* export the result to the caller */
-    bd->psi_coeffs = coeffs;
-
-#if ACC_TEST
-    {
-        double min, max, rmin, rmax;
-
-        gsl_vector_minmax(rhs_bkp, &rmin, &rmax);
-        gsl_blas_dgemv(CblasNoTrans, 1.0, mat_bkp, coeffs, -1.0, rhs_bkp);
-        gsl_vector_minmax(rhs_bkp, &min, &max);
-
-        CCTK_VInfo(CCTK_THORNSTRING, "max(|A·x - rhs|) / max(|rhs|): %16.16g",
-                   MAX(fabs(min), fabs(max)) / MAX(fabs(rmin), fabs(rmax)));
-    }
-#endif
-
-fail:
-
-#if ACC_TEST
-    gsl_matrix_free(mat_bkp);
-    gsl_vector_free(rhs_bkp);
-#endif
-
-    if (ret < 0)
-        gsl_vector_free(coeffs);
-
-    gsl_vector_free(rhs);
-    gsl_matrix_free(mat);
-
-    return ret;
-}
-
-// q function form from PHYSICAL REVIEW D 88, 103009 (2013)
-// with σ and ρ_0 hardcoded to 1 for now
-static double q_gundlach(BrillDataContext *bd, double rho, double z)
-{
-    return bd->amplitude * SQR(rho) * exp(- (SQR(rho) + SQR(z)));
-}
-
-static double laplace_q_gundlach(BrillDataContext *bd, double rho, double z)
-{
-    double r2 = SQR(rho);
-    double z2 = SQR(z);
-    double e  = exp(-r2 - z2);
-
-    return 2 * bd->amplitude * e * (1.0 + 2 * r2 * (r2 + z2 - 3));
-}
-
-static double q_eppley(BrillDataContext *bd, double rho, double z)
-{
-    double r2 = SQR(rho) + SQR(z);
-    return bd->amplitude * SQR(rho) / (1.0 + pow(r2, bd->eppley_n / 2.0));
-}
-
-static double laplace_q_eppley(BrillDataContext *bd, double rho, double z)
-{
-    const int    n   = bd->eppley_n;
-    const double r2  = SQR(rho) + SQR(z);
-    const double r2n = pow(r2, n);
-    const double rn  = sqrt(r2n);
-    const double U   = 1.0 / (1 + rn);
-
-    return bd->amplitude * (2 * U - n * (n + 4) * SQR(rho) * SQR(U) * rn / r2 + 2 * SQR(n) * SQR(rho) * SQR(U) * U * r2n /r2);
-}
-
-static BrillDataContext *init_bd(cGH *cctkGH, CCTK_REAL amplitude, int eppley_n,
-                                 int basis_order_r, int basis_order_z,
-                                 int colloc_offset_r, int colloc_offset_z,
-                                 double sf, int overdet)
-{
-    BrillDataContext *bd;
-
-    if (basis_order_r != basis_order_z)
-        CCTK_WARN(0, "Different r and z basis orders are not supported.");
-
-    bd = malloc(sizeof(*bd));
-
-    bd->q_func         = q_gundlach;
-    bd->laplace_q_func = laplace_q_gundlach;
-    //bd->q_func           = q_eppley;
-    //bd->laplace_q_func   = laplace_q_eppley;
-
-    bd->amplitude      = amplitude;
-    bd->eppley_n       = eppley_n;
-
-    bd->basis      = &sb_even_basis;
-
-    bd->nb_coeffs_x = basis_order_r;
-    bd->nb_coeffs_z = basis_order_z;
-
-    bd->nb_coeffs = bd->nb_coeffs_x * bd->nb_coeffs_z;
-
-    bd->nb_colloc_points_x = basis_order_r + overdet;
-    bd->nb_colloc_points_z = basis_order_z + overdet;
-
-    bd->nb_colloc_points = bd->nb_colloc_points_x * bd->nb_colloc_points_z;
-
-    bd->colloc_grid_order_x = basis_order_r + colloc_offset_r;
-    bd->colloc_grid_order_z = basis_order_z + colloc_offset_z;
-
-    scale_factor = sf;
-
-    return bd;
-}
-
 void brill_data(CCTK_ARGUMENTS)
 {
-    static BrillDataContext *prev_bd;
+    static BDContext *prev_bd;
 
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
 
-    BrillDataContext *bd;
-
-    long double *basis_val_r, *basis_val_z;
+    BDContext *bd;
+    double *r_val, *z_val, *psi_val, *q_val;
+    int ret;
 
     int64_t grid_size = CCTK_GFINDEX3D(cctkGH,
                                        cctk_lsh[0] - 1,
@@ -413,14 +40,22 @@ void brill_data(CCTK_ARGUMENTS)
 
     /* on the first run, solve the equation for ψ */
     if (!prev_bd) {
-        bd = init_bd(cctkGH, amplitude, eppley_n, basis_order_r, basis_order_z,
-                     colloc_grid_offset_r, colloc_grid_offset_z,
-                     scale_factor, overdet);
 
-        brill_solve(bd);
+        bd = bd_context_alloc();
+        if (!bd)
+            CCTK_WARN(0, "Memory allocation failed\n");
+
+        bd->q_func_type  = BD_Q_FUNC_GUNDLACH;
+        bd->amplitude    = amplitude;
+        bd->eppley_n     = eppley_n;
+        bd->nb_coeffs[0] = basis_order_0;
+        bd->nb_coeffs[1] = basis_order_1;
+        bd->basis_scale_factor[0] = scale_factor;
+        bd->basis_scale_factor[1] = scale_factor;
 
-        if (export_coeffs)
-            memcpy(brill_coeffs, bd->psi_coeffs->data, sizeof(*brill_coeffs) * bd->nb_coeffs);
+        ret = bd_solve(bd);
+        if (ret < 0)
+            CCTK_WARN(0, "Error solving the Brill wave initial data equation\n");
 
         prev_bd = bd;
     } else
@@ -435,46 +70,51 @@ void brill_data(CCTK_ARGUMENTS)
     memset(gxz, 0, sizeof(*kxz) * grid_size);
     memset(gyz, 0, sizeof(*kyz) * grid_size);
 
-    /* precompute the basis values on the grid points */
-    basis_val_r = malloc(sizeof(*basis_val_r) * bd->nb_coeffs_x * cctk_lsh[1] * cctk_lsh[0]);
-    basis_val_z = malloc(sizeof(*basis_val_z) * bd->nb_coeffs_z * cctk_lsh[2]);
-
-    for (int i = 0; i < cctk_lsh[2]; i++) {
-        CCTK_REAL zz = z[CCTK_GFINDEX3D(cctkGH, 0, 0, i)];
-        for (int j = 0; j < bd->nb_coeffs_z; j++)
-            basis_val_z[i * bd->nb_coeffs_z + j] = bd->basis->eval(zz, j);
-    }
-
+    /* construct the coordinate vectors to be passed to the library */
+    r_val = malloc(sizeof(*r_val) * cctk_lsh[1] * cctk_lsh[0]);
     for (int j = 0; j < cctk_lsh[1]; j++)
         for (int i = 0; i < cctk_lsh[0]; i++) {
             CCTK_REAL xx = x[CCTK_GFINDEX3D(cctkGH, i, j, 0)];
             CCTK_REAL yy = y[CCTK_GFINDEX3D(cctkGH, i, j, 0)];
             CCTK_REAL r = sqrt(SQR(xx) + SQR(yy));
 
-            for (int k = 0; k < bd->nb_coeffs_x; k++)
-                basis_val_r[(j * cctk_lsh[0] + i) * bd->nb_coeffs_x + k] = bd->basis->eval(r, k);
+            r_val[j * cctk_lsh[0] + i] = r;
         }
+    z_val = malloc(sizeof(*z_val) * cctk_lsh[2]);
+    for (int i = 0; i < cctk_lsh[2]; i++)
+        z_val[i] = z[CCTK_GFINDEX3D(cctkGH, 0, 0, i)];
+
+    psi_val = malloc(sizeof(*psi_val) * grid_size);
+    q_val   = malloc(sizeof(*q_val)   * grid_size);
 
 #pragma omp parallel for
-    for (int k = 0; k < cctk_lsh[2]; k++)
-        for (int j = 0; j < cctk_lsh[1]; j++)
-            for (int i = 0; i < cctk_lsh[0]; i++) {
-                int    index = CCTK_GFINDEX3D(cctkGH, i, j, k);
-                long double xx = x[index], yy = y[index], zz = z[index];
+    for (int j = 0; j < cctkGH->cctk_lsh[1]; j++) {
+        double *r_y   = r_val   + j * cctkGH->cctk_lsh[0];
+        double *psi_y = psi_val + j * cctkGH->cctk_lsh[0] * cctkGH->cctk_lsh[2];
+        double *q_y   = q_val   + j * cctkGH->cctk_lsh[0] * cctkGH->cctk_lsh[2];
 
-                long double r2 = SQR(xx) + SQR(yy);
-                long double r  = sqrtl(r2);
+        int err;
 
-                long double q   = bd->q_func(bd, r, zz);
-                long double e2q = expl(2 * q);
+        err = bd_eval_psi(bd, r_y, cctk_lsh[0], z_val, cctk_lsh[2],
+                          (unsigned int[2]){ 0, 0}, psi_y, cctk_lsh[0]);
+        if (err < 0)
+            CCTK_WARN(0, "Error evaluating the conformal factor");
 
-                long double psi = 1.0, psi4;
+        err = bd_eval_q(bd, r_y, cctk_lsh[0], z_val, cctk_lsh[2],
+                        (unsigned int[2]){ 0, 0}, q_y, cctk_lsh[0]);
+        if (err < 0)
+            CCTK_WARN(0, "Error evaluating the q function");
+
+        for (int k = 0; k < cctk_lsh[2]; k++)
+            for (int i = 0; i < cctk_lsh[0]; i++) {
+                int    index = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                double xx = x[index], yy = y[index];
+                double r2 = SQR(xx) + SQR(yy);
 
-                for (int n = 0; n < bd->nb_coeffs_z; n++)
-                    for (int m = 0; m < bd->nb_coeffs_x; m++)
-                        psi += bd->psi_coeffs->data[n * bd->nb_coeffs_z + m] * basis_val_r[(j * cctk_lsh[0] + i) * bd->nb_coeffs_x + m] * basis_val_z[k * bd->nb_coeffs_z + n];
+                double e2q = exp(2 * q_y[k * cctkGH->cctk_lsh[0] + i]);
+                double psi = psi_y[k * cctk_lsh[0] + i];
 
-                psi4 = SQR(SQR(psi));
+                double psi4 = SQR(SQR(psi));
 
                 if (r2 < EPS)
                     r2 = EPS;
@@ -483,5 +123,12 @@ void brill_data(CCTK_ARGUMENTS)
                 gyy[index] = psi4 * (e2q + (1 - e2q) * SQR(xx) / r2);
                 gzz[index] = psi4 * e2q;
                 gxy[index] = psi4 * (-(1.0 - e2q) * xx * yy / r2);
+
             }
+    }
+
+    free(r_val);
+    free(z_val);
+    free(psi_val);
+    free(q_val);
 }