Attempt to enforce regularity at origin.origin_regularity

probably useless, committing half-done in case it's ever needed again

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);