From 0bc654a449de9308beaad98ee00861338654a447 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 9 Feb 2018 14:55:29 +0100 Subject: Attempt to enforce regularity at origin. probably useless, committing half-done in case it's ever needed again --- nlsolve.c | 4 +-- nlsolve.h | 2 +- pssolve.c | 108 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- pssolve.h | 2 +- 4 files changed, 111 insertions(+), 5 deletions(-) diff --git a/nlsolve.c b/nlsolve.c index 5d691e0..bb7e083 100644 --- a/nlsolve.c +++ b/nlsolve.c @@ -97,7 +97,7 @@ struct NLSolvePriv { }; int tdi_nlsolve_solve(NLSolveContext *ctx, NLEqCallback eq_eval, - NLEqJacobianCallback eq_jac_eval, void *opaque, double *coeffs) + NLEqJacobianCallback eq_jac_eval, void *opaque, double *coeffs, double L) { NLSolvePriv *s = ctx->priv; int64_t start, totaltime_start; @@ -189,7 +189,7 @@ int tdi_nlsolve_solve(NLSolveContext *ctx, NLEqCallback eq_eval, } // solve for delta - ret = tdi_pssolve_solve(s->ps_ctx, s->eq_coeffs, s->rhs, s->delta); + ret = tdi_pssolve_solve(s->ps_ctx, s->eq_coeffs, s->rhs, s->delta, L); if (ret < 0) return ret; diff --git a/nlsolve.h b/nlsolve.h index 1d3fae5..0fa8e9e 100644 --- a/nlsolve.h +++ b/nlsolve.h @@ -101,7 +101,7 @@ void tdi_nlsolve_context_free(NLSolveContext **ctx); int tdi_nlsolve_solve(NLSolveContext *ctx, NLEqCallback eq_eval, NLEqJacobianCallback eq_jac_eval, - void *opaque, double *coeffs); + void *opaque, double *coeffs, double L); void tdi_nlsolve_print_stats(NLSolveContext *ctx); diff --git a/pssolve.c b/pssolve.c index f2288c2..8ba613c 100644 --- a/pssolve.c +++ b/pssolve.c @@ -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++; diff --git a/pssolve.h b/pssolve.h index 6d5d1c0..b4a5966 100644 --- a/pssolve.h +++ b/pssolve.h @@ -181,7 +181,7 @@ void tdi_pssolve_context_free(PSSolveContext **ctx); */ 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); int tdi_pssolve_diff_order(enum PSSolveDiffOrder order, unsigned int dir); -- cgit v1.2.3