diff options
Diffstat (limited to 'pssolve.c')
| -rw-r--r-- | pssolve.c | 108 |
1 files changed, 107 insertions, 1 deletions
@@ -164,7 +164,7 @@ static int lu_invert(TDLogger *logger, const int N, double *mat, double *rhs, in int tdi_pssolve_solve(PSSolveContext *ctx, const double *(**eq_coeffs)[PSSOLVE_DIFF_ORDER_NB], - const double *rhs, double *coeffs) + /*const*/ double *rhs, double *coeffs, double L) { PSSolvePriv *s = ctx->priv; double rhs_max; @@ -193,6 +193,112 @@ int tdi_pssolve_solve(PSSolveContext *ctx, } } +#if TD_POLAR && 0 + // regularity conditions at origin for the polar basis + // FIXME: this does not properly belong here, should be handled more + // generically elsewhere + + // + // d^(2k) SB_{2n} (-1)^(n + k%2) · (2n + 1)^(2 · (k > 0)) · ∑_{l=0}^{k-1} (2n + 1)^(2l) sb_r0_derivative_coeffs[k-1][l] + // -------------- (x = 0) = ------------------------------------------------------------------------------------------------------ + // dx^(2k) L^(2k) + // + static const double sb_r0_derivative_coeffs[6][6] = { + // k = 1 (2nd derivative) + { 1.0 }, + // k = 2 (4th derivative) + { 8.0, 1.0 }, + // k = 3 (6th derivative) + { 184.0, 40.0, 1.0 }, + // k = 4 (8th derivative) + { 8448.0, 2464.0, 112.0, 1.0 }, + // k = 5 (10th derivative) + { 648576.0, 229760.0, 14448.0, 240.0, 1.0 }, + // k = 6 (12th derivative) + { 74972160.0, 30633856.0, 2393600.0, 55968.0, 440.0, 1.0 }, + }; + for (int i = 0; i < ctx->nb_equations; i++) { + PSEquationContext *eq_ctx = &s->eqs[i]; + + // FIXME assuming the basis is even cosines × odd SB + // i.e. for every m the corresponding part of the expansion is + // N-1 + // cos(2mθ) ∑ SB_{2n}(r) · c_{mn} + // n=0 + // and we need to enforce the corresponding radial part to have a root + // of order 2m at r=0, i.e. kth derivatives at r=0 are zero for k from 0 + // to 2m-1: + // N-1 / \ + // ∀k ∈ [0, 2m - 1]: ∑ |c_{mn} d^{k} SB_{2n}(r=0) / dr^k| = 0 + // n=0 \ / + // odd derivatives of SB_{2n} at r=0 are automatically zero, so this is + // just m nontrivial conditions + for (int m = 0; m < eq_ctx->nb_coeffs[1]; m++) { + for (int k = 0; k < MIN(m, 1); k++) { + // we replace the conditions from the outermost collocation points + // for the given m + //int idx_grid = m * eq_ctx->nb_colloc_points[0] + (eq_ctx->nb_colloc_points[0] - 1 - k); + int idx_grid = m * eq_ctx->nb_colloc_points[0] + k; + + double *mat = s->eqs[i].mat; + ptrdiff_t mat_stride = s->nb_coeffs; + + for (int j = 0; j < ctx->nb_equations; j++) { + // the regularity conditions are applied on the matrix "diagonal" - + // for ith variable in the ith equation block + if (i == j) { + for (int coeff_angular = 0; coeff_angular < eq_ctx->nb_coeffs[1]; coeff_angular++) { + if (coeff_angular == m) { + for (int coeff_radial = 0; coeff_radial < eq_ctx->nb_coeffs[0]; coeff_radial++) { + int idx_coeff = coeff_angular * eq_ctx->nb_coeffs[0] + coeff_radial; + double val; + + if (k > 0) { + double two_np1 = 2.0 * coeff_radial + 1.0; + double tmp = 0.0; + + for (int l = 0; l < k; l++) + tmp += pow(two_np1, 2.0 * l) * sb_r0_derivative_coeffs[k - 1][l]; + + val = SQR(two_np1) * tmp; + } else { + val = 1.0; + } + + val *= pow(-1.0, coeff_radial + (k % 2)) / pow(L, 2.0 * k); + + mat[idx_grid + mat_stride * idx_coeff] = val; + } + } else { + for (int coeff_radial = 0; coeff_radial < eq_ctx->nb_coeffs[0]; coeff_radial++) { + int idx_coeff = coeff_angular * eq_ctx->nb_coeffs[0] + coeff_radial; + mat[idx_grid + mat_stride * idx_coeff] = 0.0; + } + } + } + } else { + for (int idx_coeff = 0; idx_coeff < NB_COEFFS(&s->eqs[j]); idx_coeff++) + mat[idx_grid + mat_stride * idx_coeff] = 0.0; + } + + mat += NB_COEFFS(&s->eqs[j]) * mat_stride; + } + + + rhs[(s->eqs[i].mat - s->mat) + idx_grid] = 0.0; + } + } + } +#endif + +#if 0 + for (int i = 0; i < s->nb_coeffs; i++) { + for (int j = 0; j < s->nb_coeffs; j++) + fprintf(stdout, "%16.16g\t", s->mat[i + s->nb_coeffs * j]); + fprintf(stdout, "\n"); + } +#endif + ctx->construct_matrix_time += gettime() - start; ctx->construct_matrix_count++; |