From 9a54962f3ff3d17daf2d29d311f063d88887b803 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 11 Sep 2014 13:07:47 +0200 Subject: Allow using higher order collocation grids. Default to 3, since it appears to give much better results. --- param.ccl | 10 ++++++++++ src/brill.c | 22 +++++++++++++++++----- 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); -- cgit v1.2.3