[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:41 UTC 2009
The following commit has been merged in the upstream branch:
commit 5e7820ba236219c8abc06a32f226251dd703100c
Author: Stephane Popinet <popinet at users.sf.net>
Date: Fri Aug 1 12:48:18 2008 +1000
Reimplementation of SourceVicosityExplicit
Also works for axisymmetric flows (but for 3D flows yet).
darcs-hash:20080801024818-d4795-62452a1262d6da8e29b0c10ef49769087c0b529d.gz
diff --git a/src/domain.h b/src/domain.h
index cf01164..4570812 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -72,9 +72,10 @@ struct _GfsDomain {
struct _GfsDomainClass {
GtsWGraphClass parent_class;
- void (* post_read) (GfsDomain *, GtsFile * fp);
- gdouble (* face_fraction) (const GfsDomain *, const FttCellFace *);
- gdouble (* cell_fraction) (const GfsDomain *, const FttCell *);
+ void (* post_read) (GfsDomain *, GtsFile * fp);
+ gdouble (* face_map) (const GfsDomain *, const FttCellFace *);
+ gdouble (* cell_map) (const GfsDomain *, const FttCell *);
+ gdouble (* solid_map) (const GfsDomain *, const FttCell *);
};
#define GFS_DOMAIN(obj) GTS_OBJECT_CAST (obj,\
@@ -97,6 +98,8 @@ void gfs_domain_cell_traverse (GfsDomain * domain,
gint max_depth,
FttCellTraverseFunc func,
gpointer data);
+#define gfs_domain_traverse_leaves(d,f,data) (gfs_domain_cell_traverse(d, \
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1, f,data))
void gfs_domain_cell_traverse_condition (GfsDomain * domain,
FttTraverseType order,
FttTraverseFlags flags,
@@ -305,8 +308,8 @@ static inline
gdouble gfs_domain_face_fraction (const GfsDomain * domain, const FttCellFace * face)
{
gdouble f = GFS_FACE_FRACTION (face);
- if (GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->face_fraction)
- f *= (* GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->face_fraction) (domain, face);
+ if (GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->face_map)
+ f *= (* GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->face_map) (domain, face);
return f;
}
@@ -322,11 +325,11 @@ static inline
gdouble gfs_domain_face_fraction_right (const GfsDomain * domain, const FttCellFace * face)
{
gdouble f = GFS_FACE_FRACTION_RIGHT (face);
- if (GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->face_fraction) {
+ if (GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->face_map) {
FttCellFace face1;
face1.cell = face->neighbor;
face1.d = FTT_OPPOSITE_DIRECTION (face->d);
- f *= (* GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->face_fraction) (domain, &face1);
+ f *= (* GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->face_map) (domain, &face1);
}
return f;
}
@@ -343,11 +346,27 @@ static inline
gdouble gfs_domain_cell_fraction (const GfsDomain * domain, const FttCell * cell)
{
gdouble a = GFS_IS_MIXED (cell) ? GFS_STATE (cell)->solid->a : 1.;
- if (GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->cell_fraction)
- a *= (* GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->cell_fraction) (domain, cell);
+ if (GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->cell_map)
+ a *= (* GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->cell_map) (domain, cell);
return a;
}
+/**
+ * gfs_domain_solid_map:
+ * @domain; a #GfsDomain.
+ * @cell: a mixed #FttCell.
+ *
+ * Returns: the coordinate mapping at the center of area of the solid
+ * surface contained within @cell.
+ */
+static inline
+gdouble gfs_domain_solid_map (const GfsDomain * domain, const FttCell * cell)
+{
+ if (GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->solid_map)
+ return (* GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->solid_map) (domain, cell);
+ return 1.;
+}
+
#ifdef __cplusplus
}
#endif /* __cplusplus */
diff --git a/src/fluid.c b/src/fluid.c
index fcccbfd..3f152b5 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -655,6 +655,7 @@ static void face_weighted_gradient (const FttCellFace * face,
g->b += w*(gcf.a*GFS_VARIABLE (f.cell, v) - gcf.c);
}
if (dimension > 2) {
+ /* fixme??? */
g->a /= n/2.;
g->b /= n/2.;
}
@@ -1808,6 +1809,53 @@ gdouble gfs_face_interpolated_value (const FttCellFace * face,
}
/**
+ * gfs_face_weighted_interpolated_value:
+ * @face: a #FttFace.
+ * @v: a #GfsVariable index.
+ *
+ * Computes the value of variable @v on the @face weighted by the
+ * value of the @v field of the face state vector using interpolation
+ * from the cell-centered values. The value returned is second order
+ * accurate in space and conservative, in the sense that values at a
+ * coarse/fine cell boundary are consistent.
+ *
+ * Returns: the weighted value of variable @v on the face.
+ */
+gdouble gfs_face_weighted_interpolated_value (const FttCellFace * face,
+ guint v)
+{
+ g_return_val_if_fail (face != NULL, 0.);
+
+ if (face->neighbor) {
+ if (FTT_CELL_IS_LEAF (face->neighbor)) {
+ gdouble w = GFS_STATE (face->cell)->f[face->d].v, x1 = 1., v1;
+ v1 = neighbor_value (face, v, &x1);
+ return w*((x1 - 0.5)*GFS_VARIABLE (face->cell, v) + 0.5*v1)/x1;
+ }
+ else {
+ /* neighbor is at a deeper level */
+ FttCellChildren children;
+ FttCellFace f;
+ gdouble val = 0.;
+ guint i, n;
+
+ f.d = FTT_OPPOSITE_DIRECTION (face->d);
+ n = ftt_cell_children_direction (face->neighbor, f.d, &children);
+ f.neighbor = face->cell;
+ for (i = 0; i < n; i++)
+ if ((f.cell = children.c[i])) {
+ gdouble w = GFS_STATE (f.cell)->f[f.d].v, x1 = 1., v1;
+ v1 = neighbor_value (&f, v, &x1);
+ val += w*v1;
+ }
+ return val/n;
+ }
+ }
+ else
+ return GFS_STATE (face->cell)->f[face->d].v*GFS_VARIABLE (face->cell, v);
+}
+
+/**
* gfs_normal_divergence:
* @cell: a #FttCell.
* @v: a #GfsVariable.
diff --git a/src/fluid.h b/src/fluid.h
index 44d9c4f..727a759 100644
--- a/src/fluid.h
+++ b/src/fluid.h
@@ -102,6 +102,8 @@ void gfs_cell_coarse_fine (FttCell * cell,
GfsVariable * v);
gdouble gfs_face_interpolated_value (const FttCellFace * face,
guint v);
+gdouble gfs_face_weighted_interpolated_value (const FttCellFace * face,
+ guint v);
typedef gdouble (* GfsCenterGradient) (FttCell * cell,
FttComponent c,
guint v);
diff --git a/src/poisson.c b/src/poisson.c
index 38279d7..af52f4a 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -752,16 +752,17 @@ typedef struct {
gdouble dt;
GfsVariable * dia;
GfsFunction * alpha;
+ GfsDomain * domain;
} DiffusionCoeff;
static void diffusion_coef (FttCellFace * face, DiffusionCoeff * c)
{
- GfsStateVector * s = GFS_STATE (face->cell);
- gdouble v = c->lambda2[face->d/2]*c->dt*gfs_source_diffusion_face (c->d, face);
+ gdouble v =
+ c->lambda2[face->d/2]*c->dt*
+ gfs_source_diffusion_face (c->d, face)*
+ gfs_domain_face_fraction (c->domain, face);
- if (GFS_IS_MIXED (face->cell))
- v *= s->solid->s[face->d];
- s->f[face->d].v = v;
+ GFS_STATE (face->cell)->f[face->d].v = v;
switch (ftt_face_type (face)) {
case FTT_FINE_FINE:
@@ -780,15 +781,18 @@ static void diffusion_mixed_coef (FttCell * cell, DiffusionCoeff * c)
{
reset_coeff (cell);
if (GFS_IS_MIXED (cell))
- GFS_STATE (cell)->solid->v = c->dt*gfs_source_diffusion_cell (c->d, cell);
- GFS_VARIABLE (cell, c->dia->i) = c->alpha ? 1./gfs_function_value (c->alpha, cell) : 1.;
- if (GFS_VARIABLE (cell, c->dia->i) <= 0.) {
- FttVector p;
- ftt_cell_pos (cell, &p);
- g_log (G_LOG_DOMAIN, G_LOG_LEVEL_ERROR,
- "density is negative (%g) at cell (%g,%g,%g).\n"
- "Please check your definition of alpha.",
- GFS_VARIABLE (cell, c->dia->i), p.x, p.y, p.z);
+ GFS_STATE (cell)->solid->v =
+ c->dt*gfs_domain_solid_map (c->domain, cell)*gfs_source_diffusion_cell (c->d, cell);
+ if (c->dia) {
+ GFS_VALUE (cell, c->dia) = c->alpha ? 1./gfs_function_value (c->alpha, cell) : 1.;
+ if (GFS_VALUE (cell, c->dia) <= 0.) {
+ FttVector p;
+ ftt_cell_pos (cell, &p);
+ g_log (G_LOG_DOMAIN, G_LOG_LEVEL_ERROR,
+ "density is negative (%g) at cell (%g,%g,%g).\n"
+ "Please check your definition of alpha.",
+ GFS_VALUE (cell, c->dia), p.x, p.y, p.z);
+ }
}
}
@@ -815,7 +819,6 @@ void gfs_diffusion_coefficients (GfsDomain * domain,
g_return_if_fail (domain != NULL);
g_return_if_fail (d != NULL);
- g_return_if_fail (dia != NULL);
g_return_if_fail (beta >= 0.5 && beta <= 1.);
for (i = 0; i < FTT_DIMENSION; i++) {
@@ -827,6 +830,7 @@ void gfs_diffusion_coefficients (GfsDomain * domain,
coef.dt = beta*dt;
coef.dia = dia;
coef.alpha = alpha;
+ coef.domain = domain;
gfs_domain_cell_traverse (domain,
FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) diffusion_mixed_coef, &coef);
diff --git a/src/simulation.c b/src/simulation.c
index fb42aef..e61c192 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -1792,25 +1792,32 @@ static void axi_read (GtsObject ** object, GtsFile * fp)
GFS_DOMAIN (*object)->refpos.y = 0.5;
}
-static gdouble axi_face_fraction (const GfsDomain * domain, const FttCellFace * face)
+static gdouble axi_face_map (const GfsDomain * domain, const FttCellFace * face)
{
FttVector p;
ftt_face_pos (face, &p);
return p.y;
}
-static gdouble axi_cell_fraction (const GfsDomain * domain, const FttCell * cell)
+static gdouble axi_cell_map (const GfsDomain * domain, const FttCell * cell)
{
FttVector p;
gfs_cell_cm (cell, &p);
return p.y;
}
+static gdouble axi_solid_map (const GfsDomain * domain, const FttCell * cell)
+{
+ g_assert (GFS_IS_MIXED (cell));
+ return GFS_STATE (cell)->solid->ca.y;
+}
+
static void axi_class_init (GfsSimulationClass * klass)
{
GTS_OBJECT_CLASS (klass)->read = axi_read;
- GFS_DOMAIN_CLASS (klass)->face_fraction = axi_face_fraction;
- GFS_DOMAIN_CLASS (klass)->cell_fraction = axi_cell_fraction;
+ GFS_DOMAIN_CLASS (klass)->face_map = axi_face_map;
+ GFS_DOMAIN_CLASS (klass)->cell_map = axi_cell_map;
+ GFS_DOMAIN_CLASS (klass)->solid_map = axi_solid_map;
}
GfsSimulationClass * gfs_axi_class (void)
diff --git a/src/simulation.h b/src/simulation.h
index 493f628..e14224c 100644
--- a/src/simulation.h
+++ b/src/simulation.h
@@ -137,6 +137,8 @@ GfsSimulationClass * gfs_poisson_class (void);
/* GfsAxi: Header */
+#define GFS_IS_AXI(obj) (gts_object_is_from_class (obj, gfs_axi_class ()))
+
GfsSimulationClass * gfs_axi_class (void);
#ifdef __cplusplus
diff --git a/src/source.c b/src/source.c
index 5a08cc0..94b02e7 100644
--- a/src/source.c
+++ b/src/source.c
@@ -53,23 +53,24 @@ gdouble gfs_variable_mac_source (GfsVariable * v, FttCell * cell)
return sum;
}
-static void add_sources (FttCell * cell, gpointer * data)
+typedef struct {
+ GfsVariable * v, * sv;
+ gdouble dt;
+} SourcePar;
+
+static void add_sources (FttCell * cell, SourcePar * p)
{
- GfsVariable * v = data[0];
- GfsVariable * sv = data[1];
- gdouble * dt = data[2];
- GSList * i = GTS_SLIST_CONTAINER (v->sources)->items;
+ GSList * i = GTS_SLIST_CONTAINER (p->v->sources)->items;
gdouble sum = 0;
- i = GTS_SLIST_CONTAINER (v->sources)->items;
while (i) {
GfsSourceGeneric * s = i->data;
if (s->centered_value)
- sum += (* s->centered_value) (s, cell, v);
+ sum += (* s->centered_value) (s, cell, p->v);
i = i->next;
}
- GFS_VARIABLE (cell, sv->i) += (*dt)*sum;
+ GFS_VALUE (cell, p->sv) += p->dt*sum;
}
/**
@@ -91,15 +92,45 @@ void gfs_domain_variable_centered_sources (GfsDomain * domain,
g_return_if_fail (sv != NULL);
if (v->sources) {
- gpointer data[3];
-
- data[0] = v;
- data[1] = sv;
- data[2] = &dt;
+ SourcePar p;
+ p.v = v;
+ p.sv = sv;
+ p.dt = dt;
gfs_domain_cell_traverse (domain,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) add_sources, data);
+ (FttCellTraverseFunc) add_sources, &p);
+ }
+}
+
+/**
+ * gfs_domain_variable_fluxes:
+ * @domain: a #GfsDomain.
+ * @v: a #GfsVariable.
+ * @dt: the timestep.
+ *
+ * Returns: a new temporary variable containing the fluxes or %NULL.
+ */
+GfsVariable * gfs_domain_variable_fluxes (GfsDomain * domain,
+ GfsVariable * v,
+ gdouble dt)
+{
+ GfsVariable * sv = NULL;
+
+ g_return_val_if_fail (domain != NULL, NULL);
+ g_return_val_if_fail (v != NULL, NULL);
+
+ GSList * i = GTS_SLIST_CONTAINER (v->sources)->items;
+ while (i) {
+ if (GFS_SOURCE_GENERIC (i->data)->flux) {
+ if (sv == NULL) {
+ sv = gfs_temporary_variable (domain);
+ gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) gfs_cell_reset, sv);
+ }
+ (* GFS_SOURCE_GENERIC (i->data)->flux) (i->data, domain, v, sv, dt);
+ }
+ i = i->next;
}
+ return sv;
}
/* GfsSourceGeneric: Object */
@@ -729,7 +760,6 @@ gdouble gfs_source_diffusion_cell (GfsSourceDiffusion * d, FttCell * cell)
{
g_return_val_if_fail (d != NULL, 0.);
g_return_val_if_fail (cell != NULL, 0.);
- g_return_val_if_fail (GFS_IS_MIXED (cell), 0.);
return gfs_diffusion_cell (d->D, cell);
}
@@ -875,6 +905,12 @@ static void source_viscosity_read (GtsObject ** o, GtsFile * fp)
source = GFS_SOURCE_VISCOSITY (*o);
domain = GFS_DOMAIN (gfs_object_simulation (source));
+ if (GFS_IS_AXI (domain) && !GFS_IS_SOURCE_VISCOSITY_EXPLICIT (source)) {
+ gts_file_error (fp,
+ "GfsSourceViscosity does not work for axisymmetric flows (yet)\n"
+ "Use GfsSourceViscosityExplicit instead\n");
+ return;
+ }
if (!(source->v = gfs_domain_velocity (domain))) {
gts_file_error (fp, "cannot find velocity components");
return;
@@ -964,6 +1000,31 @@ GfsSourceGenericClass * gfs_source_viscosity_class (void)
/* GfsSourceViscosityExplicit: Object */
+static void gfs_source_viscosity_explicit_destroy (GtsObject * o)
+{
+ GfsVariable ** v = GFS_SOURCE_VISCOSITY_EXPLICIT (o)->v;
+
+ if (v[0]) {
+ FttComponent c;
+ for (c = 0; c < FTT_DIMENSION; c++)
+ gts_object_destroy (GTS_OBJECT (v[c]));
+ }
+
+ (* GTS_OBJECT_CLASS (gfs_source_viscosity_explicit_class ())->parent_class->destroy) (o);
+}
+
+static void gfs_source_viscosity_explicit_read (GtsObject ** o, GtsFile * fp)
+{
+ (* GTS_OBJECT_CLASS (gfs_source_viscosity_explicit_class ())->parent_class->read) (o, fp);
+ if (fp->type == GTS_ERROR)
+ return;
+
+ FttComponent c;
+ for (c = 0; c < FTT_DIMENSION; c++)
+ GFS_SOURCE_VISCOSITY_EXPLICIT (*o)->v[c] =
+ gfs_temporary_variable (GFS_DOMAIN (gfs_object_simulation (*o)));
+}
+
static gdouble source_viscosity_stability (GfsSourceGeneric * s,
GfsSimulation * sim)
{
@@ -974,17 +1035,151 @@ static gdouble source_viscosity_stability (GfsSourceGeneric * s,
par.alpha = sim->physical_params.alpha;
gfs_domain_cell_traverse (GFS_DOMAIN (sim), FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) cell_diffusion_stability, &par);
- return par.dtmax;
+ return 0.1*par.dtmax;
+}
+
+static void compute_transverse (FttCell * cell, GfsSourceViscosityExplicit * s)
+{
+ GfsVariable ** v = GFS_SOURCE_VISCOSITY (s)->v;
+ gdouble h = ftt_cell_size (cell);
+ FttComponent c;
+ for (c = 0; c < FTT_DIMENSION; c++)
+ GFS_VALUE (cell, s->v[c]) = gfs_center_gradient (cell, (c + 1) % FTT_DIMENSION, v[c]->i)/h;
+}
+
+static gboolean gfs_source_viscosity_explicit_event (GfsEvent * event, GfsSimulation * sim)
+{
+ if ((* GFS_EVENT_CLASS (GTS_OBJECT_CLASS (gfs_source_viscosity_explicit_class ())->parent_class)->event) (event, sim)) {
+ FttComponent c;
+
+ gfs_domain_cell_traverse (GFS_DOMAIN (sim), FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) compute_transverse, event);
+ for (c = 0; c < FTT_DIMENSION; c++)
+ gfs_domain_bc (GFS_DOMAIN (sim), FTT_TRAVERSE_LEAFS, -1,
+ GFS_SOURCE_VISCOSITY_EXPLICIT (event)->v[c]);
+ return TRUE;
+ }
+ return FALSE;
}
static void source_viscosity_explicit_class_init (GfsSourceGenericClass * klass)
{
+ GTS_OBJECT_CLASS (klass)->destroy = gfs_source_viscosity_explicit_destroy;
+ GTS_OBJECT_CLASS (klass)->read = gfs_source_viscosity_explicit_read;
+ GFS_EVENT_CLASS (klass)->event = gfs_source_viscosity_explicit_event;
klass->stability = source_viscosity_stability;
}
+typedef struct {
+ GfsSourceGeneric * s;
+ GfsVariable * v, * sv;
+ gdouble dt;
+} FluxPar;
+
+/* Defining this will use the divergence-free condition to decouple
+ the diffusion equations for each component. This only works of
+ course for constant viscosity and does not work for axisymmetric flows. */
+/* #define NOTRANSVERSE */
+
+static void add_viscosity_explicit_flux (FttCell * cell, FluxPar * p)
+{
+ FttCellFace f;
+ FttCellNeighbors n;
+ GfsGradient g = { 0., 0. };
+ gdouble v0, h;
+
+ if (GFS_IS_MIXED (cell)) {
+ if (((cell)->flags & GFS_FLAG_DIRICHLET) != 0)
+ g.b = gfs_cell_dirichlet_gradient_flux (cell, p->v->i, -1., 0.);
+ }
+
+ v0 = GFS_VARIABLE (cell, p->v->i);
+ f.cell = cell;
+ ftt_cell_neighbors (cell, &n);
+ for (f.d = 0; f.d < FTT_NEIGHBORS; f.d++) {
+ GfsGradient e;
+
+ f.neighbor = n.c[f.d];
+ gfs_face_gradient_flux (&f, &e, p->v->i, -1);
+#ifndef NOTRANSVERSE
+ if (f.d/2 == p->v->component) {
+ e.a *= 2.;
+ e.b *= 2.;
+ }
+#endif
+ g.a += e.a;
+ g.b += e.b;
+ }
+ h = ftt_cell_size (cell);
+
+ gdouble transverse = 0.;
+#ifndef NOTRANSVERSE
+#if FTT_2D
+ FttComponent ortho = (p->v->component + 1) % FTT_DIMENSION;
+
+ for (f.d = 2*ortho; f.d <= 2*ortho + 1; f.d++) {
+ f.neighbor = n.c[f.d];
+#if 0
+ transverse += (FTT_FACE_DIRECT (&f) ? 1. : -1.)*
+ p->dt*
+ gfs_source_diffusion_face (GFS_SOURCE_DIFFUSION (p->s), &f)*
+ gfs_domain_face_fraction (p->v->domain, &f)*
+ gfs_face_interpolated_value (&f, GFS_SOURCE_VISCOSITY_EXPLICIT (p->s)->v[ortho]->i);
+#elif 0
+ transverse += (FTT_FACE_DIRECT (&f) ? 1. : -1.)*GFS_STATE (cell)->f[f.d].v*
+ gfs_face_interpolated_value (&f, GFS_SOURCE_VISCOSITY_EXPLICIT (p->s)->v[ortho]->i);
+#else
+ transverse += (FTT_FACE_DIRECT (&f) ? 1. : -1.)*
+ gfs_face_weighted_interpolated_value (&f, GFS_SOURCE_VISCOSITY_EXPLICIT (p->s)->v[ortho]->i);
+#endif
+ }
+#else
+ g_assert_not_implemented ();
+#endif
+#endif
+
+ GfsFunction * alpha = gfs_object_simulation (p->s)->physical_params.alpha;
+ GFS_VALUE (cell, p->sv) += (alpha ? gfs_function_value (alpha, cell) : 1.)*
+ ((g.b - g.a*v0)/h + transverse)/h;
+}
+
+static void add_axisymmetric_term (FttCell * cell, FluxPar * p)
+{
+ GfsFunction * alpha = gfs_object_simulation (p->s)->physical_params.alpha;
+ gdouble a = GFS_IS_MIXED (cell) ? GFS_STATE (cell)->solid->a : 1.;
+ GFS_VALUE (cell, p->sv) -=
+ (alpha ? gfs_function_value (alpha, cell) : 1.)*
+ 2.*gfs_source_diffusion_cell (GFS_SOURCE_DIFFUSION (p->s), cell)*
+ GFS_VALUE (cell, p->v)*
+ a*a/gfs_domain_cell_fraction (p->v->domain, cell)*
+ p->dt;
+}
+
+static void source_viscosity_explicit_flux (GfsSourceGeneric * s,
+ GfsDomain * domain,
+ GfsVariable * v, GfsVariable * sv,
+ gdouble dt)
+{
+ FluxPar p;
+
+ gfs_diffusion_coefficients (domain, GFS_SOURCE_DIFFUSION (s), dt, NULL, NULL, 1.);
+ gfs_domain_surface_bc (domain, v);
+ p.s = s;
+ p.v = v;
+ p.sv = sv;
+ p.dt = dt;
+ gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) add_viscosity_explicit_flux, &p);
+#if 1
+ if (GFS_IS_AXI (domain) && v->component == FTT_Y)
+ gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) add_axisymmetric_term, &p);
+#endif
+}
+
static void source_viscosity_explicit_init (GfsSourceGeneric * s)
{
- s->mac_value = s->centered_value = source_viscosity_value;
+ s->mac_value = NULL;
+ s->centered_value = NULL;
+ s->flux = source_viscosity_explicit_flux;
}
GfsSourceGenericClass * gfs_source_viscosity_explicit_class (void)
@@ -994,7 +1189,7 @@ GfsSourceGenericClass * gfs_source_viscosity_explicit_class (void)
if (klass == NULL) {
GtsObjectClassInfo source_viscosity_explicit_info = {
"GfsSourceViscosityExplicit",
- sizeof (GfsSourceViscosity),
+ sizeof (GfsSourceViscosityExplicit),
sizeof (GfsSourceGenericClass),
(GtsObjectClassInitFunc) source_viscosity_explicit_class_init,
(GtsObjectInitFunc) source_viscosity_explicit_init,
diff --git a/src/source.h b/src/source.h
index 9aa5767..b373a49 100644
--- a/src/source.h
+++ b/src/source.h
@@ -32,6 +32,10 @@ void gfs_domain_variable_centered_sources (GfsDomain * domain,
GfsVariable * v,
GfsVariable * sv,
gdouble dt);
+GfsVariable * gfs_domain_variable_fluxes (GfsDomain * domain,
+ GfsVariable * v,
+ gdouble dt);
+
/* GfsSourceGeneric: Header */
typedef struct _GfsSourceGeneric GfsSourceGeneric;
@@ -45,6 +49,9 @@ struct _GfsSourceGeneric {
gdouble (* mac_value) (GfsSourceGeneric *, FttCell *, GfsVariable *);
gdouble (* centered_value) (GfsSourceGeneric *, FttCell *, GfsVariable *);
gdouble (* face_value) (GfsSourceGeneric *, FttCellFace *, GfsVariable *);
+ void (* flux) (GfsSourceGeneric *, GfsDomain *,
+ GfsVariable *, GfsVariable *,
+ gdouble);
};
typedef struct _GfsSourceGenericClass GfsSourceGenericClass;
@@ -255,6 +262,19 @@ GfsSourceGenericClass * gfs_source_viscosity_class (void);
/* GfsSourceViscosityExplicit: Header */
+typedef struct _GfsSourceViscosityExplicit GfsSourceViscosityExplicit;
+
+struct _GfsSourceViscosityExplicit {
+ /*< private >*/
+ GfsSourceViscosity parent;
+
+ /*< public >*/
+ GfsVariable * v[FTT_DIMENSION];
+};
+
+#define GFS_SOURCE_VISCOSITY_EXPLICIT(obj) GTS_OBJECT_CAST (obj,\
+ GfsSourceViscosityExplicit,\
+ gfs_source_viscosity_class ())
#define GFS_IS_SOURCE_VISCOSITY_EXPLICIT(obj) (gts_object_is_from_class (obj,\
gfs_source_viscosity_explicit_class ()))
diff --git a/src/timestep.c b/src/timestep.c
index cec3e65..ac17ce5 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -701,6 +701,15 @@ static void variable_sources (GfsDomain * domain,
}
/* fixme: time should be set to t + dt/2 here for evaluation of
source terms in the call below */
+ par->fv = gfs_domain_variable_fluxes (domain, par->v, par->dt);
+ if (par->fv) {
+ GfsVariable * v = par->v;
+ par->v = sv;
+ gfs_domain_traverse_merged (domain, (GfsMergedTraverseFunc) gfs_advection_update, par);
+ par->v = v;
+ gts_object_destroy (GTS_OBJECT (par->fv));
+ par->fv = NULL;
+ }
gfs_domain_variable_centered_sources (domain, par->v, sv, par->dt);
}
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list