diff options
author | Anton Khirnov <anton@khirnov.net> | 2014-12-08 14:57:35 +0100 |
---|---|---|
committer | Anton Khirnov <anton@khirnov.net> | 2014-12-08 14:57:35 +0100 |
commit | c412ac65afca2b84b4a3a5fc94a24c8359a2c2e9 (patch) | |
tree | 7d37cf258b31396c0caf8ec99e5a721b6d44c0c5 /brill_data.c | |
parent | ff92be5824c046a2e6b290e2bb763d20a4da3805 (diff) |
Do not multiply the equation by rho.
Use a different way to avoid the coordinate singularity at origin.
Diffstat (limited to 'brill_data.c')
-rw-r--r-- | brill_data.c | 13 |
1 files 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; } } |