[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:43 UTC 2009
The following commit has been merged in the upstream branch:
commit e86db092773d659f85f14be6ff719029fed172e8
Author: Stephane Popinet <popinet at users.sf.net>
Date: Sat Sep 6 08:04:37 2008 +1000
Implicit viscosity now works for axisymmetric domains
The explicit viscosity scheme has also been simplified.
darcs-hash:20080905220437-d4795-bf5c2dcae85ca842de1e60193588705e5cd51b0f.gz
diff --git a/doc/figures/axi.tm b/doc/figures/axi.tm
index 80a46ba..9c46fcd 100644
--- a/doc/figures/axi.tm
+++ b/doc/figures/axi.tm
@@ -80,7 +80,7 @@
<tformat|<table|<row|<cell|>|<cell|c<rsup|i,j>*<frac|(v<rsub|z><rsup|\<star\>>)<rsup|i,j>-v<rsup|i,j><rsub|z>|\<Delta\>t>+(s<rsub|z>*v<rsup|2><rsub|z>)<rsup|i+1/2,j>-(s<rsub|z>*v<rsup|2><rsub|z>)<rsup|i-1/2,j>+(s<rsub|r>*v<rsub|r><rsup|>*v<rsub|z>)<rsup|i,j+1/2>-(s<rsub|r>*v<rsub|r><rsup|>*v<rsub|z>)<rsup|i,j-1/2>=>|<cell|>>|<row|<cell|>|<cell|(s<rsub|z>*S<rsub|zz>)<rsup|i+1/2,j>-(s<rsub|z>*S<rsub|zz>)<rsup|i-1/2,j>+(s<rsub|r>*S<rsub|zr>)<rsup|i,j+1/2>-(s<rsub|r>*S<rsub|zr>)<rsup|i,j-1/2>.>|<cell|>>>>
</eqnarray*>
- <subsubsection|Advection term>
+ <subsection|Advection term>
<\equation*>
u<rsup|n+1/2><rsub|d>=u<rsup|n>+<frac|h|2>*\<partial\><rsub|d>u<rsup|n>+<frac|\<Delta\>t|2>*\<partial\><rsub|t>u<rsup|n>+\<cal-O\>(h<rsup|2>,\<Delta\>t<rsup|2>),
@@ -95,12 +95,68 @@
which is the same as in the two-dimensional case.
- \;
+ <subsection|Implicit diffusion>
+
+ <subsubsection|Original implicit diffusion>
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|c*<frac|v<rsup|\<star\>>-v<rsup|n>|\<Delta\>t>-\<beta\>*\<nabla\>s*S<rsup|\<star\>><rsup|*>=src<rsup|n+1/2>+(1-\<beta\>)*\<nabla\>s*S<rsup|n>>|<cell|>>>>
+ </eqnarray*>
+
+ which can be rewritten
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|c*v<rsup|\<star\>>-\<nabla\><wide|s|^>*S<rsup|\<star\>>=c*v<rsup|n>+\<Delta\>t*src<rsup|n+1/2>+<frac|1-\<beta\>|\<beta\>>*\<nabla\><wide|s|^>*S<rsup|n>,>|<cell|>>>>
+ </eqnarray*>
+
+ with
+
+ <\equation*>
+ <wide|s|^>\<equiv\>\<beta\>*\<Delta\>t*s.
+ </equation*>
+
+ This finally gives
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|v<rsup|\<star\>>-<frac|1|c>\<nabla\><wide|s|^>*S<rsup|\<star\>>=v<rsup|n>+\<Delta\>t*<wide|src|^><rsup|n+1/2>+<frac|<wide|\<beta\>|^>|c>*\<nabla\><wide|s|^>*S<rsup|n>,>|<cell|>>>>
+ </eqnarray*>
+
+ with <math|><with|mode|math|<wide|\<beta\>|^>\<equiv\>(1-\<beta\>)/\<beta\>>.
+
+ <subsubsection|Axisymmetric version>
+
+ To account for the axisymmetric term this needs to be rewritten
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|c*<frac|v<rsup|\<star\>>-v<rsup|n>|\<Delta\>t>-\<beta\>*\<nabla\>s*S<rsup|\<star\>>+\<beta\>*a*2*\<nu\>*<frac|v<rsup|\<star\>>|r>=src<rsup|n+1/2>+(1-\<beta\>)*\<nabla\>s*S<rsup|n>-(1-\<beta\>)*a*2*\<nu\>*<frac|v<rsup|n>|r>>|<cell|>>>>
+ </eqnarray*>
+
+ which can be rewritten
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|c*v<rsup|\<star\>>*-\<nabla\><wide|s|^>*S<rsup|\<star\>>+2*\<beta\>*\<Delta\>t*a*\<nu\>*<frac|v<rsup|\<star\>>|r>=c*v<rsup|n>+\<Delta\>t*src<rsup|n+1/2>+<frac|1-\<beta\>|\<beta\>>*\<nabla\><wide|s|^>*S<rsup|n>-2*(1-\<beta\>)*\<Delta\>t*a*\<nu\>*<frac|v<rsup|n>|r>>|<cell|>>>>
+ </eqnarray*>
+
+ this finally gives
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|v<rsup|\<star\>>*-<frac|1|c>*\<nabla\><wide|s|^>*S<rsup|\<star\>>+2*\<beta\>*\<Delta\>t*\<nu\>*<frac|v<rsup|\<star\>>|r<rsup|2>>=v<rsup|n>+\<Delta\>t*<wide|src|^><rsup|n+1/2>+<frac|1-\<beta\>|\<beta\>*c>*\<nabla\><wide|s|^>*S<rsup|n>-2*(1-\<beta\>)*\<Delta\>t*\<nu\>*<frac|v<rsup|n>|r<rsup|2>>>|<cell|>>>>
+ </eqnarray*>
+
+ Posing <math|d\<equiv\>2*\<beta\>*\<Delta\>t*\<nu\>/r<rsup|2>> then gives
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|v<rsup|\<star\>>*(1+d)-<frac|1|c>*\<nabla\><wide|s|^>*S<rsup|\<star\>>=v<rsup|n>*<left|(>1-<wide|\<beta\>*|^>*d<right|)>+\<Delta\>t*<wide|src|^><rsup|n+1/2>+<frac|<wide|\<beta\>*|^>*|c>*\<nabla\><wide|s|^>*S<rsup|n>.>|<cell|>>>>
+ </eqnarray*>
</body>
<\references>
<\collection>
<associate|auto-1|<tuple|1|?>>
+ <associate|auto-2|<tuple|2|?>>
+ <associate|auto-3|<tuple|2.1|?>>
+ <associate|auto-4|<tuple|2.2|?>>
+ <associate|auto-5|<tuple|3|?>>
<associate|continuity|<tuple|3|?>>
<associate|momr|<tuple|1|?>>
<associate|momz|<tuple|2|?>>
diff --git a/src/domain.c b/src/domain.c
index c4989da..8f34bf5 100644
--- a/src/domain.c
+++ b/src/domain.c
@@ -2812,6 +2812,9 @@ void gfs_domain_solid_force (GfsDomain * domain,
g_return_if_fail (pm != NULL);
g_return_if_fail (vm != NULL);
+ if (GFS_IS_AXI (domain))
+ g_assert_not_implemented ();
+
pf->x = pf->y = pf->z = 0.;
pm->x = pm->y = pm->z = 0.;
data[0] = pf;
diff --git a/src/poisson.c b/src/poisson.c
index af52f4a..2968cc6 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -177,6 +177,8 @@ typedef struct {
guint u, rhs, dia, res;
gint maxlevel;
gdouble beta, omega;
+ FttComponent component;
+ guint axi;
} RelaxParams;
static void relax (FttCell * cell, RelaxParams * p)
@@ -750,7 +752,7 @@ typedef struct {
GfsSourceDiffusion * d;
gdouble lambda2[FTT_DIMENSION];
gdouble dt;
- GfsVariable * dia;
+ GfsVariable * rhoc, * axi;
GfsFunction * alpha;
GfsDomain * domain;
} DiffusionCoeff;
@@ -783,15 +785,23 @@ static void diffusion_mixed_coef (FttCell * cell, DiffusionCoeff * c)
if (GFS_IS_MIXED (cell))
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.) {
+ if (c->rhoc) {
+ gdouble rho = c->alpha ? 1./gfs_function_value (c->alpha, cell) : 1.;
+ if (rho <= 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);
+ rho, p.x, p.y, p.z);
+ }
+ gdouble f = gfs_domain_cell_fraction (c->domain, cell);
+ GFS_VALUE (cell, c->rhoc) = rho*f;
+
+ if (c->axi) {
+ FttVector p;
+ gfs_cell_cm (cell, &p);
+ GFS_VALUE (cell, c->axi) = 2.*c->dt*gfs_source_diffusion_cell (c->d, cell)/(rho*p.y*p.y);
}
}
}
@@ -801,7 +811,8 @@ static void diffusion_mixed_coef (FttCell * cell, DiffusionCoeff * c)
* @domain: a #GfsDomain.
* @d: a #GfsSourceDiffusion.
* @dt: the time-step.
- * @dia: where to store the diagonal weight.
+ * @rhoc: where to store the mass.
+ * @axi: where to store the axisymmetric term (or %NULL).
* @alpha: the inverse of density or %NULL.
* @beta: the implicitness parameter (0.5 Crank-Nicholson, 1. backward Euler).
*
@@ -810,7 +821,8 @@ static void diffusion_mixed_coef (FttCell * cell, DiffusionCoeff * c)
void gfs_diffusion_coefficients (GfsDomain * domain,
GfsSourceDiffusion * d,
gdouble dt,
- GfsVariable * dia,
+ GfsVariable * rhoc,
+ GfsVariable * axi,
GfsFunction * alpha,
gdouble beta)
{
@@ -828,9 +840,10 @@ void gfs_diffusion_coefficients (GfsDomain * domain,
}
coef.d = d;
coef.dt = beta*dt;
- coef.dia = dia;
+ coef.rhoc = rhoc;
coef.alpha = alpha;
coef.domain = domain;
+ coef.axi = axi;
gfs_domain_cell_traverse (domain,
FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) diffusion_mixed_coef, &coef);
@@ -843,24 +856,20 @@ void gfs_diffusion_coefficients (GfsDomain * domain,
NULL);
}
-
static void diffusion_rhs (FttCell * cell, RelaxParams * p)
{
- gdouble a, f, h, val;
+ gdouble f, h, val;
FttCellNeighbors neighbor;
FttCellFace face;
if (GFS_IS_MIXED (cell)) {
- a = GFS_STATE (cell)->solid->a*GFS_VARIABLE (cell, p->dia);
if (((cell)->flags & GFS_FLAG_DIRICHLET) != 0)
f = gfs_cell_dirichlet_gradient_flux (cell, p->u, -1, GFS_STATE (cell)->solid->fv);
else
f = GFS_STATE (cell)->solid->fv;
}
- else {
- a = GFS_VARIABLE (cell, p->dia);
+ else
f = 0.; /* Neumann condition by default */
- }
h = ftt_cell_size (cell);
val = GFS_VARIABLE (cell, p->u);
face.cell = cell;
@@ -870,9 +879,15 @@ static void diffusion_rhs (FttCell * cell, RelaxParams * p)
face.neighbor = neighbor.c[face.d];
gfs_face_gradient_flux (&face, &g, p->u, -1);
+ if (face.d/2 == p->component) {
+ g.a *= 2.;
+ g.b *= 2.;
+ }
f += g.b - g.a*val;
}
- GFS_VARIABLE (cell, p->rhs) += val + p->beta*f/(h*h*a);
+ GFS_VARIABLE (cell, p->rhs) += val + p->beta*f/(h*h*GFS_VARIABLE (cell, p->dia));
+ if (p->axi)
+ GFS_VARIABLE (cell, p->rhs) -= val*p->beta*GFS_VARIABLE (cell, p->axi);
}
/**
@@ -880,7 +895,8 @@ static void diffusion_rhs (FttCell * cell, RelaxParams * p)
* @domain: a #GfsDomain.
* @v: a #GfsVariable.
* @rhs: a #GfsVariable.
- * @dia: the diagonal weight.
+ * @rhoc: the mass.
+ * @axi: the axisymmetric term.
* @beta: the implicitness parameter (0.5 Crank-Nicholson, 1. backward Euler).
*
* Adds to the @rhs variable of @cell the right-hand side of the
@@ -890,7 +906,8 @@ static void diffusion_rhs (FttCell * cell, RelaxParams * p)
* gfs_diffusion_coefficients().
*/
void gfs_diffusion_rhs (GfsDomain * domain,
- GfsVariable * v, GfsVariable * rhs, GfsVariable * dia,
+ GfsVariable * v, GfsVariable * rhs,
+ GfsVariable * rhoc, GfsVariable * axi,
gdouble beta)
{
RelaxParams p;
@@ -898,13 +915,15 @@ void gfs_diffusion_rhs (GfsDomain * domain,
g_return_if_fail (domain != NULL);
g_return_if_fail (v != NULL);
g_return_if_fail (rhs != NULL);
- g_return_if_fail (dia != NULL);
+ g_return_if_fail (rhoc != NULL);
g_return_if_fail (beta >= 0.5 && beta <= 1.);
p.u = v->i;
p.rhs = rhs->i;
- p.dia = dia->i;
+ p.dia = rhoc->i;
p.beta = (1. - beta)/beta;
+ p.component = GFS_IS_AXI (domain) ? v->component : FTT_DIMENSION;
+ p.axi = axi ? axi->i : FALSE;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) diffusion_rhs, &p);
}
@@ -917,13 +936,9 @@ static void diffusion_relax (FttCell * cell, RelaxParams * p)
FttCellNeighbors neighbor;
FttCellFace face;
- if (GFS_IS_MIXED (cell)) {
- a = GFS_STATE (cell)->solid->a*GFS_VARIABLE (cell, p->dia);
- if (((cell)->flags & GFS_FLAG_DIRICHLET) != 0)
- g.b = gfs_cell_dirichlet_gradient_flux (cell, p->u, p->maxlevel, 0.);
- }
- else
- a = GFS_VARIABLE (cell, p->dia);
+ a = GFS_VARIABLE (cell, p->dia);
+ if (GFS_IS_MIXED (cell) && ((cell)->flags & GFS_FLAG_DIRICHLET) != 0)
+ g.b = gfs_cell_dirichlet_gradient_flux (cell, p->u, p->maxlevel, 0.);
face.cell = cell;
ftt_cell_neighbors (cell, &neighbor);
@@ -932,12 +947,18 @@ static void diffusion_relax (FttCell * cell, RelaxParams * p)
face.neighbor = neighbor.c[face.d];
gfs_face_gradient_flux (&face, &ng, p->u, p->maxlevel);
+ if (face.d/2 == p->component) {
+ ng.a *= 2.;
+ ng.b *= 2.;
+ }
g.a += ng.a;
g.b += ng.b;
}
a *= h*h;
g_assert (a > 0.);
g.a = 1. + g.a/a;
+ if (p->axi)
+ g.a += GFS_VARIABLE (cell, p->axi);
g.b = GFS_VARIABLE (cell, p->res) + g.b/a;
g_assert (g.a > 0.);
GFS_VARIABLE (cell, p->u) = g.b/g.a;
@@ -952,15 +973,13 @@ static void diffusion_residual (FttCell * cell, RelaxParams * p)
FttCellFace face;
h = ftt_cell_size (cell);
+ a = GFS_VARIABLE (cell, p->dia);
if (GFS_IS_MIXED (cell)) {
- a = GFS_STATE (cell)->solid->a*GFS_VARIABLE (cell, p->dia);
if (((cell)->flags & GFS_FLAG_DIRICHLET) != 0)
g.b = gfs_cell_dirichlet_gradient_flux (cell, p->u, -1, GFS_STATE (cell)->solid->fv);
else
g.b = GFS_STATE (cell)->solid->fv;
}
- else
- a = GFS_VARIABLE (cell, p->dia);
face.cell = cell;
ftt_cell_neighbors (cell, &neighbor);
@@ -969,12 +988,18 @@ static void diffusion_residual (FttCell * cell, RelaxParams * p)
face.neighbor = neighbor.c[face.d];
gfs_face_gradient_flux (&face, &ng, p->u, -1);
+ if (face.d/2 == p->component) {
+ ng.a *= 2.;
+ ng.b *= 2.;
+ }
g.a += ng.a;
g.b += ng.b;
}
a *= h*h;
g_assert (a > 0.);
g.a = 1. + g.a/a;
+ if (p->axi)
+ g.a += GFS_VARIABLE (cell, p->axi);
g.b = GFS_VARIABLE (cell, p->rhs) + g.b/a;
GFS_VARIABLE (cell, p->res) = g.b - g.a*GFS_VARIABLE (cell, p->u);
}
@@ -984,7 +1009,8 @@ static void diffusion_residual (FttCell * cell, RelaxParams * p)
* @domain: a #GfsDomain.
* @u: the variable to use as left-hand side.
* @rhs: the right-hand side.
- * @dia: the diagonal weight.
+ * @rhoc: the mass.
+ * @axi: the axisymmetric term.
* @res: the residual.
*
* Sets the @res variable of each leaf cell of @domain to the residual
@@ -997,7 +1023,8 @@ static void diffusion_residual (FttCell * cell, RelaxParams * p)
void gfs_diffusion_residual (GfsDomain * domain,
GfsVariable * u,
GfsVariable * rhs,
- GfsVariable * dia,
+ GfsVariable * rhoc,
+ GfsVariable * axi,
GfsVariable * res)
{
RelaxParams p;
@@ -1005,13 +1032,15 @@ void gfs_diffusion_residual (GfsDomain * domain,
g_return_if_fail (domain != NULL);
g_return_if_fail (u != NULL);
g_return_if_fail (rhs != NULL);
- g_return_if_fail (dia != NULL);
+ g_return_if_fail (rhoc != NULL);
g_return_if_fail (res != NULL);
p.u = u->i;
p.rhs = rhs->i;
- p.dia = dia->i;
+ p.dia = rhoc->i;
p.res = res->i;
+ p.component = GFS_IS_AXI (domain) ? u->component : FTT_DIMENSION;
+ p.axi = axi ? axi->i : FALSE;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) diffusion_residual, &p);
}
@@ -1024,7 +1053,8 @@ void gfs_diffusion_residual (GfsDomain * domain,
* @nrelax: the number of relaxations to apply at each level.
* @u: the variable to use as left-hand side.
* @rhs: the right-hand side.
- * @dia: the diagonal weight.
+ * @rhoc: the mass.
+ * @axi: the axisymmetric term.
* @res: the residual.
*
* Apply one multigrid iteration to the diffusion equation for @u.
@@ -1043,7 +1073,8 @@ void gfs_diffusion_cycle (GfsDomain * domain,
guint nrelax,
GfsVariable * u,
GfsVariable * rhs,
- GfsVariable * dia,
+ GfsVariable * rhoc,
+ GfsVariable * axi,
GfsVariable * res)
{
guint n;
@@ -1054,11 +1085,10 @@ void gfs_diffusion_cycle (GfsDomain * domain,
g_return_if_fail (domain != NULL);
g_return_if_fail (u != NULL);
g_return_if_fail (rhs != NULL);
- g_return_if_fail (dia != NULL);
+ g_return_if_fail (rhoc != NULL);
g_return_if_fail (res != NULL);
dp = gfs_temporary_variable (domain);
- dp->component = u->component;
/* compute residual on non-leafs cells */
gfs_domain_cell_traverse (domain,
@@ -1072,7 +1102,9 @@ void gfs_diffusion_cycle (GfsDomain * domain,
p.maxlevel = levelmin;
p.u = dp->i;
p.res = res->i;
- p.dia = dia->i;
+ p.dia = rhoc->i;
+ p.component = GFS_IS_AXI (domain) ? u->component : FTT_DIMENSION;
+ p.axi = axi ? axi->i : FALSE;
for (n = 0; n < 10*nrelax; n++) {
gfs_domain_homogeneous_bc (domain,
FTT_TRAVERSE_LEVEL | FTT_TRAVERSE_LEAFS,
@@ -1106,7 +1138,7 @@ void gfs_diffusion_cycle (GfsDomain * domain,
(FttCellTraverseFunc) correct, data);
gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, u);
/* compute new residual on leaf cells */
- gfs_diffusion_residual (domain, u, rhs, dia, res);
+ gfs_diffusion_residual (domain, u, rhs, rhoc, axi, res);
gts_object_destroy (GTS_OBJECT (dp));
}
diff --git a/src/poisson.h b/src/poisson.h
index f26e5d9..900782b 100644
--- a/src/poisson.h
+++ b/src/poisson.h
@@ -77,18 +77,21 @@ void gfs_poisson_cycle (GfsDomain * domain,
void gfs_diffusion_coefficients (GfsDomain * domain,
GfsSourceDiffusion * d,
gdouble dt,
- GfsVariable * dia,
+ GfsVariable * rhoc,
+ GfsVariable * axi,
GfsFunction * alpha,
gdouble beta);
void gfs_diffusion_rhs (GfsDomain * domain,
GfsVariable * v,
GfsVariable * rhs,
- GfsVariable * dia,
+ GfsVariable * rhoc,
+ GfsVariable * axi,
gdouble beta);
void gfs_diffusion_residual (GfsDomain * domain,
GfsVariable * u,
GfsVariable * rhs,
- GfsVariable * dia,
+ GfsVariable * rhoc,
+ GfsVariable * axi,
GfsVariable * res);
void gfs_diffusion_cycle (GfsDomain * domain,
guint levelmin,
@@ -96,7 +99,8 @@ void gfs_diffusion_cycle (GfsDomain * domain,
guint nrelax,
GfsVariable * u,
GfsVariable * rhs,
- GfsVariable * dia,
+ GfsVariable * rhoc,
+ GfsVariable * axi,
GfsVariable * res);
#ifdef __cplusplus
diff --git a/src/source.c b/src/source.c
index 94b02e7..8fa4941 100644
--- a/src/source.c
+++ b/src/source.c
@@ -119,6 +119,9 @@ GfsVariable * gfs_domain_variable_fluxes (GfsDomain * domain,
g_return_val_if_fail (domain != NULL, NULL);
g_return_val_if_fail (v != NULL, NULL);
+ if (!v->sources)
+ return NULL;
+
GSList * i = GTS_SLIST_CONTAINER (v->sources)->items;
while (i) {
if (GFS_SOURCE_GENERIC (i->data)->flux) {
@@ -892,6 +895,65 @@ GfsSourceGenericClass * gfs_source_diffusion_explicit_class (void)
/* GfsSourceViscosity: Object */
+typedef struct {
+ GfsSourceGeneric * s;
+ GfsVariable * v, * sv, * tv;
+ gdouble dt;
+} FluxPar;
+
+static void add_viscosity_transverse_flux (FttCell * cell, FluxPar * p)
+{
+ FttCellFace f;
+ FttCellNeighbors n;
+ gdouble transverse = 0.;
+
+ f.cell = cell;
+ ftt_cell_neighbors (cell, &n);
+#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];
+ transverse += (FTT_FACE_DIRECT (&f) ? 1. : -1.)*
+ gfs_face_weighted_interpolated_value (&f, p->tv->i);
+ }
+#else
+ g_assert_not_implemented ();
+#endif
+
+ GfsFunction * alpha = gfs_object_simulation (p->s)->physical_params.alpha;
+ gdouble h = ftt_cell_size (cell);
+ GFS_VALUE (cell, p->sv) += (alpha ? gfs_function_value (alpha, cell) : 1.)*transverse/h;
+}
+
+static void compute_transverse (FttCell * cell, FluxPar * p)
+{
+ GfsVariable ** v = GFS_SOURCE_VISCOSITY (p->s)->v;
+ gdouble h = ftt_cell_size (cell);
+ GFS_VALUE (cell, p->tv) = gfs_center_gradient (cell,
+ p->v->component,
+ v[(p->v->component + 1) % FTT_DIMENSION]->i)/h;
+}
+
+static void source_viscosity_transverse_flux (GfsSourceGeneric * s,
+ GfsDomain * domain,
+ GfsVariable * v, GfsVariable * sv,
+ gdouble dt)
+{
+ FluxPar p;
+
+ gfs_diffusion_coefficients (domain, GFS_SOURCE_DIFFUSION (s), dt, NULL, NULL, NULL, 1.);
+ p.s = s;
+ p.v = v;
+ p.sv = sv;
+ p.dt = dt;
+ p.tv = gfs_temporary_variable (domain);
+ gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) compute_transverse, &p);
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, p.tv);
+ gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) add_viscosity_transverse_flux, &p);
+ gts_object_destroy (GTS_OBJECT (p.tv));
+}
+
static void source_viscosity_read (GtsObject ** o, GtsFile * fp)
{
GfsSourceViscosity * source;
@@ -906,10 +968,10 @@ 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;
+ GfsSourceGeneric * s = GFS_SOURCE_GENERIC (source);
+ s->mac_value = NULL;
+ s->centered_value = NULL;
+ s->flux = source_viscosity_transverse_flux;
}
if (!(source->v = gfs_domain_velocity (domain))) {
gts_file_error (fp, "cannot find velocity components");
@@ -975,6 +1037,7 @@ static void source_viscosity_init (GfsSourceGeneric * s)
{
s->mac_value = source_viscosity_value;
s->centered_value = source_viscosity_non_diffusion_value;
+ s->flux = NULL;
}
GfsSourceGenericClass * gfs_source_viscosity_class (void)
@@ -1000,31 +1063,6 @@ 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)
{
@@ -1038,44 +1076,11 @@ static gdouble source_viscosity_stability (GfsSourceGeneric * s,
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. */
@@ -1086,7 +1091,7 @@ static void add_viscosity_explicit_flux (FttCell * cell, FluxPar * p)
FttCellFace f;
FttCellNeighbors n;
GfsGradient g = { 0., 0. };
- gdouble v0, h;
+ gdouble v0;
if (GFS_IS_MIXED (cell)) {
if (((cell)->flags & GFS_FLAG_DIRICHLET) != 0)
@@ -1110,7 +1115,6 @@ static void add_viscosity_explicit_flux (FttCell * cell, FluxPar * p)
g.a += e.a;
g.b += e.b;
}
- h = ftt_cell_size (cell);
gdouble transverse = 0.;
#ifndef NOTRANSVERSE
@@ -1119,19 +1123,8 @@ static void add_viscosity_explicit_flux (FttCell * cell, FluxPar * p)
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
+ gfs_face_weighted_interpolated_value (&f, p->tv->i);
}
#else
g_assert_not_implemented ();
@@ -1139,6 +1132,7 @@ static void add_viscosity_explicit_flux (FttCell * cell, FluxPar * p)
#endif
GfsFunction * alpha = gfs_object_simulation (p->s)->physical_params.alpha;
+ gdouble h = ftt_cell_size (cell);
GFS_VALUE (cell, p->sv) += (alpha ? gfs_function_value (alpha, cell) : 1.)*
((g.b - g.a*v0)/h + transverse)/h;
}
@@ -1162,17 +1156,19 @@ static void source_viscosity_explicit_flux (GfsSourceGeneric * s,
{
FluxPar p;
- gfs_diffusion_coefficients (domain, GFS_SOURCE_DIFFUSION (s), dt, NULL, NULL, 1.);
+ gfs_diffusion_coefficients (domain, GFS_SOURCE_DIFFUSION (s), dt, NULL, NULL, NULL, 1.);
gfs_domain_surface_bc (domain, v);
p.s = s;
p.v = v;
p.sv = sv;
p.dt = dt;
+ p.tv = gfs_temporary_variable (domain);
+ gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) compute_transverse, &p);
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, p.tv);
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
+ gts_object_destroy (GTS_OBJECT (p.tv));
}
static void source_viscosity_explicit_init (GfsSourceGeneric * s)
@@ -1189,7 +1185,7 @@ GfsSourceGenericClass * gfs_source_viscosity_explicit_class (void)
if (klass == NULL) {
GtsObjectClassInfo source_viscosity_explicit_info = {
"GfsSourceViscosityExplicit",
- sizeof (GfsSourceViscosityExplicit),
+ sizeof (GfsSourceViscosity),
sizeof (GfsSourceGenericClass),
(GtsObjectClassInitFunc) source_viscosity_explicit_class_init,
(GtsObjectInitFunc) source_viscosity_explicit_init,
diff --git a/src/source.h b/src/source.h
index b373a49..660b0f7 100644
--- a/src/source.h
+++ b/src/source.h
@@ -262,19 +262,6 @@ 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 ac17ce5..001c90f 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -580,7 +580,8 @@ void gfs_predicted_face_velocities (GfsDomain * domain,
* @par: the multilevel parameters.
* @v: a #GfsVariable.
* @rhs: the right-hand side.
- * @dia: the diagonal weight.
+ * @rhoc: the mass.
+ * @axi: the axisymmetric term.
*
* Solves a diffusion equation for variable @v using a Crank-Nicholson
* scheme with multilevel relaxations.
@@ -593,7 +594,8 @@ void gfs_diffusion (GfsDomain * domain,
GfsMultilevelParams * par,
GfsVariable * v,
GfsVariable * rhs,
- GfsVariable * dia)
+ GfsVariable * rhoc,
+ GfsVariable * axi)
{
guint minlevel, maxlevel;
GfsVariable * res;
@@ -602,7 +604,7 @@ void gfs_diffusion (GfsDomain * domain,
g_return_if_fail (par != NULL);
g_return_if_fail (v != NULL);
g_return_if_fail (rhs != NULL);
- g_return_if_fail (dia != NULL);
+ g_return_if_fail (rhoc != NULL);
res = gfs_temporary_variable (domain);
@@ -610,14 +612,14 @@ void gfs_diffusion (GfsDomain * domain,
if (par->minlevel > minlevel)
minlevel = par->minlevel;
maxlevel = gfs_domain_depth (domain);
- gfs_diffusion_residual (domain, v, rhs, dia, res);
+ gfs_diffusion_residual (domain, v, rhs, rhoc, axi, res);
par->residual_before = par->residual =
gfs_domain_norm_variable (domain, res, NULL, FTT_TRAVERSE_LEAFS, -1);
gdouble res_max_before = par->residual.infty;
par->niter = 0;
while (par->niter < par->nitermin ||
(par->residual.infty > par->tolerance && par->niter < par->nitermax)) {
- gfs_diffusion_cycle (domain, minlevel, maxlevel, par->nrelax, v, rhs, dia, res);
+ gfs_diffusion_cycle (domain, minlevel, maxlevel, par->nrelax, v, rhs, rhoc, axi, res);
par->residual = gfs_domain_norm_variable (domain, res, NULL, FTT_TRAVERSE_LEAFS, -1);
if (par->residual.infty == res_max_before) /* convergence has stopped!! */
break;
@@ -719,19 +721,22 @@ static void variable_diffusion (GfsDomain * domain,
GfsVariable * rhs,
GfsFunction * alpha)
{
- GfsVariable * dia;
+ GfsVariable * rhoc, * axi;
- dia = gfs_temporary_variable (domain);
+ rhoc = gfs_temporary_variable (domain);
+ axi = GFS_IS_AXI (domain) && par->v->component == FTT_Y ? gfs_temporary_variable (domain) : NULL;
- gfs_diffusion_coefficients (domain, d, par->dt, dia, alpha, d->D->par.beta);
+ gfs_diffusion_coefficients (domain, d, par->dt, rhoc, axi, alpha, d->D->par.beta);
gfs_domain_surface_bc (domain, par->v);
- gfs_diffusion_rhs (domain, par->v, rhs, dia, d->D->par.beta);
+ gfs_diffusion_rhs (domain, par->v, rhs, rhoc, axi, d->D->par.beta);
/* fixme: time shoud be set to t + dt here in case boundary values are
time-dependent in the call below */
gfs_domain_surface_bc (domain, par->v);
- gfs_diffusion (domain, &d->D->par, par->v, rhs, dia);
+ gfs_diffusion (domain, &d->D->par, par->v, rhs, rhoc, axi);
- gts_object_destroy (GTS_OBJECT (dia));
+ if (axi)
+ gts_object_destroy (GTS_OBJECT (axi));
+ gts_object_destroy (GTS_OBJECT (rhoc));
}
/**
diff --git a/src/timestep.h b/src/timestep.h
index db11a99..a10e118 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -65,7 +65,8 @@ void gfs_diffusion (GfsDomain * domain,
GfsMultilevelParams * par,
GfsVariable * v,
GfsVariable * rhs,
- GfsVariable * dia);
+ GfsVariable * rhoc,
+ GfsVariable * axi);
void gfs_centered_velocity_advection_diffusion (GfsDomain * domain,
guint dimension,
GfsAdvectionParams * apar,
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list