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