summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAnton Khirnov <anton@khirnov.net>2020-01-02 14:13:23 +0100
committerAnton Khirnov <anton@khirnov.net>2020-01-20 14:45:07 +0100
commit004ea10198960ce0f7395975b027fd112edf960c (patch)
treec2cfd37aacbc7b84a4d1e19dad7f545d05832eae
parent3992a96fcd8f4bd63365c1c141184b3370881de6 (diff)
-rw-r--r--ell_grid_solve.c260
-rw-r--r--ell_grid_solve.h2
-rw-r--r--mg2d.c114
-rw-r--r--mg2d.h7
-rw-r--r--step_control.h283
5 files changed, 483 insertions, 183 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]);
diff --git a/ell_grid_solve.h b/ell_grid_solve.h
index 095c591..27c2848 100644
--- a/ell_grid_solve.h
+++ b/ell_grid_solve.h
@@ -88,6 +88,8 @@ typedef struct EGSRelaxContext {
*/
double relax_multiplier;
+ int adaptive_step;
+
Timer timer_correct;
} EGSRelaxContext;
diff --git a/mg2d.c b/mg2d.c
index 06e7b38..7d01761 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -37,6 +37,7 @@
#include "mg2d_boundary_internal.h"
#include "mg2d_constants.h"
#include "ndarray.h"
+#include "step_control.h"
#include "timer.h"
#include "transfer.h"
@@ -79,6 +80,8 @@ typedef struct MG2DLevel {
struct MG2DLevel *child;
+ StepControl sc;
+
/* timings */
int64_t count_cycles;
Timer timer_solve;
@@ -228,6 +231,7 @@ static int mg_prolong_u(MG2DContext *ctx, MG2DLevel *l)
return 0;
}
+#define SUBGRID_DIVERGE 1
static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level, int allow_exact)
{
enum EGSType solve_type = (allow_exact && level->dg->nb_components == 1 &&
@@ -264,7 +268,7 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level, int allow_exact)
goto finish;
}
- for (int i = 0; i < ctx->nb_cycles; i++) {
+ for (int cycle = 0; cycle < ctx->nb_cycles; cycle++) {
double res_prev;
/* pre-restrict relaxation */
@@ -277,17 +281,19 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level, int allow_exact)
}
if (level->depth < ctx->priv->dh.nb_levels - 1) {
- /* restrict the residual as to the coarser-level rhs */
- ret = mg_restrict_rhs(ctx, level);
- if (ret < 0)
- return ret;
-
- /* solve on the coarser level */
- if (level->child) {
- ret = mg_solve_subgrid(ctx, level->child, 1);
+ do {
+ /* restrict the residual as to the coarser-level rhs */
+ ret = mg_restrict_rhs(ctx, level);
if (ret < 0)
return ret;
- }
+
+ /* solve on the coarser level */
+ if (level->child) {
+ ret = mg_solve_subgrid(ctx, level->child, 1);
+ if (ret < 0)
+ return ret;
+ }
+ } while (ret == SUBGRID_DIVERGE);
/* prolongate the coarser-level correction */
ret = mg_prolong_u(ctx, level);
@@ -313,10 +319,59 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level, int allow_exact)
}
/* post-correct relaxation */
- for (int j = 0; j < ctx->nb_relax_post; j++) {
- ret = mg_solve(ctx, level, solve_type, "post-relax", 0);
- if (ret < 0)
- return ret;
+ {
+ double conv_history[4][2] = { { 0. } };
+ int conv_hist_next = 0;
+
+ double min_fact = DBL_MAX;
+ int min_fact_idx = -1;
+
+ for (int j = 0; j < ctx->nb_relax_post; j++) {
+ conv_history[conv_hist_next][0] = level->solver->residual_max;
+
+ ret = mg_solve(ctx, level, solve_type, "post-relax", 0);
+ if (ret < 0)
+ return ret;
+
+ conv_history[conv_hist_next][1] = level->solver->residual_max;
+ conv_hist_next = (conv_hist_next + 1) % ARRAY_ELEMS(conv_history);
+ }
+
+ // find the minimum convergence factor in the history window
+ for (int i = 0; i < ARRAY_ELEMS(conv_history); i++) {
+ double *h = conv_history[i], f;
+ if (h[0] == 0.0)
+ continue;
+ f = h[0] / h[1];
+ if (f < min_fact) {
+ min_fact = f;
+ min_fact_idx = i;
+ }
+ }
+
+ // use the minimum convergence factor as input for step control
+ if (level->solver->residual_max > ctx->tol && min_fact_idx >= 0) {
+ double *h = conv_history[min_fact_idx];
+ double step;
+
+ ret = sc_step_confirm(&level->sc, level->solver->relax->relax_factor,
+ h[0], h[1]);
+ if (ret < 0)
+ return ret;
+
+ step = sc_step_get(&level->sc);
+ if (fabs(step - level->solver->relax->relax_factor) > 1e-6) {
+ fprintf(stderr, "%d step %g -> %g\n", level->depth,
+ level->solver->relax->relax_factor, step);
+ }
+ if (ret) {
+ fprintf(stderr, "%d diverge %g -> %g\n", level->depth,
+ level->solver->relax->relax_factor, step);
+ level->solver->relax->relax_factor = step;
+ return SUBGRID_DIVERGE;
+ }
+ level->solver->relax->relax_factor = step;
+ }
}
level->count_cycles++;
@@ -786,9 +841,9 @@ static int mg_levels_init(MG2DContext *ctx)
ret = mg2di_ndarray_copy(cur->solver->u, priv->u);
if (ret < 0)
return ret;
- mg2di_ndarray_free(&priv->u);
- ctx->u = cur->solver->u->data;
- ctx->u_stride = cur->solver->u->stride[0];
+ //mg2di_ndarray_free(&priv->u);
+ //ctx->u = cur->solver->u->data;
+ //ctx->u_stride = cur->solver->u->stride[0];
}
if (priv->rhs) {
@@ -865,10 +920,10 @@ static int mg_levels_init(MG2DContext *ctx)
if (ctx->fd_stencil == 1)
op_restrict = GRID_TRANSFER_FW_2;
else
- op_restrict = GRID_TRANSFER_FW_3_GENERIC;
+ op_restrict = GRID_TRANSFER_FW_1;
} else {
op_restrict = GRID_TRANSFER_LAGRANGE_1;
- op_restrict = GRID_TRANSFER_FW_3_GENERIC;
+ op_restrict = GRID_TRANSFER_FW_1_GENERIC;
}
ret = mg_setup_restrict(ctx, cur, op_restrict);
@@ -904,9 +959,12 @@ static int mg_levels_init(MG2DContext *ctx)
cur = priv->root;
while (cur) {
- cur->solver->relax->relax_factor = ctx->cfl_factor;
+ //cur->solver->relax->relax_factor = ctx->cfl_factor;
+ sc_init(&cur->sc, ctx->cfl_factor, ctx->cfl_factor * 1e-3);
+ cur->solver->relax->relax_factor = sc_step_get(&cur->sc);
cur->solver->relax->relax_multiplier = 1.0 / (diff2_max + cur->solver->step[0] * cur->solver->step[1] *
diff0_max / 8.0);
+ cur->solver->relax->adaptive_step = 0;//ctx->adaptive_step;
if (cur->depth > 0) {
ret = diff_coeffs_fixup(ctx, cur);
@@ -1058,6 +1116,7 @@ static int mg_dh_init(DomainHierarchy *dh, const DomainGeometry *root,
dh->stepsizes[depth][0] = step_root[0] * (root->domain_size[0] - 1) / (cur_dg->domain_size[0] - 1);
dh->stepsizes[depth][1] = step_root[1] * (root->domain_size[1] - 1) / (cur_dg->domain_size[1] - 1);
+ fprintf(stderr, "level %d %zux%zu\n", depth, cur_dg->domain_size[0], cur_dg->domain_size[1]);
dh->nb_levels++;
if (next_size <= 4)
@@ -1127,9 +1186,15 @@ int mg2d_solve(MG2DContext *ctx)
res_prev = res_orig;
for (int i = 0; i < ctx->maxiter; i++) {
- ret = mg_solve_subgrid(ctx, root, 1);
- if (ret < 0)
- goto fail;
+ do {
+ ret = mg_solve_subgrid(ctx, root, 1);
+ if (ret < 0)
+ goto fail;
+ if (ret == SUBGRID_DIVERGE) {
+ mg2di_ndarray_copy(s_root->u, priv->u);
+ mg2di_egs_init(s_root, 0);
+ }
+ } while (ret == SUBGRID_DIVERGE);
res_cur = s_root->residual_max;
@@ -1159,6 +1224,7 @@ int mg2d_solve(MG2DContext *ctx)
fail:
mg2di_log(&priv->logger, MG2D_LOG_ERROR, "The solver failed to converge\n");
finish:
+ ctx->residual_max = s_root->residual_max;
mg2di_timer_stop(&priv->timer_solve);
return ret;
}
@@ -1291,6 +1357,8 @@ static MG2DLevel *mg_level_alloc(const DomainGeometry *dg,
mg2di_timer_init(&level->timer_reinit);
mg2di_timer_init(&level->timer_mpi_sync);
+ sc_reset(&level->sc);
+
return level;
fail:
mg_level_free(&level);
diff --git a/mg2d.h b/mg2d.h
index f7897f6..0c12f22 100644
--- a/mg2d.h
+++ b/mg2d.h
@@ -237,6 +237,13 @@ typedef struct MG2DContext {
* Set by mg2d_solver_alloc[_mpi]().
*/
size_t local_size[2];
+
+ double residual_max;
+
+ /**
+ * Use adaptive stepping for relaxation.
+ */
+ int adaptive_step;
} MG2DContext;
/**
diff --git a/step_control.h b/step_control.h
new file mode 100644
index 0000000..187c80b
--- /dev/null
+++ b/step_control.h
@@ -0,0 +1,283 @@
+/*
+ * Copyright 2020 Anton Khirnov <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef MG2D_STEP_CONTROL_H
+#define MG2D_STEP_CONTROL_H
+
+#include <float.h>
+
+#include "common.h"
+
+#define SC_HIST_SIZE 16
+
+typedef struct StepControl {
+ /* estimated value for a good step */
+ double hint;
+ /* minimum step size, fail if this value is reached */
+ double step_min;
+
+ int steps_since_init;
+
+ double step[SC_HIST_SIZE];
+ double fact[SC_HIST_SIZE];
+
+ double conv_epsilon;
+
+ int hist_len;
+ int idx_bracket;
+ int idx_diverge;
+
+ int diverge_count;
+} StepControl;
+
+static void sci_hist_clear(StepControl *sc)
+{
+ for (int i = 0; i < SC_HIST_SIZE; i++) {
+ sc->step[i] = DBL_MAX;
+ sc->fact[i] = DBL_MAX;
+ }
+
+ sc->hist_len = 0;
+ sc->idx_bracket = -1;
+ sc->idx_diverge = -1;
+}
+
+static void sc_reset(StepControl *sc)
+{
+ sci_hist_clear(sc);
+ sc->hint = DBL_MAX;
+ sc->step_min = DBL_MAX;
+ sc->steps_since_init = 0;
+}
+
+static void sc_init(StepControl *sc, double hint, double step_min)
+{
+ sc->hint = hint;
+ sc->step_min = step_min;
+ sc->steps_since_init = 0;
+ sc->conv_epsilon = hint * 1e-3;
+}
+
+static double sc_step_get(StepControl *sc)
+{
+ // no history present, try the hint
+ if (!sc->hist_len)
+ return (sc->hint == DBL_MAX) ? 0.25 : sc->hint;
+
+ // don't have the bracket yet
+ if (sc->idx_bracket < 0) {
+ const double high_step = sc->step[sc->hist_len - 1];
+ const double high_fact = sc->fact[sc->hist_len - 1];
+
+ // don't know where divergence happens
+ // -> try a higher step
+ if (high_fact > 1.0)
+ return high_step * 1.2;
+
+ // convergence factor should decrease towards smaller timesteps
+ // so just try lower ones until we get a bracket
+ return MAX(sc->step_min, sc->step[0] * 0.8);
+ }
+
+ // TODO periodic random steps to try escaping local minima
+
+ // got a bracket, golden-section-search it until convergence
+ mg2di_assert(sc->hist_len - sc->idx_bracket >= 3);
+ mg2di_assert(sc->fact[sc->idx_bracket] < sc->fact[sc->idx_bracket + 1]);
+ mg2di_assert(sc->fact[sc->idx_bracket + 2] < sc->fact[sc->idx_bracket + 1]);
+ {
+ const double dist0 = sc->step[sc->idx_bracket + 1] - sc->step[sc->idx_bracket];
+ const double dist1 = sc->step[sc->idx_bracket + 2] - sc->step[sc->idx_bracket + 1];
+
+ // converged
+ if (dist0 < sc->conv_epsilon && dist1 < sc->conv_epsilon)
+ return sc->step[sc->idx_bracket + 1];
+
+ if (dist0 > dist1)
+ return sc->step[sc->idx_bracket] + 0.61803 * dist0;
+ else
+ return sc->step[sc->idx_bracket + 2] - 0.61803 * dist1;
+ }
+}
+
+static void sci_find_bracket(StepControl *sc)
+{
+ for (int i = sc->hist_len - 2; i >= 0; i--) {
+ if (sc->fact[i] < sc->fact[i + 1] &&
+ sc->fact[i + 2] < sc->fact[i + 1]) {
+ sc->idx_bracket = i;
+ return;
+ }
+ }
+}
+
+static int sc_step_confirm(StepControl *sc, const double step,
+ const double norm_old, const double norm_new)
+{
+ const double conv_fact = norm_old / norm_new;
+
+ // signal failure if we reached minimum step
+ if (step <= sc->step_min)
+ return -EINVAL;
+
+ // ignore the first relax step, it tends to be weird
+ sc->steps_since_init++;
+ //if (sc->steps_since_init < 2)
+ // return 0;
+
+ // divergence
+ if (conv_fact <= 1.0) {
+ //if (step < sc->hint / 2.)
+ // return 0;
+
+ // if the step is not larger then the largest known diverging step,
+ // add it into history
+ if (!sc->hist_len ||
+ sc->idx_diverge < 0 ||
+ sc->step[sc->idx_diverge] > step) {
+ int idx = sc->idx_diverge >= 0 ? sc->idx_diverge : sc->hist_len;
+
+ // history size overflow, shouldn't happen
+ if (idx >= SC_HIST_SIZE)
+ return -ENOMEM;
+
+ // largest known coverging step doesn't converge anymore,
+ // clear history
+ if (idx > 0 && sc->step[idx - 1] <= step) {
+ sci_hist_clear(sc);
+ idx = 0;
+ }
+
+ sc->step[idx] = step;
+ sc->fact[idx] = conv_fact;
+ sc->idx_diverge = idx;
+ sc->hist_len = idx + 1;
+
+ // check if we now have a bracket
+ if (sc->idx_bracket < 0 && sc->hist_len >= 3 &&
+ sc->fact[idx - 1] > conv_fact &&
+ sc->fact[idx - 1] > sc->fact[idx - 2]) {
+ sc->idx_bracket = idx - 2;
+ }
+
+ }
+
+ return 1;
+ }
+
+ // convergence
+
+ // history empty, add first element
+ if (!sc->hist_len) {
+ sc->step[0] = step;
+ sc->fact[0] = conv_fact;
+ sc->hist_len = 1;
+ return 0;
+ }
+
+ // previously diverging step is now converging
+ // drop all diverging steps from history
+ if (sc->idx_diverge >= 0 &&
+ step >= sc->step[sc->idx_diverge]) {
+ const int idx = sc->idx_diverge;
+
+ sc->idx_diverge = -1;
+ sc->step[idx] = step;
+ sc->fact[idx] = conv_fact;
+ sc->hist_len = idx + 1;
+
+ // check if we just broke bracketing
+ if (sc->idx_bracket >= 0 && sc->idx_bracket + 2 == idx &&
+ sc->fact[sc->idx_bracket + 1] <= conv_fact)
+ sc->idx_bracket = -1;
+
+ // could not have created a bracket where there was none before, since
+ // we increased the factor for topmost step
+
+ return 0;
+ }
+
+ // no bracket
+ if (sc->idx_bracket < 0) {
+ int idx_floor = -1;
+
+ for (int i = 0; i < sc->hist_len; i++) {
+ if (sc->step[i] < step)
+ idx_floor = i;
+ else
+ break;
+ }
+ if (idx_floor + 1 < SC_HIST_SIZE) {
+ memmove(sc->step + idx_floor + 2, sc->step + idx_floor + 1,
+ sizeof(*sc->step) * (SC_HIST_SIZE - (idx_floor + 2)));
+ memmove(sc->fact + idx_floor + 2, sc->fact + idx_floor + 1,
+ sizeof(*sc->fact) * (SC_HIST_SIZE - (idx_floor + 2)));
+
+ sc->step[idx_floor + 1] = step;
+ sc->fact[idx_floor + 1] = conv_fact;
+ sc->hist_len = MIN(SC_HIST_SIZE, sc->hist_len + 1);
+ if (sc->idx_diverge > idx_floor)
+ sc->idx_diverge = (sc->idx_diverge + 1 < SC_HIST_SIZE) ? sc->idx_diverge + 1 : -1;
+ }
+ sci_find_bracket(sc);
+
+ return 0;
+ }
+
+ // have bracket
+
+ // got a step outside of the bracket for some reason
+ if (step <= sc->step[sc->idx_bracket] ||
+ step >= sc->step[sc->idx_bracket + 2])
+ return 0;
+
+ if (step < sc->step[sc->idx_bracket + 1]) {
+ // step inside lower interval
+ if (conv_fact > sc->fact[sc->idx_bracket + 1]) {
+ // replace upper bound
+ sc->step[sc->idx_bracket + 2] = sc->step[sc->idx_bracket + 1];
+ sc->fact[sc->idx_bracket + 2] = sc->fact[sc->idx_bracket + 1];
+ if (sc->idx_diverge == sc->idx_bracket + 2)
+ sc->idx_diverge = -1;
+
+ sc->step[sc->idx_bracket + 1] = step;
+ sc->fact[sc->idx_bracket + 1] = conv_fact;
+ } else {
+ // replace lower bound
+ sc->step[sc->idx_bracket] = step;
+ sc->fact[sc->idx_bracket] = conv_fact;
+ }
+ } else {
+ // step inside upper interval
+ if (conv_fact > sc->fact[sc->idx_bracket + 1]) {
+ // replace lower bound
+ sc->step[sc->idx_bracket] = sc->step[sc->idx_bracket + 1];
+ sc->fact[sc->idx_bracket] = sc->fact[sc->idx_bracket + 1];
+ sc->step[sc->idx_bracket + 1] = step;
+ sc->fact[sc->idx_bracket + 1] = conv_fact;
+ } else {
+ sc->step[sc->idx_bracket + 2] = step;
+ sc->fact[sc->idx_bracket + 2] = conv_fact;
+ if (sc->idx_diverge == sc->idx_bracket + 2)
+ sc->idx_diverge = -1;
+ }
+ }
+
+ return 0;
+}
+
+#endif // MG2D_STEP_CONTROL_H