summaryrefslogtreecommitdiff
path: root/mg2d.c
diff options
context:
space:
mode:
Diffstat (limited to 'mg2d.c')
-rw-r--r--mg2d.c114
1 files changed, 91 insertions, 23 deletions
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);
generated by cgit v1.2.3 (git 2.39.1) at 2024-04-27 23:25:35 +0000