[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
s.popinet at niwa.co.nz
Fri May 15 02:52:25 UTC 2009
The following commit has been merged in the upstream branch:
commit fdc99efe4e5fbc4d01db80e38456db5a91d1818a
Author: Stephane Popinet <s.popinet at niwa.co.nz>
Date: Tue Aug 9 14:29:26 2005 +1000
Improvements to the Poisson solver
After a full review of the Poisson solver: convergence tests with pre
and post relaxations, FMG implementation etc... the following
relatively minor changes were made: (1) the prolongation operator uses
"second-order" gradient-based interpolation rather than straight
injection, (2) the number of relaxations increases exponentially for
coarser levels. None of the more complex changes seemed to improve
convergence.
darcs-hash:20050809042926-fbd8f-e1daabedf58a7f582664caddc8457d7f4cc2c061.gz
diff --git a/src/fluid.c b/src/fluid.c
index 5693026..fcb2194 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -1228,25 +1228,6 @@ gdouble gfs_cell_dirichlet_value (FttCell * cell,
}
/**
- * gfs_get_from_above:
- * @cell: a #FttCell.
- * @v: a #GfsVariable to "get from above".
- *
- * Sets the value of the variable @v of @cell to the value of this
- * variable in its parent cell.
- *
- * This function fails if @cell is the root of the cell tree.
- */
-void gfs_get_from_above (FttCell * cell, const GfsVariable * v)
-{
- g_return_if_fail (cell != NULL);
- g_return_if_fail (!FTT_CELL_IS_ROOT (cell));
- g_return_if_fail (v != NULL);
-
- GFS_VARIABLE (cell, v->i) = GFS_VARIABLE (ftt_cell_parent (cell), v->i);
-}
-
-/**
* gfs_get_from_below_intensive:
* @cell: a #FttCell.
* @v: a #GfsVariable to "get from below".
diff --git a/src/fluid.h b/src/fluid.h
index 924528d..25fc1a8 100644
--- a/src/fluid.h
+++ b/src/fluid.h
@@ -93,8 +93,6 @@ typedef enum {
void gfs_cell_cleanup (FttCell * cell);
void gfs_cell_reset (FttCell * cell,
GfsVariable * v);
-void gfs_get_from_above (FttCell * cell,
- const GfsVariable * v);
void gfs_get_from_below_intensive (FttCell * cell,
const GfsVariable * v);
void gfs_get_from_below_extensive (FttCell * cell,
diff --git a/src/ftt.c b/src/ftt.c
index 98b234d..cccd3c3 100644
--- a/src/ftt.c
+++ b/src/ftt.c
@@ -1329,7 +1329,8 @@ static void cell_traverse_level_non_leafs (FttCell * cell,
* @data: user data to pass to @func.
*
* Traverses a cell tree starting at the given root #FttCell. Calls
- * the given function for each cell visited. */
+ * the given function for each cell visited.
+ */
void ftt_cell_traverse (FttCell * root,
FttTraverseType order,
FttTraverseFlags flags,
diff --git a/src/poisson.c b/src/poisson.c
index cb5c05e..b521427 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -241,7 +241,7 @@ static void residual_set (FttCell * cell, RelaxParams * p)
for (f.d = 0; f.d < FTT_NEIGHBORS; f.d++) {
f.neighbor = neighbor.c[f.d];
if (f.neighbor) {
- FACE_GRADIENT (&f, &ng, p->u, -1);
+ FACE_GRADIENT (&f, &ng, p->u, p->maxlevel);
g.a += ng.a;
g.b += ng.b;
}
@@ -264,7 +264,7 @@ static void residual_set2D (FttCell * cell, RelaxParams * p)
for (f.d = 0; f.d < FTT_NEIGHBORS_2D; f.d++) {
f.neighbor = neighbor.c[f.d];
if (f.neighbor) {
- FACE_GRADIENT (&f, &ng, p->u, -1);
+ FACE_GRADIENT (&f, &ng, p->u, p->maxlevel);
g.a += ng.a;
g.b += ng.b;
}
@@ -309,6 +309,7 @@ void gfs_residual (GfsDomain * domain,
p.rhs = rhs->i;
p.dia = dia->i;
p.res = res->i;
+ p.maxlevel = max_depth;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, flags, max_depth,
(FttCellTraverseFunc) (d == 2 ? residual_set2D : residual_set), &p);
}
@@ -447,6 +448,45 @@ static void correct (FttCell * cell, gpointer * data)
GFS_VARIABLE (cell, u->i) += GFS_VARIABLE (cell, dp->i);
}
+static void get_from_above (FttCell * parent, GfsVariable * v)
+{
+ guint level = ftt_cell_level (parent);
+ FttCellNeighbors n;
+ FttCellChildren child;
+ FttComponent c;
+ FttVector h;
+ guint i;
+
+ ftt_cell_neighbors (parent, &n);
+ for (c = 0; c < FTT_DIMENSION; c++) {
+ FttCellFace f;
+ GfsGradient g;
+ gdouble g1, g2;
+
+ f.cell = parent;
+ f.d = 2*c;
+ f.neighbor = n.c[f.d];
+ gfs_face_gradient (&f, &g, v->i, level);
+ g1 = g.b - g.a*GFS_VARIABLE (parent, v->i);
+ f.d = 2*c + 1;
+ f.neighbor = n.c[f.d];
+ gfs_face_gradient (&f, &g, v->i, level);
+ g2 = g.b - g.a*GFS_VARIABLE (parent, v->i);
+ (&h.x)[c] = (g1 - g2)/2.;
+ }
+
+ ftt_cell_children (parent, &child);
+ for (i = 0; i < FTT_CELLS; i++)
+ if (child.c[i]) {
+ FttVector p;
+
+ GFS_VARIABLE (child.c[i], v->i) = GFS_VARIABLE (parent, v->i);
+ ftt_cell_relative_pos (child.c[i], &p);
+ for (c = 0; c < FTT_DIMENSION; c++)
+ GFS_VARIABLE (child.c[i], v->i) += (&p.x)[c]*(&h.x)[c];
+ }
+}
+
/**
* gfs_poisson_cycle:
* @domain: the domain on which to solve the Poisson equation.
@@ -502,19 +542,22 @@ void gfs_poisson_cycle (GfsDomain * domain,
gfs_domain_cell_traverse (domain,
FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, levelmin,
(FttCellTraverseFunc) gfs_cell_reset, dp);
- for (n = 0; n < 10*nrelax; n++) {
+ for (l = levelmin; l < depth; l++)
+ nrelax *= 2;
+ for (n = 0; n < nrelax; n++) {
gfs_domain_homogeneous_bc (domain,
FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
levelmin, dp, u);
gfs_relax (domain, d, levelmin, dp, res, dia);
}
+ nrelax /= 2;
/* relax from top to bottom */
- for (l = levelmin + 1; l <= depth; l++) {
+ for (l = levelmin + 1; l <= depth; l++, nrelax /= 2) {
/* get initial guess from coarser grid */
- gfs_domain_cell_traverse (domain,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, l,
- (FttCellTraverseFunc) gfs_get_from_above, dp);
+ gfs_domain_cell_traverse (domain,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_NON_LEAFS, l - 1,
+ (FttCellTraverseFunc) get_from_above, dp);
for (n = 0; n < nrelax; n++) {
gfs_domain_homogeneous_bc (domain,
FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
@@ -900,9 +943,10 @@ void gfs_diffusion_cycle (GfsDomain * domain,
/* relax from top to bottom */
for (p.maxlevel = levelmin + 1; p.maxlevel <= depth; p.maxlevel++) {
/* get initial guess from coarser grid */
- gfs_domain_cell_traverse (domain,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, p.maxlevel,
- (FttCellTraverseFunc) gfs_get_from_above, dp);
+ gfs_domain_cell_traverse (domain,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_NON_LEAFS,
+ p.maxlevel - 1,
+ (FttCellTraverseFunc) get_from_above, dp);
for (n = 0; n < nrelax; n++) {
gfs_domain_homogeneous_bc (domain,
FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
diff --git a/src/simulation.c b/src/simulation.c
index 25a8254..9620659 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -1321,6 +1321,12 @@ static void correct_div (GfsDomain * domain, GfsVariable * divu, GfsVariable * d
(FttCellTraverseFunc) add_ddiv, data);
}
+static void copy_res (FttCell * cell, gpointer * data)
+{
+ GfsVariable * res = data[0], * res1 = data[1];
+ GFS_VARIABLE (cell, res->i) = GFS_VARIABLE (cell, res1->i);
+}
+
static void poisson_run (GfsSimulation * sim)
{
GfsDomain * domain;
@@ -1349,7 +1355,7 @@ static void poisson_run (GfsSimulation * sim)
div = gfs_temporary_variable (domain);
correct_div (domain, gfs_variable_from_name (domain->variables, "Div"), div);
gfs_poisson_coefficients (domain, NULL, 0.);
- res1 = res ? res : gfs_temporary_variable (domain);
+ res1 = gfs_temporary_variable (domain);
dia = gfs_temporary_variable (domain);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) gfs_cell_reset, dia);
@@ -1363,23 +1369,34 @@ static void poisson_run (GfsSimulation * sim)
par->niter = 0;
while (sim->time.t < sim->time.end &&
sim->time.i < sim->time.iend &&
+ sim->time.i < par->nitermax &&
par->residual.infty > par->tolerance) {
gdouble tstart;
+ if (res) {
+ gpointer data[2];
+ data[0] = res;
+ data[1] = res1;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) copy_res, data);
+ }
+
i = domain->variables;
while (i) {
gfs_event_do (GFS_EVENT (i->data), sim);
i = i->next;
}
gfs_domain_cell_traverse (domain,
- FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
- (FttCellTraverseFunc) gfs_cell_coarse_init, domain);
+ FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
+ (FttCellTraverseFunc) gfs_cell_coarse_init, domain);
gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);
tstart = gfs_clock_elapsed (domain->timer);
+ gfs_domain_timer_start (domain, "poisson_cycle");
gfs_poisson_cycle (domain, par->dimension, minlevel, maxlevel, par->nrelax, p, div, dia, res1);
par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res1);
+ gfs_domain_timer_stop (domain, "poisson_cycle");
gfs_simulation_adapt (sim, NULL);
@@ -1397,8 +1414,7 @@ static void poisson_run (GfsSimulation * sim)
(GtsFunc) gts_object_destroy, NULL);
gts_object_destroy (GTS_OBJECT (dia));
gts_object_destroy (GTS_OBJECT (div));
- if (!res)
- gts_object_destroy (GTS_OBJECT (res1));
+ gts_object_destroy (GTS_OBJECT (res1));
}
static void poisson_class_init (GfsSimulationClass * klass)
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list