[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
s.popinet at niwa.co.nz
Fri May 15 02:52:36 UTC 2009
The following commit has been merged in the upstream branch:
commit cd6d1b364fceac8381a464bab9da701a8c0e25c8
Author: Stephane Popinet <s.popinet at niwa.co.nz>
Date: Wed Sep 14 15:32:59 2005 +1000
Weighted centered pressure gradient is now an option used only by the ocean models
darcs-hash:20050914053259-fbd8f-e8469ebae6773cfb0f3cedce193f1b4fa699c0f3.gz
diff --git a/src/ocean.c b/src/ocean.c
index 606285a..c9dc5ce 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -27,6 +27,143 @@
#include "vof.h"
#include "graphic.h"
+static void reset_gradients (FttCell * cell, gpointer * data)
+{
+ GfsVariable ** g = data[0];
+ guint * dimension = data[1];
+ FttComponent c;
+
+ for (c = 0; c < *dimension; c++)
+ GFS_VARIABLE (cell, g[c]->i) = 0.;
+}
+
+static void correct_normal_velocity_weighted (FttCellFace * face,
+ gpointer * data)
+{
+ GfsGradient g;
+ gdouble dp;
+ FttFaceType type;
+ GfsStateVector * s;
+ GfsVariable * p = data[0];
+ GfsVariable ** gv = data[1];
+ gdouble * dt = data[2];
+ FttComponent c;
+
+ if (GFS_FACE_FRACTION (face) == 0.)
+ return;
+
+ s = GFS_STATE (face->cell);
+ type = ftt_face_type (face);
+ c = face->d/2;
+
+ gfs_face_weighted_gradient (face, &g, p->i, -1);
+ dp = (g.b - g.a*GFS_VARIABLE (face->cell, p->i))/ftt_cell_size (face->cell);
+ if (!FTT_FACE_DIRECT (face))
+ dp = - dp;
+
+ if (s->solid && s->solid->s[face->d] > 0.)
+ dp /= s->solid->s[face->d];
+
+ GFS_FACE_NORMAL_VELOCITY_LEFT (face) -= dp*(*dt);
+ GFS_VARIABLE (face->cell, gv[c]->i) += dp*GFS_FACE_FRACTION_LEFT (face);
+
+ switch (type) {
+ case FTT_FINE_FINE:
+ GFS_FACE_NORMAL_VELOCITY_RIGHT (face) -= dp*(*dt);
+ GFS_VARIABLE (face->neighbor, gv[c]->i) += dp*GFS_FACE_FRACTION_RIGHT (face);
+ break;
+ case FTT_FINE_COARSE: {
+ dp *= GFS_FACE_FRACTION_LEFT (face)/(FTT_CELLS/2);
+ GFS_VARIABLE (face->neighbor, gv[c]->i) += dp;
+ g_assert (GFS_FACE_FRACTION_RIGHT (face) > 0.);
+ GFS_FACE_NORMAL_VELOCITY_RIGHT (face) -= dp/GFS_FACE_FRACTION_RIGHT (face)*(*dt);
+ break;
+ }
+ default:
+ g_assert_not_reached ();
+ }
+}
+
+static void scale_gradients_weighted (FttCell * cell, gpointer * data)
+{
+ GfsVariable ** g = data[0];
+ guint * dimension = data[1];
+ FttComponent c;
+
+ if (GFS_IS_MIXED (cell)) {
+ GfsSolidVector * s = GFS_STATE (cell)->solid;
+
+ for (c = 0; c < *dimension; c++) {
+ g_assert (s->s[2*c] + s->s[2*c + 1] > 0.);
+ GFS_VARIABLE (cell, g[c]->i) /= s->s[2*c] + s->s[2*c + 1];
+ }
+ }
+ else {
+ FttCellNeighbors n;
+
+ ftt_cell_neighbors (cell, &n);
+ for (c = 0; c < *dimension; c++) {
+ FttCell * c1 = n.c[2*c], * c2 = n.c[2*c + 1];
+
+ if (c1 && c2 && !GFS_CELL_IS_GRADIENT_BOUNDARY (c1) && !GFS_CELL_IS_GRADIENT_BOUNDARY (c2))
+ GFS_VARIABLE (cell, g[c]->i) /= 2.;
+ }
+ }
+}
+
+/**
+ * gfs_correct_normal_velocities_weighted:
+ * @domain: a #GfsDomain.
+ * @dimension: the number of dimensions (2 or 3).
+ * @p: the pressure field.
+ * @g: where to store the pressure gradient.
+ * @dt: the timestep.
+ * @weighted: whether to use fraction-weighting or not.
+ *
+ * Corrects the normal velocity field of @domain using @p and and @dt.
+ *
+ * Also allocates the @g variables and fills them with the centered gradient of @p.
+ */
+static void gfs_correct_normal_velocities_weighted (GfsDomain * domain,
+ guint dimension,
+ GfsVariable * p,
+ GfsVariable ** g,
+ gdouble dt,
+ gboolean weighted)
+{
+ if (!weighted)
+ gfs_correct_normal_velocities (domain, dimension, p, g, dt);
+ else {
+ gpointer data[3];
+ FttComponent c;
+
+ g_return_if_fail (domain != NULL);
+ g_return_if_fail (p != NULL);
+ g_return_if_fail (g != NULL);
+
+ for (c = 0; c < dimension; c++) {
+ g[c] = gfs_temporary_variable (domain);
+ gfs_variable_set_vector (g[c], c);
+ }
+ data[0] = g;
+ data[1] = &dimension;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) reset_gradients, data);
+ data[0] = p;
+ data[1] = g;
+ data[2] = &dt;
+ gfs_domain_face_traverse (domain, dimension == 2 ? FTT_XY : FTT_XYZ,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttFaceTraverseFunc) correct_normal_velocity_weighted, data);
+ data[0] = g;
+ data[1] = &dimension;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) scale_gradients_weighted, data);
+ for (c = 0; c < dimension; c++)
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, g[c]);
+ }
+}
+
/* GfsOcean: Object */
static void ocean_destroy (GtsObject * object)
@@ -292,7 +429,7 @@ static void gfs_free_surface_pressure (GfsDomain * domain,
(FttCellTraverseFunc) column_pressure, p);
gfs_poisson_coefficients (toplayer, NULL, 0.);
#endif /* 2D3 or 3D */
- gfs_correct_normal_velocities (domain, 2, p, g, apar->dt/2.);
+ gfs_correct_normal_velocities_weighted (domain, 2, p, g, apar->dt/2., par->weighted);
#if (!FTT_2D)
gfs_domain_cell_traverse_boundary (domain, FTT_BACK,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
@@ -372,7 +509,8 @@ static void ocean_run (GfsSimulation * sim)
gfs_domain_timer_start (domain, "correct_normal_velocities");
gfs_poisson_coefficients (domain, NULL, 0.);
- gfs_correct_normal_velocities (domain, 2, p, g, sim->advection_params.dt/2.);
+ gfs_correct_normal_velocities_weighted (domain, 2, p, g, sim->advection_params.dt/2.,
+ sim->approx_projection_params.weighted);
#if (!FTT_2D)
gfs_domain_cell_traverse_boundary (domain, FTT_BACK,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
@@ -408,11 +546,11 @@ static void ocean_run (GfsSimulation * sim)
g);
gfs_poisson_coefficients (domain, NULL, 0.);
- gfs_correct_normal_velocities (domain, 2, p, g, 0.);
+ gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., sim->approx_projection_params.weighted);
if (gfs_has_source_coriolis (domain)) {
gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt);
gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
- gfs_correct_normal_velocities (domain, 2, p, g, 0.);
+ gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., sim->approx_projection_params.weighted);
gfs_correct_centered_velocities (domain, 2, g, -sim->advection_params.dt/2.);
}
else
@@ -449,6 +587,7 @@ static void gfs_ocean_class_init (GfsSimulationClass * klass)
static void gfs_ocean_init (GfsOcean * object)
{
+ GFS_SIMULATION (object)->approx_projection_params.weighted = 1;
object->layer = g_ptr_array_new ();
new_layer (object);
}
@@ -992,19 +1131,20 @@ static void ocean1_run (GfsSimulation * sim)
gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_half_do, sim);
- gfs_correct_normal_velocities (domain, 2, p, g, 0.); // just there so that the call below
- // has sthg to free
+ gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., FALSE);
+ // just there so that the call below
+ // has sthg to free
gfs_centered_velocity_advection_diffusion (domain, 2,
&sim->advection_params,
&sim->diffusion_params,
g);
gfs_poisson_coefficients (domain, H, 0.);
- gfs_correct_normal_velocities (domain, 2, p, g, 0.);
+ gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., sim->approx_projection_params.weighted);
if (gfs_has_source_coriolis (domain)) {
gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt);
gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
- gfs_correct_normal_velocities (domain, 2, p, g, 0.);
+ gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., sim->approx_projection_params.weighted);
gfs_correct_centered_velocities (domain, 2, g, -sim->advection_params.dt/2.);
}
else
diff --git a/src/poisson.c b/src/poisson.c
index b71f777..3694aa4 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -42,12 +42,14 @@ void gfs_multilevel_params_write (GfsMultilevelParams * par, FILE * fp)
" erelax = %u\n"
" minlevel = %u\n"
" nitermax = %u\n"
+ " weighted = %d\n"
"}",
par->tolerance,
par->nrelax,
par->erelax,
par->minlevel,
- par->nitermax);
+ par->nitermax,
+ par->weighted);
}
void gfs_multilevel_params_init (GfsMultilevelParams * par)
@@ -61,6 +63,7 @@ void gfs_multilevel_params_init (GfsMultilevelParams * par)
par->nitermax = 100;
par->dimension = FTT_DIMENSION;
+ par->weighted = FALSE;
}
void gfs_multilevel_params_read (GfsMultilevelParams * par, GtsFile * fp)
@@ -71,6 +74,7 @@ void gfs_multilevel_params_read (GfsMultilevelParams * par, GtsFile * fp)
{GTS_UINT, "erelax", TRUE},
{GTS_UINT, "minlevel", TRUE},
{GTS_UINT, "nitermax", TRUE},
+ {GTS_INT, "weighted", TRUE},
{GTS_NONE}
};
@@ -82,8 +86,8 @@ void gfs_multilevel_params_read (GfsMultilevelParams * par, GtsFile * fp)
var[2].data = &par->erelax;
var[3].data = &par->minlevel;
var[4].data = &par->nitermax;
+ var[5].data = &par->weighted;
- gfs_multilevel_params_init (par);
gts_file_assign_variables (fp, var);
if (fp->type == GTS_ERROR)
return;
diff --git a/src/poisson.h b/src/poisson.h
index 243da37..9304400 100644
--- a/src/poisson.h
+++ b/src/poisson.h
@@ -40,6 +40,7 @@ struct _GfsMultilevelParams {
guint dimension;
guint niter;
guint depth;
+ gboolean weighted;
GfsNorm residual_before, residual;
};
diff --git a/src/timestep.c b/src/timestep.c
index ffca4fa..75295eb 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -62,49 +62,27 @@ static void correct_normal_velocity (FttCellFace * face,
dp /= s->solid->s[face->d];
GFS_FACE_NORMAL_VELOCITY_LEFT (face) -= dp*(*dt);
- GFS_VARIABLE (face->cell, gv[c]->i) += dp*GFS_FACE_FRACTION_LEFT (face);
-
- switch (type) {
- case FTT_FINE_FINE:
- GFS_FACE_NORMAL_VELOCITY_RIGHT (face) -= dp*(*dt);
- GFS_VARIABLE (face->neighbor, gv[c]->i) += dp*GFS_FACE_FRACTION_RIGHT (face);
- break;
- case FTT_FINE_COARSE: {
- dp *= GFS_FACE_FRACTION_LEFT (face)/(FTT_CELLS/2);
- GFS_VARIABLE (face->neighbor, gv[c]->i) += dp;
- g_assert (GFS_FACE_FRACTION_RIGHT (face) > 0.);
- GFS_FACE_NORMAL_VELOCITY_RIGHT (face) -= dp/GFS_FACE_FRACTION_RIGHT (face)*(*dt);
- break;
- }
- default:
- g_assert_not_reached ();
- }
+ GFS_VARIABLE (face->cell, gv[c]->i) += dp;
+
+ if (type == FTT_FINE_COARSE)
+ dp *= GFS_FACE_FRACTION_LEFT (face)/(GFS_FACE_FRACTION_RIGHT (face)*FTT_CELLS/2);
+ GFS_FACE_NORMAL_VELOCITY_RIGHT (face) -= dp*(*dt);
+ GFS_VARIABLE (face->neighbor, gv[c]->i) += dp;
}
static void scale_gradients (FttCell * cell, gpointer * data)
{
GfsVariable ** g = data[0];
guint * dimension = data[1];
+ FttCellNeighbors n;
FttComponent c;
- if (GFS_IS_MIXED (cell)) {
- GfsSolidVector * s = GFS_STATE (cell)->solid;
-
- for (c = 0; c < *dimension; c++) {
- g_assert (s->s[2*c] + s->s[2*c + 1] > 0.);
- GFS_VARIABLE (cell, g[c]->i) /= s->s[2*c] + s->s[2*c + 1];
- }
- }
- else {
- FttCellNeighbors n;
-
- ftt_cell_neighbors (cell, &n);
- for (c = 0; c < *dimension; c++) {
- FttCell * c1 = n.c[2*c], * c2 = n.c[2*c + 1];
-
- if (c1 && c2 && !GFS_CELL_IS_GRADIENT_BOUNDARY (c1) && !GFS_CELL_IS_GRADIENT_BOUNDARY (c2))
- GFS_VARIABLE (cell, g[c]->i) /= 2.;
- }
+ ftt_cell_neighbors (cell, &n);
+ for (c = 0; c < *dimension; c++) {
+ FttCell * c1 = n.c[2*c], * c2 = n.c[2*c + 1];
+
+ if (c1 && c2 && !GFS_CELL_IS_GRADIENT_BOUNDARY (c1) && !GFS_CELL_IS_GRADIENT_BOUNDARY (c2))
+ GFS_VARIABLE (cell, g[c]->i) /= 2.;
}
}
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list