From c412ac65afca2b84b4a3a5fc94a24c8359a2c2e9 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Mon, 8 Dec 2014 14:57:35 +0100 Subject: Do not multiply the equation by rho. Use a different way to avoid the coordinate singularity at origin. --- brill_data.c | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/brill_data.c b/brill_data.c index 40b2082..b1e2230 100644 --- a/brill_data.c +++ b/brill_data.c @@ -165,12 +165,15 @@ static int brill_construct_matrix(BDContext *bd, gsl_matrix *mat, for (idx_coeff_rho = 0; idx_coeff_rho < bd->nb_coeffs_rho; idx_coeff_rho++) { int idx_coeff = idx_coeff_z * bd->nb_coeffs_z + idx_coeff_rho; - mat->data[idx_grid + mat->tda * idx_coeff] = basis_d2val[idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] * x_physical + - basis_dval [idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] * 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_rho * bd->nb_coeffs_rho + idx_coeff_rho] * x_physical + - basis_val [idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] * 0.25 * d2q * x_physical; + mat->data[idx_grid + mat->tda * idx_coeff] = basis_d2val[idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] * 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_rho * bd->nb_coeffs_rho + idx_coeff_rho] + + basis_val [idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] * 0.25 * d2q; + if (x_physical > EPS) + mat->data[idx_grid + mat->tda * idx_coeff] += basis_dval [idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z] / x_physical; + else + mat->data[idx_grid + mat->tda * idx_coeff] += basis_d2val [idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho] * basis_val[idx_grid_z * bd->nb_coeffs_z + idx_coeff_z]; } - rhs->data[idx_grid] = -0.25 * d2q * x_physical; + rhs->data[idx_grid] = -0.25 * d2q; } } -- cgit v1.2.3