[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:53:37 UTC 2009
The following commit has been merged in the upstream branch:
commit 81195cff62702614c4d1ffead9b950977d3de6e2
Author: Stephane Popinet <popinet at users.sf.net>
Date: Mon Feb 6 00:08:32 2006 +1100
CSF uses volume-fraction-weighted tension terms
darcs-hash:20060205130832-d4795-433cf2f7616fd446e9fd4151d8e5be33ac2a223c.gz
diff --git a/configure.in b/configure.in
index 4a8ee02..41ada63 100644
--- a/configure.in
+++ b/configure.in
@@ -12,10 +12,10 @@ dnl are available for $ac_help expansion (don't we all *love* autoconf?)
# set GFS_BINARY_AGE and GFS_INTERFACE_AGE to 0.
#
GFS_MAJOR_VERSION=0
-GFS_MINOR_VERSION=9
-GFS_MICRO_VERSION=3
-GFS_INTERFACE_AGE=0
-GFS_BINARY_AGE=0
+GFS_MINOR_VERSION=8
+GFS_MICRO_VERSION=1
+GFS_INTERFACE_AGE=1
+GFS_BINARY_AGE=1
GFS_VERSION=$GFS_MAJOR_VERSION.$GFS_MINOR_VERSION.$GFS_MICRO_VERSION
GFS_COMPILATION_FLAGS=$CFLAGS
dnl
diff --git a/src/ocean.c b/src/ocean.c
index 9439d7f..9c6a7d5 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -112,7 +112,7 @@ static void gfs_correct_normal_velocities_weighted (GfsDomain * domain,
gboolean weighted)
{
if (weighted)
- gfs_correct_normal_velocities (domain, dimension, p, g, dt, NULL);
+ gfs_correct_normal_velocities (domain, dimension, p, g, dt);
else {
gpointer data[3];
FttComponent c;
@@ -658,7 +658,7 @@ static void ocean_run (GfsSimulation * sim)
gfs_domain_timer_start (domain, "correct_normal_velocities");
gfs_poisson_coefficients (domain, NULL);
- gfs_correct_normal_velocities_weighted (domain, 2, p, g, sim->advection_params.dt/2., FALSE);
+ gfs_correct_normal_velocities (domain, 2, p, g, sim->advection_params.dt/2., NULL);
gfs_domain_cell_traverse_boundary (domain, FTT_BACK,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) compute_w,
diff --git a/src/poisson.c b/src/poisson.c
index 8442152..f9d71fd 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -21,6 +21,8 @@
#include "poisson.h"
#include "solid.h"
#include "source.h"
+#include "tension.h"
+#include "vof.h"
/**
* gfs_multilevel_params_write:
@@ -445,6 +447,94 @@ void gfs_poisson_coefficients (GfsDomain * domain,
(FttCellTraverseFunc) face_coeff_from_below, NULL);
}
+static void tension_coeff (FttCellFace * face, gpointer * data)
+{
+ gdouble * lambda2 = data[0];
+ GfsStateVector * s = GFS_STATE (face->cell);
+ gdouble v = lambda2[face->d/2];
+ GfsSourceTension * t = data[1];
+ gdouble c1 = GFS_VARIABLE (face->cell, t->c->i);
+ gdouble c2 = GFS_VARIABLE (face->neighbor, t->c->i);
+ gdouble w1 = GFS_IS_FULL (c1) ? 0. : c1*(1. - c1);
+ gdouble w2 = GFS_IS_FULL (c2) ? 0. : c2*(1. - c2);
+
+ if (w1 + w2 > 0.) {
+ GfsVariable * alpha = data[2];
+
+ v *= (w1*GFS_VARIABLE (face->cell, t->k->i) +
+ w2*GFS_VARIABLE (face->neighbor, t->k->i))/(w1 + w2);
+ if (alpha)
+ v *= gfs_face_interpolated_value (face, alpha->i);
+ }
+ else
+ v = 0.;
+
+ if (GFS_IS_MIXED (face->cell))
+ v *= s->solid->s[face->d];
+ s->f[face->d].v = v;
+
+ switch (ftt_face_type (face)) {
+ case FTT_FINE_FINE:
+ GFS_STATE (face->neighbor)->f[FTT_OPPOSITE_DIRECTION (face->d)].v = v;
+ break;
+ case FTT_FINE_COARSE:
+ GFS_STATE (face->neighbor)->f[FTT_OPPOSITE_DIRECTION (face->d)].v = G_MAXDOUBLE;
+ break;
+ default:
+ g_assert_not_reached ();
+ }
+}
+
+/**
+ * gfs_source_tension_coefficients:
+ * @s: a #GfsSourceTension.
+ * @domain: a #GfsDomain.
+ * @alpha: the inverse of density or %NULL.
+ *
+ * Initializes the face coefficients with the surface tension term
+ * (interface curvature times surface tension coefficient).
+ *
+ * If @alpha is %NULL, it is taken to be unity.
+ */
+void gfs_source_tension_coefficients (GfsSourceTension * s,
+ GfsDomain * domain,
+ GfsFunction * alpha)
+{
+ gdouble lambda2[FTT_DIMENSION];
+ gpointer data[3];
+ FttComponent i;
+
+ g_return_if_fail (s != NULL);
+ g_return_if_fail (domain != NULL);
+
+ for (i = 0; i < FTT_DIMENSION; i++) {
+ gdouble lambda = (&domain->lambda.x)[i];
+
+ lambda2[i] = lambda*lambda;
+ }
+ if (alpha) {
+ data[0] = alpha;
+ data[1] = gfs_temporary_variable (domain);
+ gfs_domain_cell_traverse (domain,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) reset_alpha_coeff, data);
+ data[2] = data[1];
+ }
+ else {
+ gfs_domain_cell_traverse (domain,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) reset_coeff, NULL);
+ data[2] = NULL;
+ }
+ data[0] = lambda2;
+ data[1] = s;
+ gfs_domain_face_traverse (domain, FTT_XYZ,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttFaceTraverseFunc) tension_coeff, data);
+ if (alpha)
+ gts_object_destroy (data[2]);
+}
+
static void correct (FttCell * cell, gpointer * data)
{
GfsVariable * u = data[0];
diff --git a/src/tension.c b/src/tension.c
index 684d870..777f91a 100644
--- a/src/tension.c
+++ b/src/tension.c
@@ -227,7 +227,7 @@ static gdouble gfs_source_tension_stability (GfsSourceGeneric * s,
GfsSimulation * sim)
{
GfsSourceTension * t = GFS_SOURCE_TENSION (s);
- gdouble sigma = GFS_VARIABLE_CURVATURE (t->k)->sigma, h;
+ gdouble h, sigma = 1.;
StabilityParams p = { G_MAXDOUBLE, -G_MAXDOUBLE, 0 };
p.alpha = sim->physical_params.alpha;
diff --git a/src/tension.h b/src/tension.h
index fd2bf6d..b251734 100644
--- a/src/tension.h
+++ b/src/tension.h
@@ -66,7 +66,10 @@ struct _GfsSourceTension {
#define GFS_IS_SOURCE_TENSION(obj) (gts_object_is_from_class (obj,\
gfs_source_tension_class ()))
-GfsSourceGenericClass * gfs_source_tension_class (void);
+GfsSourceGenericClass * gfs_source_tension_class (void);
+void gfs_source_tension_coefficients (GfsSourceTension * s,
+ GfsDomain * domain,
+ GfsFunction * alpha);
#ifdef __cplusplus
}
diff --git a/src/timestep.c b/src/timestep.c
index b86b003..97fd465 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -41,7 +41,7 @@ static void correct_normal_velocity (FttCellFace * face,
GfsGradient g;
gdouble dp;
GfsStateVector * s;
- GfsVariable * p = data[0], ** gv = data[1], * w = data[3];
+ GfsVariable * p = data[0], ** gv = data[1];
gdouble * dt = data[2];
FttComponent c;
@@ -53,8 +53,6 @@ static void correct_normal_velocity (FttCellFace * face,
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 (w)
- dp *= gfs_face_interpolated_value (face, w->i);
if (!FTT_FACE_DIRECT (face))
dp = - dp;
@@ -107,7 +105,6 @@ static void scale_gradients (FttCell * cell, gpointer * data)
* @p: the pressure field.
* @g: where to store the pressure gradient or %NULL.
* @dt: the timestep.
- * @w: an optional weight to apply to the correction.
*
* Corrects the normal velocity field of @domain using @p and and @dt.
*
@@ -118,10 +115,9 @@ void gfs_correct_normal_velocities (GfsDomain * domain,
guint dimension,
GfsVariable * p,
GfsVariable ** g,
- gdouble dt,
- GfsVariable * w)
+ gdouble dt)
{
- gpointer data[4];
+ gpointer data[3];
FttComponent c;
g_return_if_fail (domain != NULL);
@@ -140,7 +136,6 @@ void gfs_correct_normal_velocities (GfsDomain * domain,
data[0] = p;
data[1] = g;
data[2] = &dt;
- data[3] = w;
gfs_domain_face_traverse (domain, dimension == 2 ? FTT_XY : FTT_XYZ,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttFaceTraverseFunc) correct_normal_velocity, data);
@@ -215,16 +210,17 @@ void gfs_mac_projection (GfsDomain * domain,
dt = apar->dt;
apar->dt /= 2.;
+ /* Add surface tension */
+ if ((s = gfs_source_find (apar->v, gfs_source_tension_class ()))) {
+ gfs_source_tension_coefficients (GFS_SOURCE_TENSION (s), domain, alpha);
+ gfs_correct_normal_velocities (domain, FTT_DIMENSION,
+ GFS_SOURCE_TENSION (s)->c,
+ NULL, apar->dt);
+ }
/* Initialize face coefficients */
gfs_poisson_coefficients (domain, alpha);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) gfs_cell_reset, dia);
- /* Add surface tension */
- if ((s = gfs_source_find (apar->v, gfs_source_tension_class ())))
- gfs_correct_normal_velocities (domain, FTT_DIMENSION,
- GFS_SOURCE_TENSION (s)->c,
- NULL, apar->dt,
- GFS_SOURCE_TENSION (s)->k);
/* compute MAC divergence */
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) gfs_normal_divergence, div);
@@ -263,7 +259,7 @@ void gfs_mac_projection (GfsDomain * domain,
gts_object_destroy (GTS_OBJECT (dia));
gts_object_destroy (GTS_OBJECT (res));
- gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, apar->dt, NULL);
+ gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, apar->dt);
#if 0
{
@@ -374,9 +370,6 @@ void gfs_approximate_projection (GfsDomain * domain,
gfs_domain_timer_start (domain, "approximate_projection");
- /* Initialize face coefficients */
- gfs_poisson_coefficients (domain, alpha);
-
/* compute MAC velocities from centered velocities */
gfs_domain_face_traverse (domain, FTT_XYZ,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
@@ -387,16 +380,18 @@ void gfs_approximate_projection (GfsDomain * domain,
gfs_domain_velocity (domain));
/* Add surface tension */
if ((s = gfs_source_find (gfs_domain_velocity (domain)[0], gfs_source_tension_class ()))) {
+ gfs_source_tension_coefficients (GFS_SOURCE_TENSION (s), domain, alpha);
gfs_correct_normal_velocities (domain, FTT_DIMENSION,
GFS_SOURCE_TENSION (s)->c,
- g, apar->dt,
- GFS_SOURCE_TENSION (s)->k);
+ g, apar->dt);
gfs_correct_centered_velocities (domain, FTT_DIMENSION, g, apar->dt);
}
div = gfs_temporary_variable (domain);
dia = gfs_temporary_variable (domain);
res1 = res ? res : gfs_temporary_variable (domain);
+ /* Initialize face coefficients */
+ gfs_poisson_coefficients (domain, alpha);
/* Initialize diagonal coefficient */
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
@@ -449,7 +444,7 @@ void gfs_approximate_projection (GfsDomain * domain,
if (!res)
gts_object_destroy (GTS_OBJECT (res1));
- gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, apar->dt, NULL);
+ gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, apar->dt);
gfs_correct_centered_velocities (domain, FTT_DIMENSION, g, apar->dt);
gfs_domain_timer_stop (domain, "approximate_projection");
diff --git a/src/timestep.h b/src/timestep.h
index d692c57..b13e2c8 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -32,8 +32,7 @@ void gfs_correct_normal_velocities (GfsDomain * domain,
guint dimension,
GfsVariable * p,
GfsVariable ** g,
- gdouble dt,
- GfsVariable * w);
+ gdouble dt);
void gfs_mac_projection (GfsDomain * domain,
GfsMultilevelParams * par,
GfsAdvectionParams * apar,
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list