diff options
Diffstat (limited to 'solve.c')
| -rw-r--r-- | solve.c | 102 |
1 files changed, 49 insertions, 53 deletions
@@ -37,80 +37,76 @@ static int brill_construct_matrix(const BDContext *bd, double *mat, double *rhs) { BDPriv *s = bd->priv; - const double sfr = bd->basis_scale_factor_rho, sfz = bd->basis_scale_factor_z; + const double *sf = bd->basis_scale_factor; double *basis_val, *basis_dval, *basis_d2val; int idx_coeff_rho, idx_coeff_z, idx_grid_rho, idx_grid_z; - double *basis_r[3]; - double *basis_z[3]; + double *basis[2][3] = { { NULL } }; + int ret = 0; /* precompute the basis values we will need */ - for (int i = 0; i < ARRAY_ELEMS(basis_r); i++) - basis_r[i] = malloc(sizeof(*basis_r[i]) * s->nb_colloc_points_rho * bd->nb_coeffs_rho); - for (int i = 0; i < s->nb_colloc_points_rho; i++) { - double coord = s->basis->colloc_point(s->colloc_grid_order_rho, i, sfr); - for (int j = 0; j < bd->nb_coeffs_rho; j++) { - basis_r[0][i * bd->nb_coeffs_rho + j] = s->basis->eval(coord, j, sfr); - basis_r[1][i * bd->nb_coeffs_rho + j] = s->basis->eval_diff1(coord, j, sfr); - basis_r[2][i * bd->nb_coeffs_rho + j] = s->basis->eval_diff2(coord, j, sfr); + for (int i = 0; i < ARRAY_ELEMS(basis); i++) { + for (int j = 0; j < ARRAY_ELEMS(basis[i]); j++) { + basis[i][j] = malloc(sizeof(*basis[i][j]) * s->nb_colloc_points[i] * s->nb_coeffs[i]); + if (!basis[i][j]) { + ret = -ENOMEM; + goto fail; + } } - } - - for (int i = 0; i < ARRAY_ELEMS(basis_z); i++) - basis_z[i] = malloc(sizeof(*basis_z[i]) * s->nb_colloc_points_z * bd->nb_coeffs_z); - for (int i = 0; i < s->nb_colloc_points_z; i++) { - double coord = s->basis->colloc_point(s->colloc_grid_order_z, i, sfr); - for (int j = 0; j < bd->nb_coeffs_z; j++) { - basis_z[0][i * bd->nb_coeffs_z + j] = s->basis->eval(coord, j, sfr); - basis_z[1][i * bd->nb_coeffs_z + j] = s->basis->eval_diff1(coord, j, sfr); - basis_z[2][i * bd->nb_coeffs_z + j] = s->basis->eval_diff2(coord, j, sfr); + for (int j = 0; j < s->nb_colloc_points[i]; j++) { + double coord = s->basis[i]->colloc_point(s->nb_coeffs[i], j, sf[i]); + for (int k = 0; k < s->nb_coeffs[i]; k++) { + basis[i][0][j * s->nb_coeffs[i] + k] = s->basis[i]->eval (coord, k, sf[i]); + basis[i][1][j * s->nb_coeffs[i] + k] = s->basis[i]->eval_diff1(coord, k, sf[i]); + basis[i][2][j * s->nb_coeffs[i] + k] = s->basis[i]->eval_diff2(coord, k, sf[i]); + } } } -#define BASIS_RHO (basis_r[0][idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho]) -#define DBASIS_RHO (basis_r[1][idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho]) -#define D2BASIS_RHO (basis_r[2][idx_grid_rho * bd->nb_coeffs_rho + idx_coeff_rho]) -#define BASIS_Z (basis_z[0][idx_grid_z * bd->nb_coeffs_z + idx_coeff_z]) -#define DBASIS_Z (basis_z[1][idx_grid_z * bd->nb_coeffs_z + idx_coeff_z]) -#define D2BASIS_Z (basis_z[2][idx_grid_z * bd->nb_coeffs_z + idx_coeff_z]) +#define BASIS_RHO (basis[0][0][idx_grid_rho * s->nb_coeffs[0] + idx_coeff_rho]) +#define DBASIS_RHO (basis[0][1][idx_grid_rho * s->nb_coeffs[0] + idx_coeff_rho]) +#define D2BASIS_RHO (basis[0][2][idx_grid_rho * s->nb_coeffs[0] + idx_coeff_rho]) +#define BASIS_Z (basis[1][0][idx_grid_z * s->nb_coeffs[1] + idx_coeff_z]) +#define DBASIS_Z (basis[1][1][idx_grid_z * s->nb_coeffs[1] + idx_coeff_z]) +#define D2BASIS_Z (basis[1][2][idx_grid_z * s->nb_coeffs[1] + idx_coeff_z]) - for (idx_grid_z = 0; idx_grid_z < s->nb_colloc_points_z; idx_grid_z++) { - double z_physical = s->basis->colloc_point(s->colloc_grid_order_z, idx_grid_z, sfz); + for (idx_grid_z = 0; idx_grid_z < s->nb_colloc_points[1]; idx_grid_z++) { + double z_val = s->basis[1]->colloc_point(s->nb_coeffs[1], idx_grid_z, sf[1]); - for (idx_grid_rho = 0; idx_grid_rho < s->nb_colloc_points_rho; idx_grid_rho++) { - double x_physical = s->basis->colloc_point(s->colloc_grid_order_rho, idx_grid_rho, sfr); - double d2q = s->q_func->d2q_rho(bd, x_physical, z_physical) + s->q_func->d2q_z(bd, x_physical, z_physical); - int idx_grid = idx_grid_z * s->nb_colloc_points_rho + idx_grid_rho; + for (idx_grid_rho = 0; idx_grid_rho < s->nb_colloc_points[0]; idx_grid_rho++) { + double x_val = s->basis[0]->colloc_point(s->nb_coeffs[0], idx_grid_rho, sf[0]); + double d2q = s->q_func->d2q_rho(bd, x_val, z_val) + s->q_func->d2q_z(bd, x_val, z_val); + int idx_grid = idx_grid_z * s->nb_colloc_points[0] + idx_grid_rho; - for (idx_coeff_z = 0; idx_coeff_z < bd->nb_coeffs_z; idx_coeff_z++) - for (idx_coeff_rho = 0; idx_coeff_rho < bd->nb_coeffs_rho; idx_coeff_rho++) { - int idx_coeff = idx_coeff_z * bd->nb_coeffs_rho + idx_coeff_rho; + for (idx_coeff_z = 0; idx_coeff_z < s->nb_coeffs[1]; idx_coeff_z++) + for (idx_coeff_rho = 0; idx_coeff_rho < s->nb_coeffs[0]; idx_coeff_rho++) { + int idx_coeff = idx_coeff_z * s->nb_coeffs[0] + idx_coeff_rho; double val = D2BASIS_RHO * BASIS_Z + D2BASIS_Z * BASIS_RHO + BASIS_RHO * BASIS_Z * 0.25 * d2q; - if (fabs(x_physical) > EPS) - val += DBASIS_RHO * BASIS_Z / fabs(x_physical); + if (fabs(x_val) > EPS) + val += DBASIS_RHO * BASIS_Z / fabs(x_val); else val += D2BASIS_RHO * BASIS_Z; - mat[idx_grid + s->nb_colloc_points * idx_coeff] = val; + mat[idx_grid + NB_COLLOC_POINTS(s) * idx_coeff] = val; } rhs[idx_grid] = -0.25 * d2q; } } - for (int i = 0; i < ARRAY_ELEMS(basis_r); i++) - free(basis_r[i]); - for (int i = 0; i < ARRAY_ELEMS(basis_z); i++) - free(basis_z[i]); +fail: + for (int i = 0; i < ARRAY_ELEMS(basis); i++) + for (int j = 0; j < ARRAY_ELEMS(basis[i]); j++) + free(basis[i][j]); - return 0; + return ret; } static int brill_solve_linear(const BDContext *bd, double *mat, double **px, double **prhs) { const BDPriv *s = bd->priv; - const int stride = s->nb_coeffs; + const int stride = NB_COEFFS(s); int *ipiv; double *mat_f; double cond, ferr, berr, rpivot; @@ -152,18 +148,18 @@ static int brill_solve_svd(const BDContext *bd, double *mat, double *x = *px; double *rhs = *prhs; - sv = malloc(sizeof(*sv) * s->nb_coeffs); + sv = malloc(sizeof(*sv) * NB_COEFFS(s)); if (!sv) return -ENOMEM; - LAPACKE_dgelsd(LAPACK_COL_MAJOR, s->nb_colloc_points, s->nb_coeffs, 1, - mat, s->nb_colloc_points, rhs, - MAX(s->nb_coeffs, s->nb_colloc_points), + LAPACKE_dgelsd(LAPACK_COL_MAJOR, NB_COLLOC_POINTS(s), NB_COEFFS(s), 1, + mat, NB_COLLOC_POINTS(s), rhs, + MAX(NB_COEFFS(s), NB_COLLOC_POINTS(s)), sv, -1, &rank); bdi_log(bd, 0, "Least squares SVD solution to a %zdx%zd matrix: " - "rank %d, condition number %16.16g\n", s->nb_colloc_points, s->nb_coeffs, - rank, sv[s->nb_coeffs - 1] / sv[0]); + "rank %d, condition number %16.16g\n", NB_COLLOC_POINTS(s), NB_COEFFS(s), + rank, sv[NB_COEFFS(s) - 1] / sv[0]); free(sv); @@ -189,7 +185,7 @@ int bdi_solve(BDContext *bd) { BDPriv *s = bd->priv; - const int vecsize = MAX(s->nb_coeffs, s->nb_colloc_points); + const int vecsize = MAX(NB_COEFFS(s), NB_COLLOC_POINTS(s)); double *mat = NULL; double *rhs = NULL, *coeffs = NULL; @@ -197,7 +193,7 @@ int bdi_solve(BDContext *bd) int ret = 0; /* allocate and fill the pseudospectral matrix */ - mat = malloc(sizeof(*mat) * s->nb_coeffs * s->nb_colloc_points); + mat = malloc(sizeof(*mat) * NB_COEFFS(s) * NB_COLLOC_POINTS(s)); coeffs = malloc(sizeof(*coeffs) * vecsize); rhs = malloc(sizeof(*rhs) * vecsize); if (!mat || !coeffs || !rhs) { |