1 files changed, 100 insertions, 160 deletions

diff --git a/ell_grid_solve.c b/ell_grid_solve.c
index ecbf6dc..daa62e6 100644
--- a/ell_grid_solve.c
+++ b/ell_grid_solve.c
@@ -39,92 +39,9 @@
 #include "mg2d_constants.h"
 #include "ndarray.h"
 #include "residual_calc.h"
+#include "step_control.h"
 #include "timer.h"
 
-typedef struct StepControl {
-    /* estimated value for a good step */
-    double hint;
-    /* minimum step size, fail if this value is reached */
-    double step_min;
-
-    double cur_converge;
-    double cur_diverge;
-    double prev_converge;
-    double prev_diverge;
-
-    int64_t diverge_count;
-} StepControl;
-
-static void sc_reset(StepControl *sc)
-{
-    sc->hint          = DBL_MAX;
-    sc->cur_converge  = DBL_MAX;
-    sc->cur_diverge   = DBL_MAX;
-    sc->prev_converge = DBL_MAX;
-    sc->prev_diverge  = DBL_MAX;
-}
-
-static void sc_init(StepControl *sc, double hint, double step_min)
-{
-    sc->hint         = hint;
-    sc->step_min     = step_min;
-    sc->cur_converge = DBL_MAX;
-    sc->cur_diverge  = DBL_MAX;
-}
-
-static void sc_step_converge(StepControl *sc, double step, double conv_fact)
-{
-    sc->cur_converge  = step;
-    sc->prev_converge = step;
-    mg2di_assert(sc->cur_diverge == DBL_MAX || step < sc->cur_diverge);
-    if (sc->prev_diverge != DBL_MAX && step >= sc->prev_diverge)
-        sc->prev_diverge = DBL_MAX;
-}
-
-static void sc_step_diverge(StepControl *sc, double step)
-{
-    sc->cur_diverge  = step;
-    sc->prev_diverge = step;
-    if (sc->cur_converge != DBL_MAX && step <= sc->cur_converge)
-        sc->cur_converge = DBL_MAX;
-    if (sc->prev_converge != DBL_MAX && step <= sc->prev_converge)
-        sc->prev_converge = DBL_MAX;
-
-    sc->diverge_count++;
-}
-
-static double sc_step_get(StepControl *sc)
-{
-    double step = 0.25;
-    double guess_converge = sc->cur_converge == DBL_MAX ? sc->prev_converge : sc->cur_converge;
-    double guess_diverge  = sc->cur_diverge  == DBL_MAX ? sc->prev_diverge  : sc->cur_diverge;
-
-    if (guess_converge != DBL_MAX && guess_diverge != DBL_MAX) {
-        if (guess_diverge - guess_converge < 0.1 * guess_diverge)
-            step = guess_converge;
-        else
-            step = 0.5 * (guess_diverge + guess_converge);
-    } else if (guess_diverge != DBL_MAX) {
-        step = 0.8 * guess_diverge;
-    } else if (guess_converge != DBL_MAX) {
-        step = 1.1 * guess_converge;
-    } else if (sc->hint != DBL_MAX)
-        step = sc->hint;
-
-    if (sc->step_min != DBL_MAX && step < sc->step_min) {
-        fprintf(stderr, "%g %g %g %g %g\n", sc->cur_converge, sc->prev_converge, sc->cur_diverge, sc->prev_diverge, sc->hint);
-        mg2di_assert(0);
-        return -1.0;
-    }
-
-    return step;
-}
-
-static int sc_step_confirm(StepControl *sc, double conv_factor)
-{
-    return 0;
-}
-
 static const double relax_factors[FD_STENCIL_MAX] = {
     [0] = 1.0 / 5,
     [1] = 1.0 / 7,
@@ -500,18 +417,21 @@ static int solve_relax_step(EGSContext *ctx, int export_res)
 {
     EGSInternal *priv = ctx->priv;
     EGSRelaxContext *r = ctx->relax;
-    int ic = priv->iteration_cur;
+
+    const int it_cur   = priv->iteration_cur       % HIST_SIZE;
+    const int it_next  = (priv->iteration_cur + 1) % HIST_SIZE;
+    const int it_next1 = (priv->iteration_cur + 2) % HIST_SIZE;
 
     if (export_res) {
-        if (priv->u_iteration[(ic + 1) % HIST_SIZE] < ic + 1) {
-            if (priv->residual_iteration < ic) {
-                residual_calc(ctx, &priv->residual_max[ic % HIST_SIZE], NULL, priv->u_interior[ic % HIST_SIZE], 0);
-                priv->residual_iteration = ic;
+        if (priv->u_iteration[it_next] < priv->iteration_cur + 1) {
+            if (priv->residual_iteration < priv->iteration_cur) {
+                residual_calc(ctx, &priv->residual_max[it_cur], NULL, priv->u_interior[it_cur], 0);
+                priv->residual_iteration = priv->iteration_cur;
             }
-            mg2di_assert(priv->residual_iteration == ic);
+            mg2di_assert(priv->residual_iteration == priv->iteration_cur);
 
-            priv->step[ic % HIST_SIZE] = sc_step_get(&priv->sc);
-            if (priv->step[ic % HIST_SIZE] < 0) {
+            priv->step[it_cur] = sc_step_get(&priv->sc);
+            if (priv->step[it_cur] < 0) {
                 mg2di_log(&ctx->logger, MG2D_LOG_ERROR,
                           "Minimum step reached.");
                 return -EDOM;
@@ -525,8 +445,8 @@ static int solve_relax_step(EGSContext *ctx, int export_res)
 
             priv->reflect_skip = 0;
 
-            boundaries_apply(ctx, priv->u_interior[(ic + 1) % HIST_SIZE], 0);
-            priv->u_iteration[(ic + 1) % HIST_SIZE] = ic + 1;
+            boundaries_apply(ctx, priv->u_interior[it_next], 0);
+            priv->u_iteration[it_next] = priv->iteration_cur + 1;
 
             priv->steps_since_sync++;
             if (priv->steps_since_sync >= priv->steps_per_sync) {
@@ -534,80 +454,95 @@ static int solve_relax_step(EGSContext *ctx, int export_res)
                 priv->steps_since_sync = 0;
             }
         }
-        mg2di_assert(priv->u_iteration[(ic + 1) % HIST_SIZE] == ic + 1);
+        mg2di_assert(priv->u_iteration[it_next] == priv->iteration_cur + 1);
 
-        residual_calc(ctx, &priv->residual_max[(ic + 1) % HIST_SIZE], NULL,
-                      priv->u_interior[(ic + 1) % HIST_SIZE], 0);
+        residual_calc(ctx, &priv->residual_max[it_next], NULL,
+                      priv->u_interior[it_next], 0);
         boundaries_sync(ctx, ctx->residual);
-        priv->residual_iteration = ic + 1;
+        priv->residual_iteration = priv->iteration_cur + 1;
     } else {
-        double res_old, res_new;
-        int diverged;
-
-        mg2di_assert(priv->u_iteration[(ic + 1) % HIST_SIZE] == ic + 1);
-        res_old  = priv->residual_max[ic % HIST_SIZE];
-
-        priv->step[(ic + 1) % HIST_SIZE] = sc_step_get(&priv->sc);
-        if (priv->step[(ic + 1) % HIST_SIZE] < 0) {
-            mg2di_log(&ctx->logger, MG2D_LOG_ERROR,
-                      "Minimum step reached.");
-            return -EDOM;
+        double step_next;
+
+        mg2di_assert(priv->u_iteration[it_next] == priv->iteration_cur + 1);
+
+        /* Adaptive relaxation stepping - this is slightly tricky, since
+         * residual computation and computing the next iteration are
+         * merged into a single step.
+         *
+         * At this point, we have already computed the next solution iteration
+         * -- u[it_next], but don't have the residual norm for it, so we haven't
+         *  verified whether the computation converged. We now generate a trial
+         *  next timestep and use it to compute the next residual+norm and
+         *  next+1 solution.
+         */
+        step_next = r->adaptive_step ? sc_step_get(&priv->sc) :
+                    r->relax_factor * r->relax_multiplier;//priv->step[0];
+        priv->r.relax_factor = step_next;
+        //fprintf(stderr, "%zu step next %g\n", ctx->domain_size[0], step_next);
+
+        residual_calc(ctx, &priv->residual_max[it_next], priv->u_interior[it_next1],
+                      priv->u_interior[it_next], step_next);
+
+        if (r->adaptive_step) {
+            const double res_old  = priv->residual_max[it_cur];
+            int step_reject;
+
+            /* Now that we have the next residual norm, we can verify whether
+             * the u[cur] --(step[cur])--> u[next] step was converging and
+             * accept or reject it.
+             *
+             * If it's rejected, we must generate new step[cur] and redo the
+             * calculations starting from u[cur].
+             */
+            do {
+                double res_new = priv->residual_max[it_next];
+                double tmp;
+
+                step_reject = sc_step_confirm(&priv->sc, priv->step[it_cur], res_old, res_new);
+                if (step_reject < 0) {
+                    mg2di_log(&ctx->logger, MG2D_LOG_ERROR, "Error confirming step: %g", step_next);
+                    return step_reject;
+                } else if (!step_reject)
+                    break;
+
+                priv->step[it_cur] = sc_step_get(&priv->sc);
+                residual_calc(ctx, &tmp, priv->u_interior[it_next],
+                              priv->u_interior[it_cur], priv->step[it_cur]);
+
+                step_next = sc_step_get(&priv->sc);
+                fprintf(stderr, "%zu retry step next %g %g\n", ctx->domain_size[0], step_next, res_old / res_new);
+                residual_calc(ctx, &priv->residual_max[it_next], priv->u_interior[it_next1],
+                              priv->u_interior[it_next], step_next);
+                res_new  = priv->residual_max[it_next];
+            } while (step_reject);
+
+            priv->step[it_next] = step_next;
         }
 
-        residual_calc(ctx, &priv->residual_max[(ic + 1) % HIST_SIZE], priv->u_interior[(ic + 2) % HIST_SIZE],
-                      priv->u_interior[(ic + 1) % HIST_SIZE], priv->step[(ic + 1) % HIST_SIZE]);
-
-        res_new  = priv->residual_max[(ic + 1) % HIST_SIZE];
-        diverged = res_new > 1e-6 && res_new / res_old - 1.0 >= 1e-2 ;
-        while (diverged) {
-            double tmp, step;
-
-            fprintf(stderr, "%zu diverge step %g fact %g\n", ctx->domain_size[0],
-                    priv->step[(ic + 1) % HIST_SIZE], res_new / res_old);
-
-            sc_step_diverge(&priv->sc, priv->step[ic % HIST_SIZE]);
-
-            priv->step[ic % HIST_SIZE] = sc_step_get(&priv->sc);
-            if (priv->step[(ic) % HIST_SIZE] < 0) {
-                mg2di_log(&ctx->logger, MG2D_LOG_ERROR,
-                          "Minimum step reached.");
-                return -EDOM;
-            }
-            residual_calc(ctx, &tmp, priv->u_interior[(ic + 1) % HIST_SIZE],
-                          priv->u_interior[ic % HIST_SIZE], priv->step[ic % HIST_SIZE]);
-
-            priv->step[(ic + 1) % HIST_SIZE] = priv->step[ic % HIST_SIZE];
-            residual_calc(ctx, &priv->residual_max[(ic + 1) % HIST_SIZE], priv->u_interior[(ic + 2) % HIST_SIZE],
-                          priv->u_interior[(ic + 1) % HIST_SIZE], priv->step[(ic + 1) % HIST_SIZE]);
-            res_new  = priv->residual_max[(ic + 1) % HIST_SIZE];
-            diverged = res_new > 1e-6 && res_new >= res_old;
-        }
-
-        sc_step_converge(&priv->sc, priv->step[ic % HIST_SIZE], res_old / res_new);
-        boundaries_apply(ctx, priv->u_interior[(ic + 2) % HIST_SIZE], 0);
-        priv->u_iteration[(ic + 2) % HIST_SIZE] = ic + 2;
+        boundaries_apply(ctx, priv->u_interior[it_next1], 0);
+        priv->u_iteration[it_next1] = priv->iteration_cur + 2;
 
         priv->steps_since_sync++;
         if (priv->steps_since_sync >= priv->steps_per_sync) {
-            boundaries_sync(ctx, priv->u_interior[(ic + 2) % HIST_SIZE]);
+            boundaries_sync(ctx, priv->u_interior[it_next1]);
             priv->steps_since_sync = 0;
         }
 
-        {
-            static int i;
-            if (!(i & 1023)) {
-                fprintf(stderr, "%zu converge %g diverge %g guess %g diverges %ld\n", ctx->domain_size[0], priv->sc.cur_converge,
-                        priv->sc.prev_diverge, ctx->relax->relax_factor * ctx->relax->relax_multiplier, priv->sc.diverge_count);
-            }
-            i++;
-        }
+        //if (r->adaptive_step) {
+        //    static int i;
+        //    if (!(i & 1023)) {
+        //        fprintf(stderr, "%zu converge %g %g diverge %g guess %g diverges %ld\n", ctx->domain_size[0], priv->sc.cur_converge,
+        //                priv->sc.cur_converge / ctx->relax->relax_multiplier,
+        //                priv->sc.prev_diverge, ctx->relax->relax_factor * ctx->relax->relax_multiplier, priv->sc.diverge_count);
+        //    }
+        //    i++;
+        //}
     }
 
-    ic++;
-    ctx->u              = priv->u_interior  [ic % HIST_SIZE];
-    ctx->u_exterior     = priv->u_exterior  [ic % HIST_SIZE];
-    ctx->residual_max   = priv->residual_max[ic % HIST_SIZE];
-    priv->iteration_cur = ic;
+    priv->iteration_cur++;
+    ctx->u              = priv->u_interior  [(priv->iteration_cur) % HIST_SIZE];
+    ctx->u_exterior     = priv->u_exterior  [(priv->iteration_cur) % HIST_SIZE];
+    ctx->residual_max   = priv->residual_max[(priv->iteration_cur) % HIST_SIZE];
 
     return 0;
 }
@@ -1482,11 +1417,16 @@ finish:
         }
         priv->steps_since_sync = 0;
 
-        priv->step[0] = sc_step_get(&priv->sc);
-        if (priv->step[0] < 0) {
-            mg2di_log(&ctx->logger, MG2D_LOG_ERROR,
-                      "Minimum step reached. 0");
-            return -EDOM;
+        if (r->adaptive_step) {
+            priv->step[0] = sc_step_get(&priv->sc);
+            if (priv->step[0] < 0) {
+                mg2di_log(&ctx->logger, MG2D_LOG_ERROR,
+                          "Minimum step reached.\n");
+                return -EDOM;
+            }
+        } else {
+            for (int i = 0; i < ARRAY_ELEMS(priv->step); i++)
+                priv->step[i] = r->relax_factor * r->relax_multiplier;
         }
 
         residual_calc(ctx, &priv->residual_max[0], priv->u_interior[1], priv->u_interior[0], priv->step[0]);