aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAnton Khirnov <anton@khirnov.net>2014-12-08 14:57:35 +0100
committerAnton Khirnov <anton@khirnov.net>2014-12-08 14:57:35 +0100
commitc412ac65afca2b84b4a3a5fc94a24c8359a2c2e9 (patch)
tree7d37cf258b31396c0caf8ec99e5a721b6d44c0c5
parentff92be5824c046a2e6b290e2bb763d20a4da3805 (diff)
Do not multiply the equation by rho.
Use a different way to avoid the coordinate singularity at origin.
-rw-r--r--brill_data.c13
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;
}
}