diff options
Diffstat (limited to 'eval.c')
| -rw-r--r-- | eval.c | 122 |
1 files changed, 2 insertions, 120 deletions
@@ -29,124 +29,6 @@ #include "internal.h" #include "qfunc.h" -static int eval_psi_radial(const BDContext *bd, - const double *rho, int nb_coords_rho, - const double *z, int nb_coords_z, - const unsigned int diff_order[2], - double *psi, unsigned int psi_stride) -{ - BDPriv *s = bd->priv; - enum BSEvalType type; - - double *basis_val[2][3] = { { NULL } }; - - double add = 0.0; - - int ret = 0; - - if (diff_order[0] == 0 && diff_order[1] == 0) - add = 1.0; - - /* precompute the basis values on the grid points */ - for (int i = 0; i < ARRAY_ELEMS(basis_val[0]); i++) { - posix_memalign((void**)&basis_val[0][i], REQ_ALIGNMENT(char), sizeof(*basis_val[0][i]) * s->nb_coeffs[0]); - posix_memalign((void**)&basis_val[1][i], REQ_ALIGNMENT(char), sizeof(*basis_val[1][i]) * s->nb_coeffs[1]); - if (!basis_val[0][i] || !basis_val[1][i]) { - ret = -ENOMEM; - goto fail; - } - } - -#define BASIS0 (basis_val[0][0][n]) -#define DBASIS0 (basis_val[0][1][n]) -#define D2BASIS0 (basis_val[0][2][n]) -#define BASIS1 (basis_val[1][0][m]) -#define DBASIS1 (basis_val[1][1][m]) -#define D2BASIS1 (basis_val[1][2][m]) - - for (int i = 0; i < nb_coords_z; i++) { - double *dst = psi + i * psi_stride; - double zz = z[i]; - - for (int j = 0; j < nb_coords_rho; j++) { - double ppsi = add; - double rrho = rho[j]; - - double r = sqrt(SQR(rrho) + SQR(zz)); - double phi = atan2(zz, rrho); - - bdi_basis_eval_multi(s->basis[0], r, basis_val[0], 0, s->nb_coeffs[0]); - bdi_basis_eval_multi(s->basis[1], phi, basis_val[1], 0, s->nb_coeffs[1]); - - if (diff_order[0] == 0 && diff_order[1] == 0) { - for (int m = 0; m < s->nb_coeffs[1]; m++) { - double tmp = 0.0; - for (int n = 0; n < s->nb_coeffs[0]; n++) { - tmp += s->coeffs[m * s->coeffs_stride + n] * BASIS0; - } - ppsi += tmp * BASIS1; - } - } else if (diff_order[0] == 2 && diff_order[1] == 0) { - for (int m = 0; m < s->nb_coeffs[1]; m++) { - for (int n = 0; n < s->nb_coeffs[0]; n++) { - double basis_val_20 = ((r > EPS) ? (SQR(rrho / r) * D2BASIS0 * BASIS1 + SQR(zz / SQR(r)) * BASIS0 * D2BASIS1 - + (1.0 - SQR(rrho / r)) / r * DBASIS0 * BASIS1 - + 2 * rrho * zz / SQR(SQR(r)) * BASIS0 * DBASIS1 - - 2 * zz * rrho / (r * SQR(r)) * DBASIS0 * DBASIS1) : 0.0); - ppsi += s->coeffs[m * s->coeffs_stride + n] * basis_val_20; - } - } - } else if (diff_order[0] == 0 && diff_order[1] == 2) { - for (int m = 0; m < s->nb_coeffs[1]; m++) { - for (int n = 0; n < s->nb_coeffs[0]; n++) { - double basis_val_02 = ((r > EPS) ? (SQR(zz / r) * D2BASIS0 * BASIS1 + SQR(rrho / SQR(r)) * BASIS0 * D2BASIS1 - + (1.0 - SQR(zz / r)) / r * DBASIS0 * BASIS1 - - 2 * rrho * zz / SQR(SQR(r)) * BASIS0 * DBASIS1 - + 2 * zz * rrho / (r * SQR(r)) * DBASIS0 * DBASIS1) : 0.0); - ppsi += s->coeffs[m * s->coeffs_stride + n] * basis_val_02; - } - } - } else if (diff_order[0] == 1 && diff_order[1] == 0) { - for (int m = 0; m < s->nb_coeffs[1]; m++) { - for (int n = 0; n < s->nb_coeffs[0]; n++) { - double basis_val_10 = ((r > EPS) ? (DBASIS0 * BASIS1 * rrho / r - BASIS0 * DBASIS1 * zz / SQR(r)) : 0.0); - - ppsi += s->coeffs[m * s->coeffs_stride + n] * basis_val_10; - } - } - } else if (diff_order[0] == 0 && diff_order[1] == 1) { - for (int m = 0; m < s->nb_coeffs[1]; m++) { - for (int n = 0; n < s->nb_coeffs[0]; n++) { - double basis_val_01 = ((r > EPS) ? (DBASIS0 * BASIS1 * zz / r + BASIS0 * DBASIS1 * rrho / SQR(r)) : 0.0); - - ppsi += s->coeffs[m * s->coeffs_stride + n] * basis_val_01; - } - } - } else if (diff_order[0] == 1 && diff_order[1] == 1) { - for (int m = 0; m < s->nb_coeffs[1]; m++) { - for (int n = 0; n < s->nb_coeffs[0]; n++) { - double basis_val_11 = ((r > EPS) ? (rrho * zz / SQR(r) * D2BASIS0 * BASIS1 - rrho * zz / SQR(SQR(r)) * BASIS0 * D2BASIS1 - - rrho * zz / (r * SQR(r)) * DBASIS0 * BASIS1 - - (1.0 - SQR(zz / r)) / SQR(r) * BASIS0 * DBASIS1 - + (SQR(rrho) - SQR(zz)) / (r * SQR(r)) * DBASIS0 * DBASIS1) : 0.0); - - ppsi += s->coeffs[m * s->coeffs_stride + n] * basis_val_11; - } - } - } - *dst++ = ppsi; - } - } - -fail: - for (int i = 0; i < ARRAY_ELEMS(basis_val[0]); i++) { - free(basis_val[0][i]); - free(basis_val[1][i]); - } - - return ret; -} - static int eval_psi_cylindric(const BDContext *bd, const double *rho, int nb_coords_rho, const double *z, int nb_coords_z, @@ -233,8 +115,8 @@ int bd_eval_psi(const BDContext *bd, case COORD_SYSTEM_CYLINDRICAL: ret = eval_psi_cylindric(bd, rho, nb_coords_rho, z, nb_coords_z, diff_order, psi, psi_stride); break; - case COORD_SYSTEM_RADIAL: ret = eval_psi_radial(bd, rho, nb_coords_rho, z, nb_coords_z, - diff_order, psi, psi_stride); + case COORD_SYSTEM_RADIAL: ret = bdi_eval_psi_radial(bd, rho, nb_coords_rho, z, nb_coords_z, + diff_order, psi, psi_stride); break; } |