diff options
Diffstat (limited to 'src/brill.c')
| -rw-r--r-- | src/brill.c | 22 |
1 files changed, 17 insertions, 5 deletions
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); |