[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
s.popinet at niwa.co.nz
Fri May 15 02:52:30 UTC 2009
The following commit has been merged in the upstream branch:
commit 0940e6a925aeb541e33eb7f78528c07086686d48
Author: Stephane Popinet <s.popinet at niwa.co.nz>
Date: Fri Aug 26 11:28:13 2005 +1000
"Exponential" relaxation is controlable ("erelax" parameter)
darcs-hash:20050826012813-fbd8f-02a8af9fbd0f1646fb129248c33089777d09148a.gz
diff --git a/src/event.c b/src/event.c
index 516ceab..798730f 100644
--- a/src/event.c
+++ b/src/event.c
@@ -707,8 +707,9 @@ static void stream_from_vorticity (GfsDomain * domain,
gdouble tolerance)
{
GfsNorm norm;
- guint maxlevel, maxit = 100;
+ guint maxit = 100;
GfsVariable * res, * dia;
+ GfsMultilevelParams par;
g_return_if_fail (domain != NULL);
@@ -722,9 +723,12 @@ static void stream_from_vorticity (GfsDomain * domain,
res = gfs_temporary_variable (domain);
gfs_residual (domain, FTT_DIMENSION, FTT_TRAVERSE_LEAFS, -1, stream, vorticity, dia, res);
norm = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res);
- maxlevel = gfs_domain_depth (domain);
+ par.depth = gfs_domain_depth (domain);
+ par.minlevel = 0;
+ par.nrelax = 4;
+ par.dimension = FTT_DIMENSION;
while (norm.infty > tolerance && maxit) {
- gfs_poisson_cycle (domain, FTT_DIMENSION, 0, maxlevel, 4, stream, vorticity, dia, res);
+ gfs_poisson_cycle (domain, &par, stream, vorticity, dia, res);
norm = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res);
maxit--;
}
diff --git a/src/ocean.c b/src/ocean.c
index 06e147d..8f930c0 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -219,7 +219,6 @@ static void gfs_free_surface_pressure (GfsDomain * domain,
GfsVariable * res,
gdouble G)
{
- guint minlevel, maxlevel;
GfsDomain * toplayer;
FreeSurfaceParams fp;
GfsVariable * res1, * g[2];
@@ -263,14 +262,12 @@ static void gfs_free_surface_pressure (GfsDomain * domain,
(FttCellTraverseFunc) scale_divergence_helmoltz, &fp);
/* solve for pressure */
- minlevel = toplayer->rootlevel;
- if (par->minlevel > minlevel)
- minlevel = par->minlevel;
- maxlevel = gfs_domain_depth (toplayer);
+ par->depth = gfs_domain_depth (toplayer);
gfs_residual (toplayer, 2, FTT_TRAVERSE_LEAFS, -1, p, fp.div, fp.dia, res1);
par->residual_before = par->residual =
gfs_domain_norm_residual (toplayer, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
par->niter = 0;
+ par->dimension = 2;
while (par->residual.infty > par->tolerance && par->niter < par->nitermax) {
#if 0
fprintf (stderr, "%d bias: %g first: %g second: %g infty: %g\n",
@@ -280,7 +277,7 @@ static void gfs_free_surface_pressure (GfsDomain * domain,
par->residual.second,
par->residual.infty);
#endif
- gfs_poisson_cycle (toplayer, 2, minlevel, maxlevel, par->nrelax, p, fp.div, fp.dia, res1);
+ gfs_poisson_cycle (toplayer, par, p, fp.div, fp.dia, res1);
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 b521427..7296a12 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -22,6 +22,84 @@
#include "solid.h"
#include "source.h"
+/**
+ * gfs_multilevel_params_write:
+ * @par: the multilevel parameters.
+ * @fp: a file pointer.
+ *
+ * Writes in @fp a text representation of the multilevel parameters
+ * @par.
+ */
+void gfs_multilevel_params_write (GfsMultilevelParams * par, FILE * fp)
+{
+ g_return_if_fail (par != NULL);
+ g_return_if_fail (fp != NULL);
+
+ fprintf (fp,
+ "{\n"
+ " tolerance = %g\n"
+ " nrelax = %u\n"
+ " erelax = %u\n"
+ " minlevel = %u\n"
+ " nitermax = %u\n"
+ "}",
+ par->tolerance,
+ par->nrelax,
+ par->erelax,
+ par->minlevel,
+ par->nitermax);
+}
+
+void gfs_multilevel_params_init (GfsMultilevelParams * par)
+{
+ g_return_if_fail (par != NULL);
+
+ par->tolerance = 1e-3;
+ par->nrelax = 4;
+ par->erelax = 2;
+ par->minlevel = 0;
+ par->nitermax = 100;
+
+ par->dimension = FTT_DIMENSION;
+}
+
+void gfs_multilevel_params_read (GfsMultilevelParams * par, GtsFile * fp)
+{
+ GtsFileVariable var[] = {
+ {GTS_DOUBLE, "tolerance", TRUE},
+ {GTS_UINT, "nrelax", TRUE},
+ {GTS_UINT, "erelax", TRUE},
+ {GTS_UINT, "minlevel", TRUE},
+ {GTS_UINT, "nitermax", TRUE},
+ {GTS_NONE}
+ };
+
+ g_return_if_fail (par != NULL);
+ g_return_if_fail (fp != NULL);
+
+ var[0].data = &par->tolerance;
+ var[1].data = &par->nrelax;
+ var[2].data = &par->erelax;
+ var[3].data = &par->minlevel;
+ var[4].data = &par->nitermax;
+
+ gfs_multilevel_params_init (par);
+ gts_file_assign_variables (fp, var);
+ if (fp->type == GTS_ERROR)
+ return;
+
+ if (par->tolerance <= 0.) {
+ gts_file_variable_error (fp, var, "tolerance",
+ "tolerance `%g' must be strictly positive",
+ par->tolerance);
+ return;
+ }
+ if (par->nrelax == 0)
+ gts_file_variable_error (fp, var, "nrelax", "nrelax must be non zero");
+ if (par->erelax == 0)
+ gts_file_variable_error (fp, var, "erelax", "erelax must be non zero");
+}
+
#define FACE_GRADIENT(x, y, z, u) (gfs_face_weighted_gradient (x, y, z, u))
//#define FACE_GRADIENT(x, y, z, u) (gfs_face_gradient_flux_centered (x, y, z, u))
@@ -490,10 +568,7 @@ static void get_from_above (FttCell * parent, GfsVariable * v)
/**
* gfs_poisson_cycle:
* @domain: the domain on which to solve the Poisson equation.
- * @d: number of dimensions (2 or 3).
- * @levelmin: the top level of the multigrid hierarchy.
- * @depth: the total depth of the domain.
- * @nrelax: the number of relaxations to apply at each level.
+ * @p: the #GfsMultilevelParams.
* @u: the variable to use as left-hand side.
* @rhs: the variable to use as right-hand side.
* @dia: the diagonal weight.
@@ -511,27 +586,26 @@ static void get_from_above (FttCell * parent, GfsVariable * v)
* of @res (i.e. the cell tree is ready for another iteration).
*/
void gfs_poisson_cycle (GfsDomain * domain,
- guint d,
- guint levelmin,
- guint depth,
- guint nrelax,
+ GfsMultilevelParams * p,
GfsVariable * u,
GfsVariable * rhs,
GfsVariable * dia,
GfsVariable * res)
{
- guint n, l;
+ guint n, l, nrelax, minlevel;
GfsVariable * dp;
gpointer data[2];
g_return_if_fail (domain != NULL);
- g_return_if_fail (d > 1 && d <= 3);
+ g_return_if_fail (p != NULL);
+ g_return_if_fail (p->dimension > 1 && p->dimension <= 3);
g_return_if_fail (u != NULL);
g_return_if_fail (rhs != NULL);
g_return_if_fail (dia != NULL);
g_return_if_fail (res != NULL);
dp = gfs_temporary_variable (domain);
+ minlevel = MAX (domain->rootlevel, p->minlevel);
/* compute residual on non-leafs cells */
gfs_domain_cell_traverse (domain,
@@ -540,20 +614,21 @@ void gfs_poisson_cycle (GfsDomain * domain,
/* relax top level */
gfs_domain_cell_traverse (domain,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, levelmin,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, minlevel,
(FttCellTraverseFunc) gfs_cell_reset, dp);
- for (l = levelmin; l < depth; l++)
- nrelax *= 2;
+ nrelax = p->nrelax;
+ for (l = minlevel; l < p->depth; l++)
+ nrelax *= p->erelax;
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);
+ minlevel, dp, u);
+ gfs_relax (domain, p->dimension, minlevel, dp, res, dia);
}
- nrelax /= 2;
+ nrelax /= p->erelax;
/* relax from top to bottom */
- for (l = levelmin + 1; l <= depth; l++, nrelax /= 2) {
+ for (l = minlevel + 1; l <= p->depth; l++, nrelax /= p->erelax) {
/* get initial guess from coarser grid */
gfs_domain_cell_traverse (domain,
FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_NON_LEAFS, l - 1,
@@ -562,7 +637,7 @@ void gfs_poisson_cycle (GfsDomain * domain,
gfs_domain_homogeneous_bc (domain,
FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
l, dp, u);
- gfs_relax (domain, d, l, dp, res, dia);
+ gfs_relax (domain, p->dimension, l, dp, res, dia);
}
}
/* correct on leaf cells */
@@ -572,7 +647,7 @@ void gfs_poisson_cycle (GfsDomain * domain,
(FttCellTraverseFunc) correct, data);
gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, u);
/* compute new residual on leaf cells */
- gfs_residual (domain, d, FTT_TRAVERSE_LEAFS, -1, u, rhs, dia, res);
+ gfs_residual (domain, p->dimension, FTT_TRAVERSE_LEAFS, -1, u, rhs, dia, res);
gts_object_destroy (GTS_OBJECT (dp));
}
diff --git a/src/poisson.h b/src/poisson.h
index b852377..243da37 100644
--- a/src/poisson.h
+++ b/src/poisson.h
@@ -29,6 +29,26 @@ extern "C" {
#include "domain.h"
+typedef struct _GfsMultilevelParams GfsMultilevelParams;
+
+struct _GfsMultilevelParams {
+ gdouble tolerance;
+ guint nrelax, erelax;
+ guint minlevel;
+ guint nitermax;
+
+ guint dimension;
+ guint niter;
+ guint depth;
+ GfsNorm residual_before, residual;
+};
+
+void gfs_multilevel_params_init (GfsMultilevelParams * par);
+void gfs_multilevel_params_write (GfsMultilevelParams * par,
+ FILE * fp);
+void gfs_multilevel_params_read (GfsMultilevelParams * par,
+ GtsFile * fp);
+
void gfs_relax (GfsDomain * domain,
guint d,
gint max_depth,
@@ -47,10 +67,7 @@ void gfs_poisson_coefficients (GfsDomain * domain,
GfsVariable * c,
gdouble rho);
void gfs_poisson_cycle (GfsDomain * domain,
- guint d,
- guint levelmin,
- guint depth,
- guint nrelax,
+ GfsMultilevelParams * p,
GfsVariable * u,
GfsVariable * rhs,
GfsVariable * dia,
diff --git a/src/simulation.c b/src/simulation.c
index 9620659..d42ed0c 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -1331,7 +1331,6 @@ static void poisson_run (GfsSimulation * sim)
{
GfsDomain * domain;
GfsVariable * dia, * div, * res = NULL, * res1, * p;
- guint minlevel, maxlevel;
GfsMultilevelParams * par = &sim->approx_projection_params;
GSList * i;
@@ -1359,10 +1358,7 @@ static void poisson_run (GfsSimulation * sim)
dia = gfs_temporary_variable (domain);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) gfs_cell_reset, dia);
- minlevel = domain->rootlevel;
- if (par->minlevel > minlevel)
- minlevel = par->minlevel;
- maxlevel = gfs_domain_depth (domain);
+ par->depth = gfs_domain_depth (domain);
gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res1);
par->residual_before = par->residual =
gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res1);
@@ -1394,7 +1390,7 @@ static void poisson_run (GfsSimulation * 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);
+ gfs_poisson_cycle (domain, par, p, div, dia, res1);
par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res1);
gfs_domain_timer_stop (domain, "poisson_cycle");
diff --git a/src/timestep.c b/src/timestep.c
index 7757d45..78c62a3 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -24,77 +24,6 @@
#include "source.h"
#include "solid.h"
-/**
- * gfs_multilevel_params_write:
- * @par: the multilevel parameters.
- * @fp: a file pointer.
- *
- * Writes in @fp a text representation of the multilevel parameters
- * @par.
- */
-void gfs_multilevel_params_write (GfsMultilevelParams * par, FILE * fp)
-{
- g_return_if_fail (par != NULL);
- g_return_if_fail (fp != NULL);
-
- fprintf (fp,
- "{\n"
- " tolerance = %g\n"
- " nrelax = %u\n"
- " minlevel = %u\n"
- " nitermax = %u\n"
- "}",
- par->tolerance,
- par->nrelax,
- par->minlevel,
- par->nitermax);
-}
-
-void gfs_multilevel_params_init (GfsMultilevelParams * par)
-{
- g_return_if_fail (par != NULL);
-
- par->tolerance = 1e-3;
- par->nrelax = 4;
- par->minlevel = 0;
- par->nitermax = 100;
-
- par->dimension = FTT_DIMENSION;
-}
-
-void gfs_multilevel_params_read (GfsMultilevelParams * par, GtsFile * fp)
-{
- GtsFileVariable var[] = {
- {GTS_DOUBLE, "tolerance", TRUE},
- {GTS_UINT, "nrelax", TRUE},
- {GTS_UINT, "minlevel", TRUE},
- {GTS_UINT, "nitermax", TRUE},
- {GTS_NONE}
- };
-
- g_return_if_fail (par != NULL);
- g_return_if_fail (fp != NULL);
-
- var[0].data = &par->tolerance;
- var[1].data = &par->nrelax;
- var[2].data = &par->minlevel;
- var[3].data = &par->nitermax;
-
- gfs_multilevel_params_init (par);
- gts_file_assign_variables (fp, var);
- if (fp->type == GTS_ERROR)
- return;
-
- if (par->tolerance <= 0.) {
- gts_file_variable_error (fp, var, "tolerance",
- "tolerance `%g' must be strictly positive",
- par->tolerance);
- return;
- }
- if (par->nrelax == 0)
- gts_file_variable_error (fp, var, "nrelax", "nrelax must be non zero");
-}
-
static void reset_gradients (FttCell * cell, gpointer * data)
{
GfsVariable ** g = data[0];
@@ -262,7 +191,6 @@ void gfs_mac_projection (GfsDomain * domain,
GfsVariable * p,
GfsVariable ** g)
{
- guint minlevel, maxlevel;
gdouble dt;
gpointer data[2];
GfsVariable * div, * dia, * res;
@@ -310,17 +238,14 @@ void gfs_mac_projection (GfsDomain * domain,
#endif
/* solve for pressure */
- minlevel = domain->rootlevel;
- if (par->minlevel > minlevel)
- minlevel = par->minlevel;
- maxlevel = gfs_domain_depth (domain);
+ par->depth = gfs_domain_depth (domain);
gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res);
par->residual_before = par->residual =
gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res);
par->niter = 0;
while (par->residual.infty > par->tolerance &&
par->niter < par->nitermax) {
- gfs_poisson_cycle (domain, par->dimension, minlevel, maxlevel, par->nrelax, p, div, dia, res);
+ gfs_poisson_cycle (domain, par, p, div, dia, res);
par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res);
par->niter++;
}
@@ -425,7 +350,6 @@ void gfs_approximate_projection (GfsDomain * domain,
GfsVariable * p,
GfsVariable * res)
{
- guint minlevel, maxlevel;
gpointer data[2];
GfsVariable * dia, * div, * g[FTT_DIMENSION], * res1;
@@ -476,10 +400,7 @@ void gfs_approximate_projection (GfsDomain * domain,
#endif
/* solve for pressure */
- minlevel = domain->rootlevel;
- if (par->minlevel > minlevel)
- minlevel = par->minlevel;
- maxlevel = gfs_domain_depth (domain);
+ par->depth = gfs_domain_depth (domain);
gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res1);
par->residual_before = par->residual =
gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
@@ -494,7 +415,7 @@ void gfs_approximate_projection (GfsDomain * domain,
par->residual.second,
par->residual.infty);
#endif
- gfs_poisson_cycle (domain, par->dimension, minlevel, maxlevel, par->nrelax, p, div, dia, res1);
+ gfs_poisson_cycle (domain, par, p, div, dia, res1);
par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
par->niter++;
}
diff --git a/src/timestep.h b/src/timestep.h
index 06f1177..756d68e 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -26,27 +26,8 @@ extern "C" {
#include "advection.h"
#include "poisson.h"
-
-typedef struct _GfsMultilevelParams GfsMultilevelParams;
-
-struct _GfsMultilevelParams {
- gdouble tolerance;
- guint nrelax;
- guint minlevel;
- guint nitermax;
-
- guint dimension;
- guint niter;
- GfsNorm residual_before, residual;
-};
-
#include "variable.h"
-void gfs_multilevel_params_init (GfsMultilevelParams * par);
-void gfs_multilevel_params_write (GfsMultilevelParams * par,
- FILE * fp);
-void gfs_multilevel_params_read (GfsMultilevelParams * par,
- GtsFile * fp);
void gfs_correct_normal_velocities (GfsDomain * domain,
guint dimension,
GfsVariable * p,
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list