[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:53:43 UTC 2009
The following commit has been merged in the upstream branch:
commit 6e66e72551a4292385c9594bc80c6f0759168bc6
Author: Stephane Popinet <popinet at users.sf.net>
Date: Thu Oct 19 14:16:12 2006 +1000
SourceTension and SourceTensionCSS are derived from a new SourceTensionGeneric object
Also SourceTensionCSS "works" again thanks to the new function
gfs_youngs_gradient (actually a renamed version of the old
gfs_youngs_normal implementation).
darcs-hash:20061019041612-d4795-88be66033a317c8d7791d246f2b718051ab92e6f.gz
diff --git a/src/poisson.c b/src/poisson.c
index 605f588..871a4d7 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -22,7 +22,6 @@
#include "solid.h"
#include "source.h"
#include "tension.h"
-#include "vof.h"
/**
* gfs_multilevel_params_write:
@@ -451,26 +450,26 @@ static void tension_coeff (FttCellFace * face, gpointer * data)
{
gdouble * lambda2 = data[0];
GfsStateVector * s = GFS_STATE (face->cell);
- GfsSourceTension * t = data[1];
+ GfsSourceTensionGeneric * t = data[1];
gdouble v = lambda2[face->d/2]*t->sigma;
- GfsVariable * alpha = data[2];
+ GfsVariable * alpha = data[2], * kappa = GFS_SOURCE_TENSION (data[1])->k;
gdouble c1 = GFS_VARIABLE (face->cell, t->c->i);
gdouble c2 = GFS_VARIABLE (face->neighbor, t->c->i);
gdouble w1 = c1*(1. - c1);
gdouble w2 = c2*(1. - c2);
if (w1 + w2 > 0.)
- v *= (w1*GFS_VARIABLE (face->cell, t->k->i) +
- w2*GFS_VARIABLE (face->neighbor, t->k->i))/(w1 + w2);
+ v *= (w1*GFS_VARIABLE (face->cell, kappa->i) +
+ w2*GFS_VARIABLE (face->neighbor, kappa->i))/(w1 + w2);
else {
- if (GFS_VARIABLE (face->cell, t->k->i) < G_MAXDOUBLE) {
- if (GFS_VARIABLE (face->neighbor, t->k->i) < G_MAXDOUBLE)
- v *= (GFS_VARIABLE (face->cell, t->k->i) + GFS_VARIABLE (face->neighbor, t->k->i))/2.;
+ if (GFS_VARIABLE (face->cell, kappa->i) < G_MAXDOUBLE) {
+ if (GFS_VARIABLE (face->neighbor, kappa->i) < G_MAXDOUBLE)
+ v *= (GFS_VARIABLE (face->cell, kappa->i) + GFS_VARIABLE (face->neighbor, kappa->i))/2.;
else
- v *= GFS_VARIABLE (face->cell, t->k->i);
+ v *= GFS_VARIABLE (face->cell, kappa->i);
}
- else if (GFS_VARIABLE (face->neighbor, t->k->i) < G_MAXDOUBLE)
- v *= GFS_VARIABLE (face->neighbor, t->k->i);
+ else if (GFS_VARIABLE (face->neighbor, kappa->i) < G_MAXDOUBLE)
+ v *= GFS_VARIABLE (face->neighbor, kappa->i);
else
v = 1e6;
}
diff --git a/src/tension.c b/src/tension.c
index dba6f9e..32a9383 100644
--- a/src/tension.c
+++ b/src/tension.c
@@ -23,15 +23,14 @@
#include "tension.h"
#include "vof.h"
-/* GfsSourceTensionCSS: Object */
+/* GfsSourceTensionGeneric: Object */
-static void gfs_source_tension_css_read (GtsObject ** o, GtsFile * fp)
+static void gfs_source_tension_generic_read (GtsObject ** o, GtsFile * fp)
{
- GfsSourceTensionCSS * s = GFS_SOURCE_TENSION_CSS (*o);
+ GfsSourceTensionGeneric * s = GFS_SOURCE_TENSION_GENERIC (*o);
GfsDomain * domain = GFS_DOMAIN (gfs_object_simulation (*o));
- FttComponent c;
- (* GTS_OBJECT_CLASS (gfs_source_tension_css_class ())->parent_class->read) (o, fp);
+ (* GTS_OBJECT_CLASS (gfs_source_tension_generic_class ())->parent_class->read) (o, fp);
if (fp->type == GTS_ERROR)
return;
@@ -45,12 +44,6 @@ static void gfs_source_tension_css_read (GtsObject ** o, GtsFile * fp)
}
gts_file_next_token (fp);
- for (c = 0; c < FTT_DIMENSION; c++) {
- static gchar * name[3] = {"_Tx", "_Ty", "_Tz"};
- if ((s->t[c] = gfs_variable_from_name (domain->variables, name[c])) == NULL)
- s->t[c] = gfs_domain_add_variable (domain, name[c], NULL);
- }
-
if (fp->type != GTS_INT && fp->type != GTS_FLOAT) {
gts_file_error (fp, "expecting a number (sigma)");
return;
@@ -60,20 +53,110 @@ static void gfs_source_tension_css_read (GtsObject ** o, GtsFile * fp)
gts_file_next_token (fp);
}
-static void gfs_source_tension_css_write (GtsObject * o, FILE * fp)
+static void gfs_source_tension_generic_write (GtsObject * o, FILE * fp)
+{
+ GfsSourceTensionGeneric * t = GFS_SOURCE_TENSION_GENERIC (o);
+ (* GTS_OBJECT_CLASS (gfs_source_tension_generic_class ())->parent_class->write) (o, fp);
+ fprintf (fp, " %s %g", t->c->name, t->sigma);
+}
+
+typedef struct {
+ gdouble amin, amax;
+ guint depth;
+ GfsFunction * alpha;
+ GfsVariable * c;
+} StabilityParams;
+
+static void min_max_alpha (FttCell * cell, StabilityParams * p)
+{
+ guint level = ftt_cell_level (cell);
+ if (level > p->depth &&
+ GFS_VARIABLE (cell, p->c->i) > 1e-3 &&
+ GFS_VARIABLE (cell, p->c->i) < 1. - 1.e-3)
+ p->depth = level;
+ if (p->alpha) {
+ gdouble a = gfs_function_value (p->alpha, cell);
+ if (a < p->amin) p->amin = a;
+ if (a > p->amax) p->amax = a;
+ }
+}
+
+static gdouble gfs_source_tension_generic_stability (GfsSourceGeneric * s,
+ GfsSimulation * sim)
+{
+ GfsSourceTensionGeneric * t = GFS_SOURCE_TENSION_GENERIC (s);
+ gdouble h;
+ StabilityParams p = { G_MAXDOUBLE, -G_MAXDOUBLE, 0 };
+
+ p.alpha = sim->physical_params.alpha;
+ p.c = t->c;
+ gfs_domain_cell_traverse (GFS_DOMAIN (sim), FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) min_max_alpha, &p);
+ h = ftt_level_size (p.depth);
+ if (p.alpha) {
+ gdouble rhom = (1./p.amin + 1./p.amax)/2.;
+ return sqrt (rhom*h*h*h/(2.*M_PI*t->sigma));
+ }
+ else
+ return sqrt (h*h*h/(2.*M_PI*t->sigma));
+}
+
+static void gfs_source_tension_generic_class_init (GfsSourceGenericClass * klass)
+{
+ GTS_OBJECT_CLASS (klass)->read = gfs_source_tension_generic_read;
+ GTS_OBJECT_CLASS (klass)->write = gfs_source_tension_generic_write;
+ klass->stability = gfs_source_tension_generic_stability;
+}
+
+GfsSourceGenericClass * gfs_source_tension_generic_class (void)
+{
+ static GfsSourceGenericClass * klass = NULL;
+
+ if (klass == NULL) {
+ GtsObjectClassInfo gfs_source_tension_generic_info = {
+ "GfsSourceTensionGeneric",
+ sizeof (GfsSourceTensionGeneric),
+ sizeof (GfsSourceGenericClass),
+ (GtsObjectClassInitFunc) gfs_source_tension_generic_class_init,
+ (GtsObjectInitFunc) NULL,
+ (GtsArgSetFunc) NULL,
+ (GtsArgGetFunc) NULL
+ };
+ klass =
+ gts_object_class_new (GTS_OBJECT_CLASS (gfs_source_velocity_class ()),
+ &gfs_source_tension_generic_info);
+ }
+
+ return klass;
+}
+
+/* GfsSourceTensionCSS: Object */
+
+static void gfs_source_tension_css_read (GtsObject ** o, GtsFile * fp)
{
- (* GTS_OBJECT_CLASS (gfs_source_tension_css_class ())->parent_class->write) (o, fp);
- fprintf (fp, " %s %g", GFS_SOURCE_TENSION_CSS (o)->c->name, GFS_SOURCE_TENSION_CSS (o)->sigma);
+ GfsSourceTensionCSS * s = GFS_SOURCE_TENSION_CSS (*o);
+ GfsDomain * domain = GFS_DOMAIN (gfs_object_simulation (*o));
+ FttComponent c;
+
+ (* GTS_OBJECT_CLASS (gfs_source_tension_css_class ())->parent_class->read) (o, fp);
+ if (fp->type == GTS_ERROR)
+ return;
+
+ for (c = 0; c < FTT_DIMENSION; c++) {
+ static gchar * name[3] = {"_Tx", "_Ty", "_Tz"};
+ if ((s->t[c] = gfs_variable_from_name (domain->variables, name[c])) == NULL)
+ s->t[c] = gfs_domain_add_variable (domain, name[c], NULL);
+ }
}
static void foreach_cell_normal (FttCell * cell, GfsSourceTensionCSS * s)
{
FttVector n;
gdouble nn = 0.;
- gdouble sigh = s->sigma/ftt_cell_size (cell);
+ gdouble sigh = GFS_SOURCE_TENSION_GENERIC (s)->sigma/ftt_cell_size (cell);
FttComponent c;
- gfs_youngs_normal (cell, s->c, &n);
+ gfs_youngs_normal (cell, GFS_SOURCE_TENSION_GENERIC (s)->c, &n);
for (c = 0; c < FTT_DIMENSION; c++)
nn += (&n.x)[c]*(&n.x)[c];
nn = sqrt (nn + 1e-50);
@@ -90,9 +173,9 @@ static void foreach_cell_tension_css (FttCell * cell, GfsSourceTensionCSS * s)
gdouble alpha = sim->physical_params.alpha ?
gfs_function_value (sim->physical_params.alpha, cell) : 1.;
- gfs_youngs_normal (cell, s->g[0], &nx);
- gfs_youngs_normal (cell, s->g[1], &ny);
- gfs_youngs_normal (cell, s->g[2], &nxy);
+ gfs_youngs_gradient (cell, s->g[0], &nx);
+ gfs_youngs_gradient (cell, s->g[1], &ny);
+ gfs_youngs_gradient (cell, s->g[2], &nxy);
GFS_VARIABLE (cell, s->t[0]->i) = alpha*(ny.x - nxy.y)/h;
GFS_VARIABLE (cell, s->t[1]->i) = alpha*(nx.y - nxy.x)/h;
@@ -132,7 +215,6 @@ static gdouble gfs_source_tension_css_value (GfsSourceGeneric * s,
static void gfs_source_tension_css_class_init (GfsSourceGenericClass * klass)
{
GTS_OBJECT_CLASS (klass)->read = gfs_source_tension_css_read;
- GTS_OBJECT_CLASS (klass)->write = gfs_source_tension_css_write;
GFS_EVENT_CLASS (klass)->event_half = gfs_source_tension_css_event;
klass->centered_value = gfs_source_tension_css_value;
}
@@ -152,7 +234,7 @@ GfsSourceGenericClass * gfs_source_tension_css_class (void)
(GtsArgGetFunc) NULL
};
klass =
- gts_object_class_new (GTS_OBJECT_CLASS (gfs_source_velocity_class ()),
+ gts_object_class_new (GTS_OBJECT_CLASS (gfs_source_tension_generic_class ()),
&gfs_source_tension_css_info);
}
@@ -171,16 +253,6 @@ static void gfs_source_tension_read (GtsObject ** o, GtsFile * fp)
return;
if (fp->type != GTS_STRING) {
- gts_file_error (fp, "expecting a variable (C)");
- return;
- }
- if ((s->c = gfs_variable_from_name (domain->variables, fp->token->str)) == NULL) {
- gts_file_error (fp, "unknown variable `%s'", fp->token->str);
- return;
- }
- gts_file_next_token (fp);
-
- if (fp->type != GTS_STRING) {
gts_file_error (fp, "expecting a variable (Kappa)");
return;
}
@@ -193,68 +265,18 @@ static void gfs_source_tension_read (GtsObject ** o, GtsFile * fp)
return;
}
gts_file_next_token (fp);
-
- if (fp->type != GTS_INT && fp->type != GTS_FLOAT) {
- gts_file_error (fp, "expecting a number (sigma)");
- return;
- }
- s->sigma = atof (fp->token->str);
- gts_file_next_token (fp);
}
static void gfs_source_tension_write (GtsObject * o, FILE * fp)
{
- GfsSourceTension * s = GFS_SOURCE_TENSION (o);
(* GTS_OBJECT_CLASS (gfs_source_tension_class ())->parent_class->write) (o, fp);
- fprintf (fp, " %s %s %g", s->c->name, s->k->name, s->sigma);
-}
-
-typedef struct {
- gdouble amin, amax;
- guint depth;
- GfsFunction * alpha;
- GfsVariable * c;
-} StabilityParams;
-
-static void min_max_alpha (FttCell * cell, StabilityParams * p)
-{
- guint level = ftt_cell_level (cell);
- if (level > p->depth &&
- GFS_VARIABLE (cell, p->c->i) > 1e-3 &&
- GFS_VARIABLE (cell, p->c->i) < 1. - 1.e-3)
- p->depth = level;
- if (p->alpha) {
- gdouble a = gfs_function_value (p->alpha, cell);
- if (a < p->amin) p->amin = a;
- if (a > p->amax) p->amax = a;
- }
-}
-
-static gdouble gfs_source_tension_stability (GfsSourceGeneric * s,
- GfsSimulation * sim)
-{
- GfsSourceTension * t = GFS_SOURCE_TENSION (s);
- gdouble h;
- StabilityParams p = { G_MAXDOUBLE, -G_MAXDOUBLE, 0 };
-
- p.alpha = sim->physical_params.alpha;
- p.c = t->c;
- gfs_domain_cell_traverse (GFS_DOMAIN (sim), FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) min_max_alpha, &p);
- h = ftt_level_size (p.depth);
- if (p.alpha) {
- gdouble rhom = (1./p.amin + 1./p.amax)/2.;
- return sqrt (rhom*h*h*h/(2.*M_PI*t->sigma));
- }
- else
- return sqrt (h*h*h/(2.*M_PI*t->sigma));
+ fprintf (fp, " %s", GFS_SOURCE_TENSION (o)->k->name);
}
static void gfs_source_tension_class_init (GfsSourceGenericClass * klass)
{
GTS_OBJECT_CLASS (klass)->read = gfs_source_tension_read;
GTS_OBJECT_CLASS (klass)->write = gfs_source_tension_write;
- klass->stability = gfs_source_tension_stability;
}
GfsSourceGenericClass * gfs_source_tension_class (void)
@@ -272,7 +294,7 @@ GfsSourceGenericClass * gfs_source_tension_class (void)
(GtsArgGetFunc) NULL
};
klass =
- gts_object_class_new (GTS_OBJECT_CLASS (gfs_source_velocity_class ()),
+ gts_object_class_new (GTS_OBJECT_CLASS (gfs_source_tension_generic_class ()),
&gfs_source_tension_info);
}
diff --git a/src/tension.h b/src/tension.h
index 8299356..79f92a8 100644
--- a/src/tension.h
+++ b/src/tension.h
@@ -26,18 +26,38 @@ extern "C" {
#include "source.h"
+/* GfsSourceTensionGeneric: Header */
+
+typedef struct _GfsSourceTensionGeneric GfsSourceTensionGeneric;
+
+struct _GfsSourceTensionGeneric {
+ /*< private >*/
+ GfsSourceVelocity parent;
+
+ /*< public >*/
+ GfsVariable * c;
+ gdouble sigma;
+};
+
+#define GFS_SOURCE_TENSION_GENERIC(obj) GTS_OBJECT_CAST (obj,\
+ GfsSourceTensionGeneric,\
+ gfs_source_tension_generic_class ())
+#define GFS_IS_SOURCE_TENSION_GENERIC(obj) (gts_object_is_from_class (obj,\
+ gfs_source_tension_generic_class ()))
+
+GfsSourceGenericClass * gfs_source_tension_generic_class (void);
+
/* GfsSourceTensionCSS: Header */
typedef struct _GfsSourceTensionCSS GfsSourceTensionCSS;
struct _GfsSourceTensionCSS {
/*< private >*/
- GfsSourceVelocity parent;
+ GfsSourceTensionGeneric parent;
GfsVariable * g[3];
/*< public >*/
- GfsVariable * c, * t[FTT_DIMENSION];
- gdouble sigma;
+ GfsVariable * t[FTT_DIMENSION];
};
#define GFS_SOURCE_TENSION_CSS(obj) GTS_OBJECT_CAST (obj,\
@@ -54,11 +74,10 @@ typedef struct _GfsSourceTension GfsSourceTension;
struct _GfsSourceTension {
/*< private >*/
- GfsSourceVelocity parent;
+ GfsSourceTensionGeneric parent;
/*< public >*/
- GfsVariable * c, * k;
- gdouble sigma;
+ GfsVariable * k;
};
#define GFS_SOURCE_TENSION(obj) GTS_OBJECT_CAST (obj,\
diff --git a/src/timestep.c b/src/timestep.c
index 2b83404..c75bf1b 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -214,7 +214,7 @@ void gfs_mac_projection (GfsDomain * domain,
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,
+ GFS_SOURCE_TENSION_GENERIC (s)->c,
NULL, apar->dt);
}
/* Initialize face coefficients */
@@ -382,7 +382,7 @@ void gfs_approximate_projection (GfsDomain * domain,
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,
+ GFS_SOURCE_TENSION_GENERIC (s)->c,
g, apar->dt);
gfs_correct_centered_velocities (domain, FTT_DIMENSION, g, apar->dt);
}
diff --git a/src/vof.c b/src/vof.c
index b81cfcf..3078a6e 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -374,6 +374,47 @@ gdouble gfs_plane_alpha (FttVector * m, gdouble c)
}
#endif /* 3D */
+/**
+ * gfs_youngs_gradient:
+ * @cell: a #FttCell.
+ * @v: a #GfsVariable.
+ * @g: a #FttVector.
+ *
+ * Fills @g with the Youngs-averaged gradients of @v
+ * normalised by the size of @cell.
+ */
+void gfs_youngs_gradient (FttCell * cell, GfsVariable * v, FttVector * g)
+{
+ static FttDirection d[(FTT_DIMENSION - 1)*4][FTT_DIMENSION] = {
+#if FTT_2D
+ {FTT_RIGHT, FTT_TOP}, {FTT_LEFT, FTT_TOP}, {FTT_LEFT, FTT_BOTTOM}, {FTT_RIGHT, FTT_BOTTOM}
+#else /* 3D */
+ {FTT_RIGHT, FTT_TOP, FTT_FRONT}, {FTT_LEFT, FTT_TOP, FTT_FRONT},
+ {FTT_LEFT, FTT_BOTTOM, FTT_FRONT}, {FTT_RIGHT, FTT_BOTTOM, FTT_FRONT},
+ {FTT_RIGHT, FTT_TOP, FTT_BACK}, {FTT_LEFT, FTT_TOP, FTT_BACK},
+ {FTT_LEFT, FTT_BOTTOM, FTT_BACK}, {FTT_RIGHT, FTT_BOTTOM, FTT_BACK},
+#endif /* 3D */
+ };
+ gdouble u[(FTT_DIMENSION - 1)*4];
+ guint i;
+
+ g_return_if_fail (cell != NULL);
+ g_return_if_fail (v != NULL);
+ g_return_if_fail (g != NULL);
+
+ for (i = 0; i < (FTT_DIMENSION - 1)*4; i++)
+ u[i] = gfs_cell_corner_value (cell, d[i], v, -1);
+
+#if FTT_2D
+ g->x = (u[0] + u[3] - u[1] - u[2])/2.;
+ g->y = (u[0] + u[1] - u[2] - u[3])/2.;
+#else /* 3D */
+ g->x = (u[0] + u[3] + u[4] + u[7] - u[1] - u[2] - u[5] - u[6])/4.;
+ g->y = (u[0] + u[1] + u[4] + u[5] - u[2] - u[3] - u[6] - u[7])/4.;
+ g->z = (u[0] + u[1] + u[2] + u[3] - u[4] - u[5] - u[6] - u[7])/4.;
+#endif /* 3D */
+}
+
static FttCell * domain_and_boundary_locate (GfsDomain * domain, FttVector p, guint level)
{
FttCell * cell = gfs_domain_locate (domain, p, level);
@@ -428,8 +469,8 @@ static void stencil (FttCell * cell, GfsVariable * v, gdouble f[3][3])
* @v: a #GfsVariable.
* @n: a #FttVector.
*
- * Fills @n with the Youngs-averaged gradients of @v
- * normalised by the size of @cell.
+ * Fills @n with the components of the normal to the interface defined
+ * by volume fraction @v. Youngs' method is used.
*/
void gfs_youngs_normal (FttCell * cell, GfsVariable * v, FttVector * n)
{
diff --git a/src/vof.h b/src/vof.h
index 897660d..b359d97 100644
--- a/src/vof.h
+++ b/src/vof.h
@@ -49,6 +49,9 @@ void gfs_plane_center (FttVector * m,
gdouble a,
FttVector * p);
#endif /* 3D */
+void gfs_youngs_gradient (FttCell * cell,
+ GfsVariable * v,
+ FttVector * g);
void gfs_youngs_normal (FttCell * cell,
GfsVariable * v,
FttVector * n);
diff --git a/test/capwave/capwave.gfs b/test/capwave/capwave.gfs
index 49dc601..a6f51a4 100644
--- a/test/capwave/capwave.gfs
+++ b/test/capwave/capwave.gfs
@@ -49,7 +49,7 @@
Refine LEVEL
VariableTracer {} T { scheme = vof }
VariableCurvature {} K T
- SourceTension {} T K 1
+ SourceTension {} T 1 K
AdvectionParams { scheme = none }
SourceDiffusion {} U 0.0182571749236
SourceDiffusion {} V 0.0182571749236
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list