[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:03 UTC 2009
The following commit has been merged in the upstream branch:
commit 49e728dbdd8fc11b76d2472ac3662a86374d6945
Author: Stephane Popinet <popinet at users.sf.net>
Date: Tue Mar 4 13:05:43 2008 +1100
GfsSource terms for velocity components are consistent with projection
i.e. they are computed in a manner consistent with the pressure gradient
discretisation (in a manner similar to the surface tension). This is necessary
in particular to guarantee exact hydrostatic balance with non-linear pressure
distributions (e.g. the 'hydrostatic/quadratic' test case).
darcs-hash:20080304020543-d4795-a9d9737e4dcdfeb74c3caf421b36c8a13533cbb4.gz
diff --git a/src/ocean.c b/src/ocean.c
index 312c684..009122e 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -931,9 +931,12 @@ static void gfs_source_hydrostatic_class_init (GfsSourceGenericClass * klass)
GFS_EVENT_CLASS (klass)->event = gfs_source_hydrostatic_event;
GFS_EVENT_CLASS (klass)->event_half = gfs_source_hydrostatic_event_half;
+}
- klass->mac_value = gfs_source_hydrostatic_mac_value;
- klass->centered_value = gfs_source_hydrostatic_centered_value;
+static void gfs_source_hydrostatic_init (GfsSourceGeneric * s)
+{
+ s->mac_value = gfs_source_hydrostatic_mac_value;
+ s->centered_value = gfs_source_hydrostatic_centered_value;
}
GfsSourceGenericClass * gfs_source_hydrostatic_class (void)
@@ -946,7 +949,7 @@ GfsSourceGenericClass * gfs_source_hydrostatic_class (void)
sizeof (GfsSourceHydrostatic),
sizeof (GfsSourceGenericClass),
(GtsObjectClassInitFunc) gfs_source_hydrostatic_class_init,
- (GtsObjectInitFunc) NULL,
+ (GtsObjectInitFunc) gfs_source_hydrostatic_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
@@ -1042,7 +1045,11 @@ static void gfs_source_friction_class_init (GfsSourceGenericClass * klass)
GTS_OBJECT_CLASS (klass)->read = gfs_source_friction_read;
GTS_OBJECT_CLASS (klass)->write = gfs_source_friction_write;
GFS_EVENT_CLASS (klass)->event = gfs_source_friction_event;
- klass->mac_value = klass->centered_value = gfs_source_friction_saved_value;
+}
+
+static void gfs_source_friction_init (GfsSourceGeneric * s)
+{
+ s->mac_value = s->centered_value = gfs_source_friction_saved_value;
}
GfsSourceGenericClass * gfs_source_friction_class (void)
@@ -1055,7 +1062,7 @@ GfsSourceGenericClass * gfs_source_friction_class (void)
sizeof (GfsSourceFriction),
sizeof (GfsSourceGenericClass),
(GtsObjectClassInitFunc) gfs_source_friction_class_init,
- (GtsObjectInitFunc) NULL,
+ (GtsObjectInitFunc) gfs_source_friction_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
diff --git a/src/source.c b/src/source.c
index 27bb395..51dcbfd 100644
--- a/src/source.c
+++ b/src/source.c
@@ -44,11 +44,10 @@ gdouble gfs_variable_mac_source (GfsVariable * v, FttCell * cell)
sum = 0.;
i = GTS_SLIST_CONTAINER (v->sources)->items;
while (i) {
- GtsObject * o = i->data;
+ GfsSourceGeneric * s = i->data;
- if (GFS_SOURCE_GENERIC_CLASS (o->klass)->mac_value)
- sum += (* GFS_SOURCE_GENERIC_CLASS (o->klass)->mac_value)
- (GFS_SOURCE_GENERIC (o), cell, v);
+ if (s->mac_value)
+ sum += (* s->mac_value) (s, cell, v);
i = i->next;
}
return sum;
@@ -64,11 +63,10 @@ static void add_sources (FttCell * cell, gpointer * data)
i = GTS_SLIST_CONTAINER (v->sources)->items;
while (i) {
- GtsObject * o = i->data;
+ GfsSourceGeneric * s = i->data;
- if (GFS_SOURCE_GENERIC_CLASS (o->klass)->centered_value)
- sum += (* GFS_SOURCE_GENERIC_CLASS (o->klass)->centered_value)
- (GFS_SOURCE_GENERIC (o), cell, v);
+ if (s->centered_value)
+ sum += (* s->centered_value) (s, cell, v);
i = i->next;
}
GFS_VARIABLE (cell, sv->i) += (*dt)*sum;
@@ -293,6 +291,13 @@ static void source_destroy (GtsObject * o)
(* GTS_OBJECT_CLASS (gfs_source_class ())->parent_class->destroy) (o);
}
+static gdouble source_face_value (GfsSourceGeneric * s,
+ FttCellFace * face,
+ GfsVariable * v)
+{
+ return gfs_function_face_value (GFS_SOURCE (s)->intensity, face);
+}
+
static void source_read (GtsObject ** o, GtsFile * fp)
{
if (GTS_OBJECT_CLASS (gfs_source_class ())->parent_class->read)
@@ -303,6 +308,14 @@ static void source_read (GtsObject ** o, GtsFile * fp)
GFS_SOURCE (*o)->intensity = gfs_function_new (gfs_function_class (), 0.);
gfs_function_read (GFS_SOURCE (*o)->intensity, gfs_object_simulation (*o), fp);
+ if (fp->type != GTS_ERROR) {
+ GfsSourceGeneric * s = GFS_SOURCE_GENERIC (*o);
+ gchar * name = GFS_SOURCE_SCALAR (s)->v->name;
+ if (!strcmp (name, "U") || !strcmp (name, "V") || !strcmp (name, "W")) {
+ s->mac_value = s->centered_value = NULL;
+ s->face_value = source_face_value;
+ }
+ }
}
static void source_write (GtsObject * o, FILE * fp)
@@ -313,6 +326,13 @@ static void source_write (GtsObject * o, FILE * fp)
gfs_function_write (GFS_SOURCE (o)->intensity, fp);
}
+static void source_class_init (GfsSourceGenericClass * klass)
+{
+ GTS_OBJECT_CLASS (klass)->destroy = source_destroy;
+ GTS_OBJECT_CLASS (klass)->read = source_read;
+ GTS_OBJECT_CLASS (klass)->write = source_write;
+}
+
static gdouble source_value (GfsSourceGeneric * s,
FttCell * cell,
GfsVariable * v)
@@ -320,13 +340,9 @@ static gdouble source_value (GfsSourceGeneric * s,
return gfs_function_value (GFS_SOURCE (s)->intensity, cell);
}
-static void source_class_init (GfsSourceGenericClass * klass)
+static void source_init (GfsSourceGeneric * s)
{
- GTS_OBJECT_CLASS (klass)->destroy = source_destroy;
- GTS_OBJECT_CLASS (klass)->read = source_read;
- GTS_OBJECT_CLASS (klass)->write = source_write;
-
- klass->mac_value = klass->centered_value = source_value;
+ s->mac_value = s->centered_value = source_value;
}
GfsSourceGenericClass * gfs_source_class (void)
@@ -339,7 +355,7 @@ GfsSourceGenericClass * gfs_source_class (void)
sizeof (GfsSource),
sizeof (GfsSourceGenericClass),
(GtsObjectClassInitFunc) source_class_init,
- (GtsObjectInitFunc) NULL,
+ (GtsObjectInitFunc) source_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
@@ -366,6 +382,11 @@ static gboolean source_control_event (GfsEvent * event, GfsSimulation * sim)
return FALSE;
}
+static void source_control_class_init (GfsSourceGenericClass * klass)
+{
+ GFS_EVENT_CLASS (klass)->event = source_control_event;
+}
+
static gdouble source_control_value (GfsSourceGeneric * s,
FttCell * cell,
GfsVariable * v)
@@ -373,10 +394,9 @@ static gdouble source_control_value (GfsSourceGeneric * s,
return GFS_SOURCE_CONTROL (s)->s;
}
-static void source_control_class_init (GfsSourceGenericClass * klass)
+static void source_control_init (GfsSourceGeneric * s)
{
- GFS_EVENT_CLASS (klass)->event = source_control_event;
- klass->mac_value = klass->centered_value = source_control_value;
+ s->mac_value = s->centered_value = source_control_value;
}
GfsSourceGenericClass * gfs_source_control_class (void)
@@ -389,7 +409,7 @@ GfsSourceGenericClass * gfs_source_control_class (void)
sizeof (GfsSourceControl),
sizeof (GfsSourceGenericClass),
(GtsObjectClassInitFunc) source_control_class_init,
- (GtsObjectInitFunc) NULL,
+ (GtsObjectInitFunc) source_control_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
@@ -639,13 +659,12 @@ static void source_diffusion_class_init (GfsSourceGenericClass * klass)
GTS_OBJECT_CLASS (klass)->write = source_diffusion_write;
GFS_EVENT_CLASS (klass)->event = source_diffusion_event;
-
- klass->mac_value = source_diffusion_value;
}
static void source_diffusion_init (GfsSourceDiffusion * d)
{
d->D = GFS_DIFFUSION (gts_object_new (GTS_OBJECT_CLASS (gfs_diffusion_class ())));
+ GFS_SOURCE_GENERIC (d)->mac_value = source_diffusion_value;
}
GfsSourceGenericClass * gfs_source_diffusion_class (void)
@@ -783,10 +802,14 @@ static void gfs_source_diffusion_explicit_class_init (GfsSourceGenericClass * kl
GFS_EVENT_CLASS (klass)->event = gfs_source_diffusion_explicit_event;
GTS_OBJECT_CLASS (klass)->read = gfs_source_diffusion_explicit_read;
GTS_OBJECT_CLASS (klass)->destroy = gfs_source_diffusion_explicit_destroy;
- klass->mac_value = klass->centered_value = source_diffusion_explicit_value;
klass->stability = source_diffusion_stability;
}
+static void gfs_source_diffusion_explicit_init (GfsSourceGeneric * s)
+{
+ s->mac_value = s->centered_value = source_diffusion_explicit_value;
+}
+
GfsSourceGenericClass * gfs_source_diffusion_explicit_class (void)
{
static GfsSourceGenericClass * klass = NULL;
@@ -797,7 +820,7 @@ GfsSourceGenericClass * gfs_source_diffusion_explicit_class (void)
sizeof (GfsSourceDiffusionExplicit),
sizeof (GfsSourceGenericClass),
(GtsObjectClassInitFunc) gfs_source_diffusion_explicit_class_init,
- (GtsObjectInitFunc) NULL,
+ (GtsObjectInitFunc) gfs_source_diffusion_explicit_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
@@ -881,9 +904,12 @@ static void source_viscosity_class_init (GfsSourceGenericClass * klass)
{
GTS_OBJECT_CLASS (klass)->read = source_viscosity_read;
GTS_OBJECT_CLASS (klass)->write = source_viscosity_write;
+}
- klass->mac_value = source_viscosity_value;
- klass->centered_value = source_viscosity_non_diffusion_value;
+static void source_viscosity_init (GfsSourceGeneric * s)
+{
+ s->mac_value = source_viscosity_value;
+ s->centered_value = source_viscosity_non_diffusion_value;
}
GfsSourceGenericClass * gfs_source_viscosity_class (void)
@@ -896,7 +922,7 @@ GfsSourceGenericClass * gfs_source_viscosity_class (void)
sizeof (GfsSourceViscosity),
sizeof (GfsSourceGenericClass),
(GtsObjectClassInitFunc) source_viscosity_class_init,
- (GtsObjectInitFunc) NULL,
+ (GtsObjectInitFunc) source_viscosity_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
@@ -924,10 +950,14 @@ static gdouble source_viscosity_stability (GfsSourceGeneric * s,
static void source_viscosity_explicit_class_init (GfsSourceGenericClass * klass)
{
- klass->mac_value = klass->centered_value = source_viscosity_value;
klass->stability = source_viscosity_stability;
}
+static void source_viscosity_explicit_init (GfsSourceGeneric * s)
+{
+ s->mac_value = s->centered_value = source_viscosity_value;
+}
+
GfsSourceGenericClass * gfs_source_viscosity_explicit_class (void)
{
static GfsSourceGenericClass * klass = NULL;
@@ -938,7 +968,7 @@ GfsSourceGenericClass * gfs_source_viscosity_explicit_class (void)
sizeof (GfsSourceViscosity),
sizeof (GfsSourceGenericClass),
(GtsObjectClassInitFunc) source_viscosity_explicit_class_init,
- (GtsObjectInitFunc) NULL,
+ (GtsObjectInitFunc) source_viscosity_explicit_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
@@ -1071,8 +1101,12 @@ static void gfs_source_coriolis_class_init (GfsSourceGenericClass * klass)
GTS_OBJECT_CLASS (klass)->write = gfs_source_coriolis_write;
GFS_EVENT_CLASS (klass)->event = gfs_source_coriolis_event;
- klass->mac_value = gfs_source_coriolis_mac_value;
- klass->centered_value = gfs_source_coriolis_centered_value;
+}
+
+static void gfs_source_coriolis_init (GfsSourceGeneric * s)
+{
+ s->mac_value = gfs_source_coriolis_mac_value;
+ s->centered_value = gfs_source_coriolis_centered_value;
}
GfsSourceGenericClass * gfs_source_coriolis_class (void)
@@ -1085,7 +1119,7 @@ GfsSourceGenericClass * gfs_source_coriolis_class (void)
sizeof (GfsSourceCoriolis),
sizeof (GfsSourceGenericClass),
(GtsObjectClassInitFunc) gfs_source_coriolis_class_init,
- (GtsObjectInitFunc) NULL,
+ (GtsObjectInitFunc) gfs_source_coriolis_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
diff --git a/src/source.h b/src/source.h
index b1c49db..903035a 100644
--- a/src/source.h
+++ b/src/source.h
@@ -42,6 +42,9 @@ struct _GfsSourceGeneric {
/*< public >*/
GfsVariable * v;
+ gdouble (* mac_value) (GfsSourceGeneric *, FttCell *, GfsVariable *);
+ gdouble (* centered_value) (GfsSourceGeneric *, FttCell *, GfsVariable *);
+ gdouble (* face_value) (GfsSourceGeneric *, FttCellFace *, GfsVariable *);
};
typedef struct _GfsSourceGenericClass GfsSourceGenericClass;
@@ -51,8 +54,6 @@ struct _GfsSourceGenericClass {
GfsEventClass parent_class;
/*< public >*/
- gdouble (* mac_value) (GfsSourceGeneric *, FttCell *, GfsVariable *);
- gdouble (* centered_value) (GfsSourceGeneric *, FttCell *, GfsVariable *);
gdouble (* stability) (GfsSourceGeneric *, GfsSimulation *);
};
diff --git a/src/tension.c b/src/tension.c
index 697f7fc..5936c1d 100644
--- a/src/tension.c
+++ b/src/tension.c
@@ -221,7 +221,11 @@ static void gfs_source_tension_css_class_init (GfsSourceGenericClass * klass)
{
GTS_OBJECT_CLASS (klass)->read = gfs_source_tension_css_read;
GFS_EVENT_CLASS (klass)->event = gfs_source_tension_css_event;
- klass->centered_value = gfs_source_tension_css_value;
+}
+
+static void gfs_source_tension_css_init (GfsSourceGeneric * s)
+{
+ s->centered_value = gfs_source_tension_css_value;
}
GfsSourceGenericClass * gfs_source_tension_css_class (void)
@@ -234,7 +238,7 @@ GfsSourceGenericClass * gfs_source_tension_css_class (void)
sizeof (GfsSourceTensionCSS),
sizeof (GfsSourceGenericClass),
(GtsObjectClassInitFunc) gfs_source_tension_css_class_init,
- (GtsObjectInitFunc) NULL,
+ (GtsObjectInitFunc) gfs_source_tension_css_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
diff --git a/src/timestep.c b/src/timestep.c
index b4be6a0..8be40c2 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -174,14 +174,62 @@ static void scale_divergence (FttCell * cell, gpointer * data)
GFS_VARIABLE (cell, div->i) /= *dt;
}
-static void surface_tension (GfsDomain * domain,
- GfsVariable * u,
- gdouble dt,
- GfsFunction * alpha,
- GfsVariable ** g)
-{
- if (u->sources) {
- GSList * i = GTS_SLIST_CONTAINER (u->sources)->items;
+typedef struct {
+ GfsSourceGeneric * s;
+ GfsVariable * v, ** g;
+ gdouble dt;
+} FaceSource;
+
+static void add_face_source (FttCellFace * face,
+ FaceSource * f)
+{
+ gdouble dp;
+ FttComponent c;
+
+ if (GFS_FACE_FRACTION (face) == 0.)
+ return;
+
+ c = face->d/2;
+ dp = (* f->s->face_value) (f->s, face, f->v);
+ GFS_FACE_NORMAL_VELOCITY_LEFT (face) -= dp*f->dt;
+ if (f->g)
+ GFS_VARIABLE (face->cell, f->g[c]->i) += dp*GFS_FACE_FRACTION_LEFT (face);
+
+ if (ftt_face_type (face) == FTT_FINE_COARSE)
+ dp *= GFS_FACE_FRACTION_LEFT (face)/(GFS_FACE_FRACTION_RIGHT (face)*FTT_CELLS/2);
+ GFS_FACE_NORMAL_VELOCITY_RIGHT (face) -= dp*f->dt;
+ if (f->g)
+ GFS_VARIABLE (face->neighbor, f->g[c]->i) += dp*GFS_FACE_FRACTION_RIGHT (face);
+}
+
+static void velocity_face_sources (GfsDomain * domain,
+ GfsVariable ** u,
+ gdouble dt,
+ GfsFunction * alpha,
+ GfsVariable ** g)
+{
+ FttComponent c;
+ for (c = 0; c < FTT_DIMENSION; c++)
+ if (u[c]->sources) {
+ GSList * i = GTS_SLIST_CONTAINER (u[c]->sources)->items;
+
+ while (i) {
+ GfsSourceGeneric * s = i->data;
+ if (s->face_value) {
+ FaceSource f;
+ f.s = s;
+ f.v = u[c];
+ f.g = g;
+ f.dt = dt;
+ gfs_domain_face_traverse (domain, c,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttFaceTraverseFunc) add_face_source, &f);
+ }
+ i = i->next;
+ }
+ }
+ if (u[0]->sources) {
+ GSList * i = GTS_SLIST_CONTAINER (u[0]->sources)->items;
while (i) {
if (GFS_IS_SOURCE_TENSION (i->data)) {
@@ -214,9 +262,9 @@ void gfs_update_gradients (GfsDomain * domain,
g_return_if_fail (p != NULL);
g_return_if_fail (g != NULL);
- /* Add surface tension */
+ /* Add face sources */
reset_gradients (domain, FTT_DIMENSION, g);
- surface_tension (domain, gfs_domain_velocity (domain)[0], 0., alpha, g);
+ velocity_face_sources (domain, gfs_domain_velocity (domain), 0., alpha, g);
/* Initialize face coefficients */
gfs_poisson_coefficients (domain, alpha);
/* Add pressure gradient */
@@ -232,9 +280,9 @@ static void mac_projection (GfsDomain * domain,
GfsVariable * res,
GfsVariable ** g)
{
- /* Add surface tension */
+ /* Add face sources */
reset_gradients (domain, FTT_DIMENSION, g);
- surface_tension (domain, gfs_domain_velocity (domain)[0], apar->dt, alpha, g);
+ velocity_face_sources (domain, gfs_domain_velocity (domain), apar->dt, alpha, g);
GfsVariable * div = gfs_temporary_variable (domain);
GfsVariable * dia = gfs_temporary_variable (domain);
diff --git a/test/hydrostatic/hydrostatic.gfs b/test/hydrostatic/hydrostatic.gfs
index 8090ea9..0391f5c 100644
--- a/test/hydrostatic/hydrostatic.gfs
+++ b/test/hydrostatic/hydrostatic.gfs
@@ -18,12 +18,8 @@
ApproxProjectionParams { tolerance = 1e-12 }
ProjectionParams { tolerance = 1e-12 }
- # The pressure correction needs to be included in the
- # Crank-Nicholson scheme in order to balance the body force
- AdvectionParams { gc = 1 }
-
OutputScalarNorm { istep = 1 } v { v = V }
- EventScript { start = end } {
+ EventScript { start = end } {
if awk '{if ($9 > 1e-12) exit (1);}' < v ; then
return 0;
else
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list