From 004ea10198960ce0f7395975b027fd112edf960c Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 2 Jan 2020 14:13:23 +0100 Subject: tmp --- ell_grid_solve.c | 260 +++++++++++++++++++++---------------------------------- 1 file changed, 100 insertions(+), 160 deletions(-) (limited to 'ell_grid_solve.c') 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]); -- cgit v1.2.3