[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:35 UTC 2009
The following commit has been merged in the upstream branch:
commit 149225e729d8b1588c58426956e12d037773cb3e
Author: Stephane Popinet <popinet at users.sf.net>
Date: Sat Jul 19 11:31:26 2008 +1000
Axisymmetric Euler solver
darcs-hash:20080719013126-d4795-81325b1553a6ae0017f2276f729b7e0d42cf0ce3.gz
diff --git a/doc/figures/axi.tm b/doc/figures/axi.tm
new file mode 100644
index 0000000..aebc376
--- /dev/null
+++ b/doc/figures/axi.tm
@@ -0,0 +1,107 @@
+<TeXmacs|1.0.6.10>
+
+<style|article>
+
+<\body>
+ The incompressible Navier--Stokes equations in cylindrical coordinates are
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|<frac|\<partial\>v<rsub|r>|\<partial\>t>+<frac|1|r>*<frac|\<partial\>(r*v<rsub|r><rsup|2>)|\<partial\>r>+<frac|\<partial\>(v<rsub|r*>*v<rsub|z>)|\<partial\>z>=-<frac|\<partial\>\<phi\>|\<partial\>r>+<frac|1|r>*<frac|\<partial\>(r*S<rsub|rr>)|\<partial\>r>+<frac|\<partial\>S<rsub|zr>|\<partial\>z>-<frac|S<rsub|\<theta\>\<theta\>>|r>,>|<cell|<eq-number><label|momr>>>|<row|<cell|>|<cell|<frac|\<partial\>v<rsub|z>|\<partial\>t>+<frac|1|r>*<frac|\<partial\>(r*v<rsub|r>*v<rsub|z>)|\<partial\>r>+<frac|\<partial\>(v<rsup|2><rsub|z>)|\<partial\>z>=-<frac|\<partial\>\<phi\>|\<partial\>z>+<frac|1|r>*<frac|\<partial\>(r*S<rsub|zr>)|\<partial\>r>+<frac|\<partial\>S<rsub|zz>|\<partial\>z>,>|<cell|<eq-number><label|momz>>>|<row|<cell|>|<cell|<frac|1|r>*<frac|\<partial\>(r*v<rsub|r>)|\<partial\>r>+<frac|\<partial\>v<rsub|z>|\<partial\>z>=0,>|<cell|<eq-number><label|continuity>>>>>
+ </eqnarray*>
+
+ with <math|\<phi\>=p/\<rho\>> and the stress tensor
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|S<rsub|rr>>|<cell|=>|<cell|2*\<nu\>*<frac|\<partial\>v<rsub|r>|\<partial\>r>,>>|<row|<cell|S<rsub|\<theta\>\<theta\>>>|<cell|=>|<cell|2*\<nu\>*<frac|v<rsub|r>|r>,>>|<row|<cell|S<rsub|zz>>|<cell|=>|<cell|2*\<nu\>*<frac|\<partial\>v<rsub|z>|\<partial\>z>,>>|<row|<cell|S<rsub|zr>>|<cell|=>|<cell|\<nu\>*<left|(><frac|\<partial\>v<rsub|z>|\<partial\>r>+<frac|\<partial\>v<rsub|r>|\<partial\>z><right|)>.>>>>
+ </eqnarray*>
+
+ Considering a control volume <math|\<Omega\>> with boundary
+ <math|\<partial\>\<Omega\>>, the integral equations can then be written
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|\<partial\><rsub|t><big|int><rsub|\<Omega\>>v<rsub|r>*r*dr*dz+<big|int><rsub|\<partial\>\<Omega\>>r*v<rsub|r>*v<rsub|z>*dr+<big|int><rsub|\<partial\>\<Omega\>>r*v<rsub|r>*v<rsub|r>*dz=>|<cell|>>|<row|<cell|>|<cell|-<big|int><rsub|\<partial\>\<Omega\>>\<phi\>*r*dz+<big|int><rsub|\<Omega\>>\<phi\>*dr*dz+<big|int><rsub|\<partial\>\<Omega\>>r*S<rsub|zr>*dr+<big|int><rsub|\<partial\>\<Omega\>>r*S<rsub|rr>*dz-<big|int><rsub|\<Omega\>>S<rsub|\<theta\>\<theta\>>*dr*dz,>|<cell|>>|<row|<cell|>|<cell|\<partial\><rsub|t><big|int><rsub|\<Omega\>>v<rsub|z>*r*dr*dz+<big|int><rsub|\<partial\>\<Omega\>>r*v<rsub|z>*v<rsub|z>*dr+<big|int><rsub|\<partial\>\<Omega\>>r*v<rsub|z>*v<rsub|r>*dz=>|<cell|>>|<row|<cell|>|<cell|-<big|int><rsub|\<partial\>\<Omega\>>\<phi\>*r*dr+<big|int><rsub|\<partial\>\<Omega\>>r*S<rsub|zz>*dr+<big|int><rsub|\<partial\>\<Omega\>>r*S<rsub|zr>*dz,>|<cell|>>|<row|<cell|>|<cell|<big|int><rsub|\<partial\>\<Omega\>>r*v<rsub|z>*dr+<big|int><rsub|\<partial\>\<Omega\>>r*v<rsub|r>*dz=0.>|<cell|>>|<row|<cell|>|<cell|>|<cell|>>>>
+ </eqnarray*>
+
+ Posing
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|<with|mode|text|<math|a<rsup|i,j>>>>|<cell|\<equiv\>>|<cell|<big|int><rsub|\<Omega\>>dr*dz,>>|<row|<cell|c<rsup|i,j>>|<cell|\<equiv\>>|<cell|<big|int><rsub|\<Omega\>>r*dr*dz\<approx\>r<rsup|j>*a<rsup|i,j>,>>|<row|<cell|s<rsub|z><rsup|i,j>>|<cell|\<equiv\>>|<cell|<big|int><rsub|\<partial\>\<Omega\><rsup|j>>r*dr\<approx\>r<rsup|j>*<big|int><rsub|\<partial\>\<Omega\><rsup|j>>dr,>>|<row|<cell|s<rsup|i,j-1/2><rsub|r>>|<cell|\<equiv\>>|<cell|<big|int><rsub|\<partial\>\<Omega\><rsup|j-1/2>>r*dz<rsub|>=r<rsup|j-1/2>*<big|int><rsub|\<partial\>\<Omega\><rsup|j-1/2>>dz,>>>>
+ </eqnarray*>
+
+ we get the discrete form of the equations as
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|\<partial\><rsub|t>(c*v<rsub|r>)<rsup|i,j>+(s<rsub|z>*v<rsub|r>*v<rsub|z>)<rsup|i+1/2,j>-(s<rsub|z>*v<rsub|r>*v<rsub|z>)<rsup|i-1/2,j>+(s<rsub|r>*v<rsub|r><rsup|2>)<rsup|i,j+1/2>-(s<rsub|r>*v<rsub|r><rsup|2>)<rsup|i,j-1/2>=>|<cell|>>|<row|<cell|>|<cell|(s<rsub|r>\<phi\>)<rsup|i,j-1/2>-(s<rsub|r>\<phi\>)<rsup|i,j+1/2>+<with|color|blue|<with|math-font-series|bold|(a*\<phi\>)<rsup|i,j>>>+>|<cell|>>|<row|<cell|>|<cell|(s<rsub|z>*S<rsub|zr>)<rsup|i+1/2,j>-(s<rsub|z>*S<rsub|zr>)<rsup|i-1/2,j>+(s<rsub|r>*S<rsub|rr>)<rsup|i,j+1/2>-(s<rsub|r>*S<rsub|rr>)<rsup|i,j-1/2>-<with|color|blue|<with|math-font-series|bold|(a*S<rsub|\<theta\>\<theta\>>)<rsup|i,j>>>,>|<cell|>>>>
+ </eqnarray*>
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|\<partial\><rsub|t>(c*v<rsub|z>)<rsup|i,j>+(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>\<phi\>)<rsup|i-1/2,j>-(s<rsub|z>\<phi\>)<rsup|i+1/2,j>+>|<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*>
+
+ <\equation*>
+ (s<rsub|z>*v<rsub|z>)<rsup|i+1/2,j>-(s<rsub|z>*v<rsub|z>)<rsup|i-1/2,j>+(s<rsub|r>*v<rsub|r>)<rsup|i,j+1/2>-(s<rsub|r>*v<rsub|r>)<rsup|i,j-1/2>=0,
+ </equation*>
+
+ where only the bold terms differ from the discretisation of the N--S
+ equations in Cartesian coordinates in \ two dimensions.
+
+ The projection method relies on a decomposition of the velocity as
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|v<rsup|i,j><rsub|r>>|<cell|=>|<cell|(v<rsub|r><rsup|\<star\>>)<rsup|i,j>+<frac|\<Delta\>t|c<rsup|i,j>><left|[>(s<rsub|r>\<phi\>)<rsup|i,j-1/2>-(s<rsub|r>\<phi\>)<rsup|i,j+1/2>+<with|color|blue|<with|math-font-series|bold|(a*\<phi\>)<rsup|i,j>>><right|]>,>>|<row|<cell|v<rsup|i,j><rsub|z>>|<cell|=>|<cell|(v<rsub|z><rsup|\<star\>>)<rsup|i,j>+<frac|\<Delta\>t|c<rsup|i,j>><left|[>(s<rsub|z>\<phi\>)<rsup|i-1/2,j>-(s<rsub|z>\<phi\>)<rsup|i+1/2,j><right|]>,>>>>
+ </eqnarray*>
+
+ wich leads to the Poisson-like equation
+
+ <\equation*>
+ <frac|1|\<Delta\>t>*<left|[>s<rsup|i+1/2,j><rsub|z>**(v<rsub|z><rsup|\<star\>>)<rsup|i+1/2,j>-s<rsup|i-1/2,j><rsub|z>*(v<rsub|z><rsup|\<star\>>)<rsup|i-1/2,j>+s<rsup|i,j+1/2><rsub|r>*(v<rsub|r><rsup|\<star\>>)<rsup|i,j+1/2>-s<rsup|i,j-1/2><rsub|r>*(v<rsub|r><rsup|\<star\>>)<rsup|i,j-1/2><right|]>+\<phi\><rsup|i,j>*<left|(>s<rsup|i,j><rsub|z>*<left|[>(s<rsub|z>/c)<rsup|i+1/2,j>+(s<rsub|z>/c)<rsup|i-1/2,j><right|]>+s<rsup|i,j><rsub|r>*<left|[>(s<rsub|r>/c)<rsup|i,j+1/2>+(s<rsub|r>/c)<rsup|i,j-1/2><right|]><right|)>+s<rsup|i,j+1/2><rsub|r>*(a*\<phi\>)<rsup|i,j+1/2>-s<rsup|i,j-1/2><rsub|r>*(a*\<phi\>)<rsup|i,j-1/2>-(s<rsub|z>/c)<rsup|i+1/2,j>*(s<rsub|z>\<phi\>)<rsup|i+1,j>-(s<rsub|z>/c)<rsup|i-1/2,j>*(s<rsub|z>\<phi\>)<rsup|i-1,j>-(s<rsub|r>/c)<rsup|i,j+1/2>*(s<rsub|r>\<phi\>)<rsup|i,j+1>-(s<rsub|r>/c)<rsup|i,j-1/2>*(s<rsub|r>\<phi\>)<rsup|i,j-1>=0.
+ </equation*>
+
+ This looks complicated though... What about just splitting the velocity as
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|v<rsup|i,j><rsub|r>>|<cell|=>|<cell|(v<rsub|r><rsup|\<star\>>)<rsup|i,j>+\<Delta\>t*<left|(>\<phi\><rsup|i,j-1/2>-\<phi\><rsup|i,j+1/2><right|)>,>>|<row|<cell|v<rsup|i,j><rsub|z>>|<cell|=>|<cell|(v<rsub|z><rsup|\<star\>>)<rsup|i,j>+\<Delta\>t*<left|(>\<phi\><rsup|i-1/2,j>-\<phi\><rsup|i+1/2,j><right|)>,>>>>
+ </eqnarray*>
+
+ which is still consistent with equations (<reference|momr>) and
+ (<reference|momz>). This gives
+
+ <\equation*>
+ <frac|1|\<Delta\>t>*<left|[>s<rsup|i+1/2,j><rsub|z>**(v<rsub|z><rsup|\<star\>>)<rsup|i+1/2,j>-s<rsup|i-1/2,j><rsub|z>*(v<rsub|z><rsup|\<star\>>)<rsup|i-1/2,j>+s<rsup|i,j+1/2><rsub|r>*(v<rsub|r><rsup|\<star\>>)<rsup|i,j+1/2>-s<rsup|i,j-1/2><rsub|r>*(v<rsub|r><rsup|\<star\>>)<rsup|i,j-1/2><right|]>+\<phi\><rsup|i,j>*(s<rsup|i+1/2,j><rsub|z>+s<rsup|i-1/2,j><rsub|z>+s<rsup|i,j+1/2><rsub|r>+s<rsup|i,j-1/2><rsub|r>)-s<rsub|z><rsup|i+1/2,j>*\<phi\><rsup|i+1,j>-s<rsub|z><rsup|i-1/2,j>*\<phi\><rsup|i-1,j>-s<rsub|r><rsup|i,j+1/2>*\<phi\><rsup|i,j+1>-s<rsub|r><rsup|i,j-1/2>*\<phi\><rsup|i,j-1>=0,
+ </equation*>
+
+ together with
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|c<rsup|i,j>*<frac|(v<rsub|r><rsup|\<star\>>)<rsup|i,j>-v<rsup|i,j><rsub|r>|\<Delta\>t>+(s<rsub|z>*v<rsub|r>*v<rsub|z>)<rsup|i+1/2,j>-(s<rsub|z>*v<rsub|r>*v<rsub|z>)<rsup|i-1/2,j>+(s<rsub|r>*v<rsub|r><rsup|2>)<rsup|i,j+1/2>-(s<rsub|r>*v<rsub|r><rsup|2>)<rsup|i,j-1/2>=>|<cell|>>|<row|<cell|>|<cell|(s<rsub|z>*S<rsub|zr>)<rsup|i+1/2,j>-(s<rsub|z>*S<rsub|zr>)<rsup|i-1/2,j>+(s<rsub|r>*S<rsub|rr>)<rsup|i,j+1/2>-(s<rsub|r>*S<rsub|rr>)<rsup|i,j-1/2>-<with|color|blue|<with|math-font-series|bold|(a*S<rsub|\<theta\>\<theta\>>)<rsup|i,j>>>,>|<cell|>>>>
+ </eqnarray*>
+
+ <\eqnarray*>
+ <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>
+
+ <\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>),
+ </equation*>
+
+ using equations (<reference|momr>), (<reference|momz>) and
+ (<reference|continuity>) then gives
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|v<rsup|n+1/2><rsub|r>=v<rsub|r>+<frac|h|2>*\<partial\><rsub|r>v<rsub|r><rsup|n>+<frac|\<Delta\>t|2>*<left|(>-<frac|1|r>*<frac|\<partial\>(r*v<rsub|r><rsup|2>)|\<partial\>r>-<frac|\<partial\>(v<rsub|r*>*v<rsub|z>)|\<partial\>z>+src<rsup|n><left|)>>|<cell|>|<cell|>>>>
+ </eqnarray*>
+
+ \;
+</body>
+
+<\references>
+ <\collection>
+ <associate|auto-1|<tuple|1|?>>
+ <associate|continuity|<tuple|3|?>>
+ <associate|momr|<tuple|1|?>>
+ <associate|momz|<tuple|2|?>>
+ <associate|mon|<tuple|?|?>>
+ </collection>
+</references>
\ No newline at end of file
diff --git a/src/advection.c b/src/advection.c
index 0266330..ff07584 100644
--- a/src/advection.c
+++ b/src/advection.c
@@ -392,7 +392,7 @@ void gfs_face_advection_flux (const FttCellFace * face,
g_return_if_fail (face != NULL);
g_return_if_fail (par != NULL);
- flux = GFS_FACE_FRACTION (face)*GFS_FACE_NORMAL_VELOCITY (face)*par->dt*
+ flux = gfs_domain_face_fraction (par->v->domain, face)*GFS_FACE_NORMAL_VELOCITY (face)*par->dt*
gfs_face_upwinded_value (face, GFS_FACE_UPWINDING, NULL)/ftt_cell_size (face->cell);
if (!FTT_FACE_DIRECT (face))
flux = - flux;
@@ -437,7 +437,7 @@ void gfs_face_velocity_advection_flux (const FttCellFace * face,
c = par->v->component;
g_return_if_fail (c >= 0 && c < FTT_DIMENSION);
- flux = GFS_FACE_FRACTION (face)*GFS_FACE_NORMAL_VELOCITY (face)*par->dt
+ flux = gfs_domain_face_fraction (par->v->domain, face)*GFS_FACE_NORMAL_VELOCITY (face)*par->dt
/ftt_cell_size (face->cell);
#if 0
if (c == face->d/2) /* normal component */
@@ -487,7 +487,7 @@ void gfs_face_velocity_convective_flux (const FttCellFace * face,
g_return_if_fail (face != NULL);
g_return_if_fail (par != NULL);
- g_return_if_fail (GFS_FACE_FRACTION (face) == 1.);
+ g_return_if_fail (gfs_domain_face_fraction (par->v->domain, face) == 1.);
c = par->v->component;
g_return_if_fail (c >= 0 && c < FTT_DIMENSION);
@@ -560,8 +560,8 @@ void gfs_face_advected_normal_velocity (const FttCellFace * face,
break;
case FTT_FINE_COARSE:
GFS_FACE_NORMAL_VELOCITY_RIGHT (face) +=
- u*GFS_FACE_FRACTION_LEFT (face)/(GFS_FACE_FRACTION_RIGHT (face)*
- FTT_CELLS_DIRECTION (face->d));
+ u*gfs_domain_face_fraction (par->v->domain, face)/
+ (gfs_domain_face_fraction_right (par->v->domain, face)*FTT_CELLS_DIRECTION (face->d));
break;
default:
g_assert_not_reached ();
@@ -594,8 +594,8 @@ void gfs_face_interpolated_normal_velocity (const FttCellFace * face, GfsVariabl
break;
case FTT_FINE_COARSE:
GFS_FACE_NORMAL_VELOCITY_RIGHT (face) +=
- u*GFS_FACE_FRACTION_LEFT (face)/(GFS_FACE_FRACTION_RIGHT (face)*
- FTT_CELLS_DIRECTION (face->d));
+ u*gfs_domain_face_fraction (v[0]->domain, face)/
+ (gfs_domain_face_fraction_right (v[0]->domain, face)*FTT_CELLS_DIRECTION (face->d));
break;
default:
g_assert_not_reached ();
diff --git a/src/domain.h b/src/domain.h
index b0862d5..72f5fe9 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -72,7 +72,8 @@ struct _GfsDomain {
struct _GfsDomainClass {
GtsWGraphClass parent_class;
- void (* post_read) (GfsDomain *, GtsFile * fp);
+ void (* post_read) (GfsDomain *, GtsFile * fp);
+ gdouble (* face_fraction) (const GfsDomain *, const FttCellFace *);
};
#define GFS_DOMAIN(obj) GTS_OBJECT_CAST (obj,\
@@ -291,6 +292,43 @@ void gfs_domain_sum (GfsDomain * domain,
void gfs_domain_filter (GfsDomain * domain,
GfsVariable * v,
GfsVariable * fv);
+/**
+ * gfs_domain_face_fraction:
+ * @domain; a #GfsDomain.
+ * @face: a #FttCellFace.
+ *
+ * Returns: the surface fraction of @face taking into account any
+ * orthogonal mapping of @domain.
+ */
+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);
+ return f;
+}
+
+/**
+ * gfs_domain_face_fraction_right:
+ * @domain; a #GfsDomain.
+ * @face: a #FttCellFace.
+ *
+ * Returns: the surface fraction "to the right" of @face taking into account any
+ * orthogonal mapping of @domain.
+ */
+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) {
+ 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);
+ }
+ return f;
+}
#ifdef __cplusplus
}
diff --git a/src/fluid.c b/src/fluid.c
index cc18d8f..fcccbfd 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -1818,29 +1818,16 @@ gdouble gfs_face_interpolated_value (const FttCellFace * face,
void gfs_normal_divergence (FttCell * cell,
GfsVariable * v)
{
- FttComponent c;
+ FttCellFace face;
gdouble div = 0.;
g_return_if_fail (cell != NULL);
g_return_if_fail (v != NULL);
- if (GFS_IS_MIXED (cell)) {
- GfsSolidVector * solid = GFS_STATE (cell)->solid;
-
- for (c = 0; c < FTT_DIMENSION; c++) {
- FttDirection d = 2*c;
-
- div += (solid->s[d]*GFS_STATE (cell)->f[d].un -
- solid->s[d + 1]*GFS_STATE (cell)->f[d + 1].un);
- }
- }
- else
- for (c = 0; c < FTT_DIMENSION; c++) {
- FttDirection d = 2*c;
-
- div += (GFS_STATE (cell)->f[d].un -
- GFS_STATE (cell)->f[d + 1].un);
- }
+ face.cell = cell;
+ for (face.d = 0; face.d < FTT_NEIGHBORS; face.d++)
+ div += (FTT_FACE_DIRECT (&face) ? 1. : -1.)*
+ GFS_STATE (cell)->f[face.d].un*gfs_domain_face_fraction (v->domain, &face);
GFS_VARIABLE (cell, v->i) = div*ftt_cell_size (cell);
}
diff --git a/src/init.c b/src/init.c
index 5e0ca28..c427119 100644
--- a/src/init.c
+++ b/src/init.c
@@ -100,6 +100,7 @@ GtsObjectClass ** gfs_classes (void)
gfs_ocean_class (),
gfs_advection_class (),
gfs_poisson_class (),
+ gfs_axi_class (),
gfs_wave_class (),
gfs_surface_bc_class (),
diff --git a/src/poisson.c b/src/poisson.c
index 18b30bc..38279d7 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -369,25 +369,17 @@ static void reset_coeff (FttCell * cell)
f[d].v = 0.;
}
-static void poisson_coeff (FttCell * cell, gdouble * lambda2)
-{
- GfsStateVector * s = GFS_STATE (cell);
- FttDirection d;
-
- for (d = 0; d < FTT_NEIGHBORS; d++)
- s->f[d].v = lambda2[d/2];
- if (s->solid)
- for (d = 0; d < FTT_NEIGHBORS; d++)
- s->f[d].v *= s->solid->s[d];
-}
+typedef struct {
+ gdouble lambda2[FTT_DIMENSION];
+ GfsFunction * alpha;
+ GfsDomain * domain;
+} PoissonCoeff;
-static void poisson_alpha_coeff (FttCellFace * face,
- gpointer * data)
+static void poisson_coeff (FttCellFace * face,
+ PoissonCoeff * p)
{
- gdouble * lambda2 = data[0];
- gdouble alpha = gfs_function_face_value (data[1], face);
- gdouble v = lambda2[face->d/2]*alpha;
- GfsStateVector * s = GFS_STATE (face->cell);
+ gdouble alpha = p->alpha ? gfs_function_face_value (p->alpha, face) : 1.;
+ gdouble v = p->lambda2[face->d/2]*alpha*gfs_domain_face_fraction (p->domain, face);
if (alpha <= 0.) {
FttVector p;
@@ -397,10 +389,8 @@ static void poisson_alpha_coeff (FttCellFace * face,
"Please check your definition.",
alpha, p.x, p.y, p.z);
}
- 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:
GFS_STATE (face->neighbor)->f[FTT_OPPOSITE_DIRECTION (face->d)].v = v;
@@ -454,7 +444,7 @@ static void face_coeff_from_below (FttCell * cell)
void gfs_poisson_coefficients (GfsDomain * domain,
GfsFunction * alpha)
{
- gdouble lambda2[FTT_DIMENSION];
+ PoissonCoeff p;
FttComponent i;
g_return_if_fail (domain != NULL);
@@ -462,24 +452,17 @@ void gfs_poisson_coefficients (GfsDomain * domain,
for (i = 0; i < FTT_DIMENSION; i++) {
gdouble lambda = (&domain->lambda.x)[i];
- lambda2[i] = lambda*lambda;
- }
- if (alpha == NULL)
- gfs_domain_cell_traverse (domain,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) poisson_coeff, lambda2);
- else {
- gfs_domain_cell_traverse (domain,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) reset_coeff, NULL);
- gpointer data[2];
- data[0] = lambda2;
- data[1] = alpha;
- gfs_domain_face_traverse (domain, FTT_XYZ,
- FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttFaceTraverseFunc) poisson_alpha_coeff, data);
+ p.lambda2[i] = lambda*lambda;
}
gfs_domain_cell_traverse (domain,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) reset_coeff, NULL);
+ p.alpha = alpha;
+ p.domain = domain;
+ gfs_domain_face_traverse (domain, FTT_XYZ,
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttFaceTraverseFunc) poisson_coeff, &p);
+ gfs_domain_cell_traverse (domain,
FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
(FttCellTraverseFunc) face_coeff_from_below, NULL);
}
@@ -527,6 +510,7 @@ static void tension_coeff (FttCellFace * face, gpointer * data)
alpha, p.x, p.y, p.z);
}
v *= alpha;
+ /* fixme: mapping */
if (GFS_IS_MIXED (face->cell))
v *= s->solid->s[face->d];
s->f[face->d].v = v;
diff --git a/src/simulation.c b/src/simulation.c
index e77a807..d135c25 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -1780,3 +1780,47 @@ GfsSimulationClass * gfs_poisson_class (void)
return klass;
}
+
+/* GfsAxi: Object */
+
+static void axi_read (GtsObject ** object, GtsFile * fp)
+{
+ (* GTS_OBJECT_CLASS (gfs_axi_class ())->parent_class->read) (object, fp);
+ if (fp->type == GTS_ERROR)
+ return;
+
+ GFS_DOMAIN (*object)->refpos.y = 0.5;
+}
+
+static gdouble axi_face_fraction (const GfsDomain * domain, const FttCellFace * face)
+{
+ FttVector p;
+ ftt_face_pos (face, &p);
+ return p.y;
+}
+
+static void axi_class_init (GfsSimulationClass * klass)
+{
+ GTS_OBJECT_CLASS (klass)->read = axi_read;
+ GFS_DOMAIN_CLASS (klass)->face_fraction = axi_face_fraction;
+}
+
+GfsSimulationClass * gfs_axi_class (void)
+{
+ static GfsSimulationClass * klass = NULL;
+
+ if (klass == NULL) {
+ GtsObjectClassInfo gfs_axi_info = {
+ "GfsAxi",
+ sizeof (GfsSimulation),
+ sizeof (GfsSimulationClass),
+ (GtsObjectClassInitFunc) axi_class_init,
+ (GtsObjectInitFunc) NULL,
+ (GtsArgSetFunc) NULL,
+ (GtsArgGetFunc) NULL
+ };
+ klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_simulation_class ()), &gfs_axi_info);
+ }
+
+ return klass;
+}
diff --git a/src/simulation.h b/src/simulation.h
index 9292e55..493f628 100644
--- a/src/simulation.h
+++ b/src/simulation.h
@@ -135,6 +135,10 @@ GfsSimulationClass * gfs_advection_class (void);
GfsSimulationClass * gfs_poisson_class (void);
+/* GfsAxi: Header */
+
+GfsSimulationClass * gfs_axi_class (void);
+
#ifdef __cplusplus
}
#endif /* __cplusplus */
diff --git a/src/timestep.c b/src/timestep.c
index 9670cf3..f9978a0 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -50,6 +50,7 @@ static void scale_cell_gradients (FttCell * cell, gpointer * data)
guint * dimension = data[1];
FttComponent c;
+ /* fixme: mapping??? */
if (GFS_IS_MIXED (cell)) {
GfsSolidVector * s = GFS_STATE (cell)->solid;
@@ -116,9 +117,7 @@ static void correct_normal_velocity (FttCellFace * face,
dp = (g.b - g.a*GFS_VARIABLE (face->cell, p->i))/ftt_cell_size (face->cell);
if (!FTT_FACE_DIRECT (face))
dp = - dp;
-
- if (s->solid && s->solid->s[face->d] > 0.)
- dp /= s->solid->s[face->d];
+ dp /= gfs_domain_face_fraction (p->domain, face);
GFS_FACE_NORMAL_VELOCITY_LEFT (face) -= dp*(*dt);
if (gv)
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list