[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:07 UTC 2009
The following commit has been merged in the upstream branch:
commit fc80f2fff753079da724d72eabbd977402c97ede
Author: Stephane Popinet <popinet at users.sf.net>
Date: Tue Mar 11 14:20:04 2008 +1100
Minimal Residual Smoothing implementation
This guarantees that the Poisson solver does not diverge (but does not
guarantee convergence). It can also accelerate convergence depending on the
problem. Coupled with other smaller changes included in this patch ('minlevel'
tuning in particular) this improves the robustness of the solver.
darcs-hash:20080311032004-d4795-ecaa4c108c57f9dd1230943a04fbb6aa1f047a4d.gz
diff --git a/src/event.c b/src/event.c
index a24d240..c234066 100644
--- a/src/event.c
+++ b/src/event.c
@@ -769,7 +769,7 @@ static void stream_from_vorticity (GfsDomain * domain,
gfs_multilevel_params_init (&par);
par.depth = gfs_domain_depth (domain);
while (norm.infty > tolerance && maxit) {
- gfs_poisson_cycle (domain, &par, stream, vorticity, dia, res);
+ gfs_poisson_cycle (domain, &par, stream, vorticity, dia, res, NULL, NULL);
norm = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res);
maxit--;
}
diff --git a/src/ocean.c b/src/ocean.c
index 009122e..af7c0a7 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -233,7 +233,7 @@ static void gfs_free_surface_pressure (GfsDomain * toplayer,
par->residual.second,
par->residual.infty);
#endif
- gfs_poisson_cycle (toplayer, par, p, fp.div, fp.dia, res1);
+ gfs_poisson_cycle (toplayer, par, p, fp.div, fp.dia, res1, NULL, NULL);
par->residual = gfs_domain_norm_residual (toplayer, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
par->niter++;
}
diff --git a/src/poisson.c b/src/poisson.c
index 2eb7cd3..b74e33e 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -634,6 +634,49 @@ static void get_from_below_2D (FttCell * cell, const GfsVariable * v)
GFS_VARIABLE (cell, v->i) = val;
}
+typedef struct {
+ GfsVariable * s, * r, * u, * v;
+ gdouble srs, rs2, beta;
+} MRSData;
+
+static void compute_beta (FttCell * cell, MRSData * data)
+{
+ gdouble rs = GFS_VALUE (cell, data->r) - GFS_VALUE (cell, data->s);
+ data->rs2 += rs*rs;
+ data->srs -= GFS_VALUE (cell, data->s)*rs;
+}
+
+static void update_sv (FttCell * cell, MRSData * data)
+{
+ GFS_VALUE (cell, data->s) += data->beta*(GFS_VALUE (cell, data->r) - GFS_VALUE (cell, data->s));
+ GFS_VALUE (cell, data->v) += data->beta*(GFS_VALUE (cell, data->u) - GFS_VALUE (cell, data->v));
+ GFS_VALUE (cell, data->r) = GFS_VALUE (cell, data->s);
+ GFS_VALUE (cell, data->u) = GFS_VALUE (cell, data->v);
+}
+
+/* See Jun Zhang, "Multigrid acceleration techniques and applications
+ to the numerical solutions of partial differential equations", PhD
+ Thesis, George Washington University, section 4.2 */
+static void minimal_residual_smoothing (GfsDomain * domain,
+ GfsVariable * v, GfsVariable * s,
+ GfsVariable * u, GfsVariable * res)
+{
+ MRSData data;
+ data.s = s;
+ data.r = res;
+ data.u = u;
+ data.v = v;
+ data.srs = data.rs2 = 0.;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) compute_beta, &data);
+ if (data.rs2 > 0.) {
+ data.beta = data.srs/data.rs2;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) update_sv, &data);
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, u);
+ }
+}
+
/**
* gfs_poisson_cycle:
* @domain: the domain on which to solve the Poisson equation.
@@ -659,7 +702,9 @@ void gfs_poisson_cycle (GfsDomain * domain,
GfsVariable * u,
GfsVariable * rhs,
GfsVariable * dia,
- GfsVariable * res)
+ GfsVariable * res,
+ GfsVariable * v,
+ GfsVariable * s)
{
guint n, l, nrelax, minlevel;
GfsVariable * dp;
@@ -719,6 +764,8 @@ void gfs_poisson_cycle (GfsDomain * domain,
gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, u);
/* compute new residual on leaf cells */
gfs_residual (domain, p->dimension, FTT_TRAVERSE_LEAFS, -1, u, rhs, dia, res);
+ if (v && s)
+ minimal_residual_smoothing (domain, v, s, u, res);
gts_object_destroy (GTS_OBJECT (dp));
}
diff --git a/src/poisson.h b/src/poisson.h
index 8c65153..965fc5b 100644
--- a/src/poisson.h
+++ b/src/poisson.h
@@ -73,7 +73,9 @@ void gfs_poisson_cycle (GfsDomain * domain,
GfsVariable * u,
GfsVariable * rhs,
GfsVariable * dia,
- GfsVariable * res);
+ GfsVariable * res,
+ GfsVariable * v,
+ GfsVariable * s);
void gfs_diffusion_coefficients (GfsDomain * domain,
GfsSourceDiffusion * d,
diff --git a/src/simulation.c b/src/simulation.c
index f70d13e..0ab1504 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -1451,10 +1451,20 @@ static void copy_res (FttCell * cell, gpointer * data)
GFS_VARIABLE (cell, res->i) = GFS_VARIABLE (cell, res1->i);
}
+typedef struct {
+ GfsVariable * u, * r, * s, * v;
+} MRSData;
+
+static void init_mrs (FttCell * cell, MRSData * d)
+{
+ GFS_VALUE (cell, d->v) = GFS_VALUE (cell, d->u);
+ GFS_VALUE (cell, d->s) = GFS_VALUE (cell, d->r);
+}
+
static void poisson_run (GfsSimulation * sim)
{
GfsDomain * domain = GFS_DOMAIN (sim);
- GfsVariable * dia, * div, * res = NULL, * res1, * p;
+ GfsVariable * dia, * div, * res = NULL, * res1, * p, * s, * v;
GfsMultilevelParams * par = &sim->approx_projection_params;
GSList * i;
@@ -1474,11 +1484,18 @@ static void poisson_run (GfsSimulation * sim)
gfs_poisson_coefficients (domain, NULL);
res1 = gfs_temporary_variable (domain);
dia = gfs_temporary_variable (domain);
+ s = gfs_temporary_variable (domain);
+ v = gfs_temporary_variable (domain);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) gfs_cell_reset, dia);
/* compute residual */
par->depth = gfs_domain_depth (domain);
gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res1);
+ /* initialize Minimal Residual Smoothing */
+ MRSData mrs;
+ mrs.u = p; mrs.r = res1; mrs.v = v; mrs.s = s;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) init_mrs, &mrs);
/* solve for pressure */
par->residual_before = par->residual =
gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res1);
@@ -1500,7 +1517,7 @@ static void poisson_run (GfsSimulation * sim)
gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);
gfs_domain_timer_start (domain, "poisson_cycle");
- gfs_poisson_cycle (domain, par, p, div, dia, res1);
+ gfs_poisson_cycle (domain, par, p, div, dia, res1, v, s);
par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res1);
gfs_domain_timer_stop (domain, "poisson_cycle");
@@ -1524,6 +1541,8 @@ static void poisson_run (GfsSimulation * sim)
gts_object_destroy (GTS_OBJECT (dia));
gts_object_destroy (GTS_OBJECT (div));
gts_object_destroy (GTS_OBJECT (res1));
+ gts_object_destroy (GTS_OBJECT (s));
+ gts_object_destroy (GTS_OBJECT (v));
}
static void poisson_class_init (GfsSimulationClass * klass)
diff --git a/src/timestep.c b/src/timestep.c
index 10b0a92..6c9d121 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -272,6 +272,16 @@ void gfs_update_gradients (GfsDomain * domain,
gfs_scale_gradients (domain, FTT_DIMENSION, g);
}
+typedef struct {
+ GfsVariable * u, * r, * s, * v;
+} MRSData;
+
+static void init_mrs (FttCell * cell, MRSData * d)
+{
+ GFS_VALUE (cell, d->v) = GFS_VALUE (cell, d->u);
+ GFS_VALUE (cell, d->s) = GFS_VALUE (cell, d->r);
+}
+
static void mac_projection (GfsDomain * domain,
GfsMultilevelParams * par,
GfsAdvectionParams * apar,
@@ -287,6 +297,8 @@ static void mac_projection (GfsDomain * domain,
GfsVariable * div = gfs_temporary_variable (domain);
GfsVariable * dia = gfs_temporary_variable (domain);
GfsVariable * res1 = res ? res : gfs_temporary_variable (domain);
+ GfsVariable * s = gfs_temporary_variable (domain);
+ GfsVariable * v = gfs_temporary_variable (domain);
/* Initialize face coefficients */
gfs_poisson_coefficients (domain, alpha);
@@ -319,9 +331,16 @@ static void mac_projection (GfsDomain * domain,
/* compute residual */
par->depth = gfs_domain_depth (domain);
gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res1);
+ /* initialize Minimal Residual Smoothing */
+ MRSData mrs;
+ mrs.u = p; mrs.r = res1; mrs.v = v; mrs.s = s;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) init_mrs, &mrs);
/* solve for pressure */
par->residual_before = par->residual =
gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
+ gdouble res_max_before = par->residual.infty;
+ guint minlevel = par->minlevel, erelax = par->erelax;
par->niter = 0;
while (par->niter < par->nitermin ||
(par->residual.infty > par->tolerance && par->niter < par->nitermax)) {
@@ -333,15 +352,28 @@ static void mac_projection (GfsDomain * domain,
par->residual.second,
par->residual.infty);
#endif
- gfs_poisson_cycle (domain, par, p, div, dia, res1);
+#if 1
+ gfs_poisson_cycle (domain, par, p, div, dia, res1, v, s);
+#else
+ gfs_poisson_cycle (domain, par, p, div, dia, res1, NULL, NULL);
+#endif
par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
+ if (par->residual.infty == res_max_before) /* convergence has stopped!! */
+ break;
+ if (par->residual.infty > res_max_before/1.1 && par->minlevel < par->depth)
+ par->minlevel++;
+ res_max_before = par->residual.infty;
par->niter++;
}
+ par->minlevel = minlevel;
+ par->erelax = erelax;
gts_object_destroy (GTS_OBJECT (div));
gts_object_destroy (GTS_OBJECT (dia));
if (!res)
gts_object_destroy (GTS_OBJECT (res1));
+ gts_object_destroy (GTS_OBJECT (s));
+ gts_object_destroy (GTS_OBJECT (v));
gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, apar->dt);
gfs_scale_gradients (domain, FTT_DIMENSION, g);
@@ -407,7 +439,8 @@ void gfs_mac_projection (GfsDomain * domain,
gfs_domain_timer_stop (domain, "mac_projection");
- g_assert (par->residual.infty <= par->tolerance);
+ if (par->residual.infty > par->tolerance)
+ g_warning ("MAC projection: max residual %g > %g", par->residual.infty, par->tolerance);
}
static void correct (FttCell * cell, gpointer * data)
@@ -513,7 +546,8 @@ void gfs_approximate_projection (GfsDomain * domain,
gfs_domain_timer_stop (domain, "approximate_projection");
- g_assert (par->residual.infty <= par->tolerance);
+ if (par->residual.infty > par->tolerance)
+ g_warning ("approx projection: max residual %g > %g", par->residual.infty, par->tolerance);
}
/**
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list