Allow using higher order collocation grids.

Default to 3, since it appears to give much better results.

2 files changed, 27 insertions, 5 deletions

diff --git a/param.ccl b/param.ccl
index 6984c80..9724012 100644
--- a/param.ccl
+++ b/param.ccl
@@ -36,3 +36,13 @@ CCTK_INT overdet "Solve the system overdetermined by this many equations"
 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/src/brill.c b/src/brill.c
index 952287e..b48980d 100644
--- a/src/brill.c
+++ b/src/brill.c
@@ -82,6 +82,8 @@ static long double sb_even_eval_diff2(long double coord, int idx)
 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);
@@ -113,6 +115,9 @@ typedef struct BrillDataContext {
     int nb_colloc_points_z;
     int nb_colloc_points;
 
+    int colloc_grid_order_x;
+    int colloc_grid_order_z;
+
     gsl_vector *psi_coeffs;
 } BrillDataContext;
 
@@ -128,7 +133,7 @@ static int brill_construct_matrix(BrillDataContext *bd, gsl_matrix *mat,
     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->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);
@@ -137,10 +142,10 @@ static int brill_construct_matrix(BrillDataContext *bd, gsl_matrix *mat,
     }
 
     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->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->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;
 
@@ -326,7 +331,9 @@ static double laplace_q_gundlach(BrillDataContext *bd, double rho, double z)
 }
 
 static BrillDataContext *init_bd(cGH *cctkGH, CCTK_REAL amplitude,
-                                 int basis_order_r, int basis_order_z, double sf, int overdet)
+                                 int basis_order_r, int basis_order_z,
+                                 int colloc_offset_r, int colloc_offset_z,
+                                 double sf, int overdet)
 {
     BrillDataContext *bd;
 
@@ -351,6 +358,9 @@ static BrillDataContext *init_bd(cGH *cctkGH, CCTK_REAL amplitude,
 
     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;
@@ -374,7 +384,9 @@ void brill_data(CCTK_ARGUMENTS)
 
     /* on the first run, solve the equation for ψ */
     if (!prev_bd) {
-        bd = init_bd(cctkGH, amplitude, basis_order_r, basis_order_z, scale_factor, overdet);
+        bd = init_bd(cctkGH, amplitude, basis_order_r, basis_order_z,
+                     colloc_grid_offset_r, colloc_grid_offset_z,
+                     scale_factor, overdet);
 
         brill_solve(bd);