[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8
Stephane Popinet
popinet at users.sf.net
Tue Nov 24 12:25:12 UTC 2009
The following commit has been merged in the upstream branch:
commit 633df0dda71d90b1f46a861868a93d3d19130183
Author: Stephane Popinet <popinet at users.sf.net>
Date: Fri Sep 25 20:36:30 2009 +1000
Longitude-latitude metric
darcs-hash:20090925103630-d4795-fad6e66ec1da0522035aacab35aeb2a9f62eca91.gz
diff --git a/doc/figures/lonlat.fig b/doc/figures/lonlat.fig
new file mode 100644
index 0000000..5fb5b18
--- /dev/null
+++ b/doc/figures/lonlat.fig
@@ -0,0 +1,34 @@
+#FIG 3.2 Produced by xfig version 3.2.5
+Landscape
+Center
+Metric
+A4
+100.00
+Single
+-2
+1200 2
+5 1 0 2 0 7 50 -1 -1 0.000 0 1 0 0 3825.000 1117.500 2925 4050 3825 4185 4725 4050
+5 1 0 2 0 7 50 -1 -1 0.000 0 1 0 0 6652.500 5062.500 2925 4050 2790 5085 2925 6075
+5 1 0 2 0 7 50 -1 -1 0.000 0 1 0 0 8452.500 5062.500 4725 4050 4590 5040 4725 6075
+5 1 0 2 0 7 50 -1 -1 0.000 0 1 0 0 3825.000 3150.000 2925 6075 3780 6210 4725 6075
+5 1 0 1 0 7 50 -1 -1 0.000 0 1 0 0 5971.811 3405.259 5310 3150 5265 3465 5400 3825
+5 1 0 1 0 7 50 -1 -1 0.000 0 1 0 0 5038.176 2797.906 4680 3015 4905 3195 5265 3150
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 3
+ 2925 4050 5850 2250 2925 6075
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 7 0 0 2
+ 5850 2250 4725 4050
+2 1 1 1 0 7 50 -1 -1 4.000 0 0 7 0 0 2
+ 5850 2250 4725 6075
+4 0 0 50 -1 1 24 0.0000 4 180 165 3825 3330 r\001
+4 0 0 50 -1 1 24 0.0000 4 180 165 2970 6615 r\001
+4 0 0 50 -1 0 24 0.0000 4 180 555 3240 6615 cos\001
+4 0 0 50 -1 32 24 0.0000 4 285 210 3825 6660 q\001
+4 0 0 50 -1 1 24 0.0000 4 285 210 4140 6615 d\001
+4 0 0 50 -1 1 24 0.0000 4 180 165 1980 5175 r\001
+4 0 0 50 -1 32 24 0.0000 4 285 210 2340 5175 q\001
+4 0 0 50 -1 32 24 0.0000 4 300 225 4635 3465 l\001
+4 0 0 50 -1 32 24 0.0000 4 285 210 5085 4050 q\001
+4 0 0 50 -1 1 24 0.0000 4 285 210 2160 5175 d\001
+4 0 0 50 -1 1 24 0.0000 4 285 210 4455 3465 d\001
+4 0 0 50 -1 1 24 0.0000 4 285 210 4905 4050 d\001
+4 0 0 50 -1 32 24 0.0000 4 300 225 4320 6615 l\001
diff --git a/doc/figures/lonlat.tm b/doc/figures/lonlat.tm
new file mode 100644
index 0000000..d476007
--- /dev/null
+++ b/doc/figures/lonlat.tm
@@ -0,0 +1,361 @@
+<TeXmacs|1.0.6.11>
+
+<style|article>
+
+<\body>
+ <doc-data|<doc-title|Derivation of metric coefficients for a
+ <next-line>longitude-latitude coordinate
+ system>|<doc-author-data|<author-name|Stéphane
+ Popinet>>|<doc-date|<date>>|>
+
+ A system of conservation laws of the form
+
+ <\equation>
+ \<partial\><rsub|t><big|int><rsub|\<Omega\>>\<b-U\>*dV+<big|int><rsub|\<partial\>\<Omega\>>\<b-F\>(\<b-U\>)\<cdot\>\<b-n\>*dS=0<label|conservation>
+ </equation>
+
+ is discretised in Gerris as
+
+ <\equation*>
+ h<rsup|2>*c<rsub|i,j>*\<partial\><rsub|t>\<b-U\>+<big|sum><rsub|faces>s<rsub|f>*h*\<b-F\><rsub|f>(\<b-U\>)=0
+ </equation*>
+
+ or
+
+ <\equation>
+ h*c<rsub|i,j>*\<partial\><rsub|t>\<b-U\>+<big|sum><rsub|faces>s<rsub|f>\<b-F\><rsub|f>(\<b-U\>)*=0<label|gerris>
+ </equation>
+
+ with <math|h> the dimensional cell size. The figure below illustrates the
+ lengths of the faces of a surface element in a longitude-latitude
+ coordinate system with <math|\<lambda\>> the longitude and <math|\<theta\>>
+ the latitude.
+
+ <big-figure|<postscript|lonlat.fig|0.3par|||||>|Elementary lengths for a
+ surface element.>
+
+ This leads to the discrete form for (<reference|conservation>)
+
+ <\equation*>
+ r<rsup|2>*cos \<theta\>*d\<theta\>*d\<lambda\>*\<partial\><rsub|t>\<b-U\>+<big|sum><rsub|f<rsub|\<theta\>>>r*d\<theta\>*\<b-F\><rsub|f<rsub|\<theta\>>>(\<b-U\>)+<big|sum><rsub|f<rsub|\<lambda\>>>r*cos
+ \<theta\>*d\<lambda\>*\<b-F\><rsub|f<rsub|\<lambda\>>>(\<b-U\>)=0,
+ </equation*>
+
+ or
+
+ <\equation>
+ r*cos \<theta\>*d\<theta\>*d\<lambda\>*\<partial\><rsub|t>\<b-U\>+<big|sum><rsub|f<rsub|\<theta\>>>d\<theta\>*\<b-F\><rsub|f<rsub|\<theta\>>>(\<b-U\>)+<big|sum><rsub|f<rsub|\<lambda\>>>cos
+ \<theta\>*d\<lambda\>*\<b-F\><rsub|f<rsub|\<lambda\>>>(\<b-U\>)=0.<label|lonlat>
+ </equation>
+
+ Equating the terms with (<reference|gerris>) gives
+
+ <\with|color|blue>
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|c<rsub|i,j>>|<cell|\<equiv\>>|<cell|<frac|r|h>*cos
+ \<theta\> *d\<theta\>*d\<lambda\>,>>|<row|<cell|s<rsub|f<rsub|\<theta\>>>>|<cell|\<equiv\>>|<cell|d\<theta\>,>>|<row|<cell|s<rsub|f<rsub|\<lambda\>>>>|<cell|\<equiv\>>|<cell|cos
+ \<theta\>* d\<lambda\>.>>>>
+ </eqnarray*>
+ </with>
+
+ <subsection|Application to the Saint-Venant equations>
+
+ For the Saint-Venant equations <math|\<b-U\>> and <math|\<b-F\>(\<b-U\>)>
+ are given by
+
+ <\equation*>
+ \<b-U\>\<equiv\><matrix|<tformat|<table|<row|<cell|\<phi\>>>|<row|<cell|\<phi\>*u>>|<row|<cell|\<phi\>*v>>>>>,
+ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \<b-F\><rsub|x>(\<b-U\>)\<equiv\><matrix|<tformat|<table|<row|<cell|\<phi\>*u>>|<row|<cell|\<phi\>*u<rsup|2>+g*<frac|\<phi\><rsup|2>|2>>>|<row|<cell|\<phi\>*u*v>>>>>,
+ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \<b-F\><rsub|y>(\<b-U\>)\<equiv\><matrix|<tformat|<table|<row|<cell|\<phi\>*v>>|<row|<cell|\<phi\>*u*v>>|<row|<cell|\<phi\>*v<rsup|2>+g*<frac|\<phi\><rsup|2>|2>>>>>>.
+ </equation*>
+
+ Using (<reference|lonlat>), the equations in longitude-latitude coordinates
+ can be written
+
+ <\equation>
+ r*cos \<theta\>*\<partial\><rsub|t>\<b-U\>+<frac|1|d\<lambda\>><big|sum><rsub|f<rsub|x>>*\<b-F\><rsub|x>(\<b-U\>)+<frac|1|d\<theta\>><big|sum><rsub|f<rsub|y>>cos
+ \<theta\>*\<b-F\><rsub|y>(\<b-U\>)=0.<label|lonlat>
+ </equation>
+
+ The differential form can be recovered by taking the limit as
+ <math|d\<lambda\>> and <math|d\<theta\>> tend to zero
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|r*cos
+ \<theta\>*\<partial\><rsub|t>\<phi\>+\<partial\><rsub|\<lambda\>>(\<phi\>*u)+\<partial\><rsub|\<theta\>>(\<phi\>*v*cos
+ \<theta\>)=0,>|<cell|>>|<row|<cell|>|<cell|r*cos
+ \<theta\>*\<partial\><rsub|t>(\<phi\>*u)+\<partial\><rsub|\<lambda\>><left|(>\<phi\>*u<rsup|2>+g*<frac|\<phi\><rsup|2>|2><right|)>+\<partial\><rsub|\<theta\>><left|(>\<phi\>*u*v*cos
+ \<theta\><right|)>=0,>|<cell|>>|<row|<cell|>|<cell|r*cos
+ \<theta\>*\<partial\><rsub|t>(\<phi\>*v)+\<partial\><rsub|\<lambda\>>(\<phi\>*u*v)+\<partial\><rsub|\<theta\>><left|[><left|(>\<phi\>*v<rsup|2>+g*<frac|\<phi\><rsup|2>|2><right|)>*cos
+ \<theta\><right|]>=0.>|<cell|>>>>
+ </eqnarray*>
+
+ This can be rewritten in advective form as
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|\<partial\><rsub|t>\<phi\>+<frac|1|r*cos
+ \<theta\>*>*<left|[>\<partial\><rsub|\<lambda\>>(\<phi\>*u)+\<partial\><rsub|\<theta\>>(\<phi\>*v*cos
+ \<theta\>)<right|]>=0,>|<cell|>>|<row|<cell|>|<cell|\<partial\><rsub|t>u+<frac|1|r*cos
+ \<theta\>*>*(u*\<partial\><rsub|\<lambda\>>u+v*cos
+ \<theta\>*\<partial\><rsub|\<theta\>>
+ u)<with|color|blue|<with|color|green|-<frac|u*v|r>*tan
+ \<theta\>>>+<frac|g|r*cos \<theta\>*>*\<partial\><rsub|\<lambda\>>\<phi\>=0,>|<cell|>>|<row|<cell|>|<cell|\<partial\><rsub|t>v+<frac|1|r*cos
+ \<theta\>*>*(u*\<partial\><rsub|\<lambda\>>v+v*cos \<theta\>
+ \<partial\><rsub|\<theta\>> v)<with|color|green|+<frac|u<rsup|2>|r>*tan
+ \<theta\>>+<frac|g|r>*\<partial\><rsub|\<theta\>>\<phi\><with|color|red|-<frac|g*\<phi\>|2*r>*tan
+ \<theta\>>=0,>|<cell|>>>>
+ </eqnarray*>
+
+ where additional terms (compared to the standard form) are indicated in red
+ and missing terms in green.
+
+ <subsection|Derivation using the ``manifold approach'' of Rossmanith et
+ al., 2004>
+
+ System of ``balance laws'' on general manifolds (equation (95))
+
+ <\equation>
+ \<partial\><rsub|t>\<b-q\>+<frac|1|<sqrt|h>>*\<partial\><rsub|x<rsup|k>>\<b-f\><rsup|k>=<with|math-font-series|bold|\<psi\>><rsub|c>,<label|manifold>
+ </equation>
+
+ with <math|h> the determinant of the metric tensor and
+ <with|mode|math|<with|math-font-series|bold|\<psi\>><rsub|c>> the geometric
+ source terms. For the shallow-water equations
+
+ <\equation*>
+ \<b-q\>(\<b-x\>,t)\<equiv\><matrix|<tformat|<table|<row|<cell|\<phi\>>>|<row|<cell|\<phi\>*u<rsup|1>>>|<row|<cell|\<phi\>*u<rsup|2>>>>>>,<hspace|10pt>\<b-f\><rsup|k>(\<b-q\>,\<b-x\>)\<equiv\><sqrt|h>*<matrix|<tformat|<table|<row|<cell|\<phi\>*u<rsup|k>>>|<row|<cell|T<rsup|k1>>>|<row|<cell|T<rsup|k2>>>>>>,<hspace|10pt><with|math-font-series|bold|\<psi\>><rsub|c>(\<b-q\>,\<b-x\>)\<equiv\><matrix|<tformat|<table|<row|<cell|0>>|<row|<cell|-\<Gamma\><rsub|m
+ n><rsup|1>*T<rsup|n m>>>|<row|<cell|-\<Gamma\><rsub|m n><rsup|2>*T<rsup|n
+ m>>>>>>,
+ </equation*>
+
+ and
+
+ <\equation*>
+ T<rsup|n m>\<equiv\>\<phi\>*u<rsup|n>*u<rsup|m>+<frac|1|2>*g*\<phi\><rsup|2>*h<rsup|n
+ m>.
+ </equation*>
+
+ For a spherical coordinate system with <math|x<rsup|1>=\<lambda\>> the
+ longitude and <math|x<rsup|2>=\<theta\>> the latitude, the metric is
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|\<b-h\>>|<cell|\<equiv\>>|<cell|<matrix|<tformat|<table|<row|<cell|r<rsup|2>*cos<rsup|2>
+ \<theta\>>|<cell|0>>|<row|<cell|0>|<cell|r<rsup|2>>>>>>,>>|<row|<cell|<sqrt|h>>|<cell|\<equiv\>>|<cell|r<rsup|2>*cos
+ \<theta\>,>>>>
+ </eqnarray*>
+
+ so that
+
+ <\equation*>
+ h<rsup|11>=<frac|1|h<rsub|11>>=<frac|1|r<rsup|2>*cos<rsup|2>
+ \<theta\>>,<hspace|10pt>h<rsup|22>=<frac|1|h<rsub|22>>=<frac|1|r<rsup|2>>,<hspace|10pt>h<rsup|12>=h<rsup|21>=0.
+ </equation*>
+
+ The corresponding Christofell symbols are given by (for a general
+ orthogonal metric)
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|\<Gamma\><rsub|m
+ n><rsup|k>>|<cell|\<equiv\>>|<cell|<frac|1|2>*h<rsup|\<alpha\>k>(\<partial\><rsub|x<rsup|n>>h<rsub|\<alpha\>m>+\<partial\><rsub|x<rsup|m>>h<rsub|\<alpha\>n>-\<partial\><rsub|x<rsup|\<alpha\>>>h<rsub|m
+ n>),>>|<row|<cell|\<Gamma\><rsup|1>>|<cell|\<equiv\>>|<cell|<frac|1|2*h<rsub|11>>*<matrix|<tformat|<table|<row|<cell|\<partial\><rsub|x<rsup|1>>h<rsub|11>>|<cell|\<partial\><rsub|x<rsup|2>>h<rsub|11>>>|<row|<cell|\<partial\><rsub|x<rsup|2>>h<rsub|11>>|<cell|-\<partial\><rsub|x<rsup|1>>h<rsub|2
+ 2>>>>>>=<matrix|<tformat|<table|<row|<cell|0>|<cell|-tan
+ \<theta\>*>>|<row|<cell|-tan \<theta\>*>|<cell|0>>>>>,>>|<row|<cell|\<Gamma\><rsup|2>>|<cell|\<equiv\>>|<cell|<frac|1|2*h<rsub|22>>*<matrix|<tformat|<table|<row|<cell|-\<partial\><rsub|x<rsup|2>>h<rsub|11>>|<cell|\<partial\><rsub|x<rsup|1>>h<rsub|22>>>|<row|<cell|\<partial\><rsub|x<rsup|1>>h<rsub|22>>|<cell|\<partial\><rsub|x<rsup|2>>h<rsub|2
+ 2>>>>>>=<matrix|<tformat|<table|<row|<cell|*sin \<theta\>*cos
+ \<theta\>>|<cell|0>>|<row|<cell|0>|<cell|0>>>>>.>>>>
+ </eqnarray*>
+
+ The geometric source term is then
+
+ <\equation*>
+ <hspace|10pt><with|math-font-series|bold|\<psi\>><rsub|c>(\<b-q\>,\<b-x\>)\<equiv\><matrix|<tformat|<table|<row|<cell|0>>|<row|<cell|tan
+ \<theta\>*(T<rsup|1 2>+*T<rsup|2 1>)>>|<row|<cell|-sin \<theta\>*cos
+ \<theta\>*T<rsup|11>>>>>>=<matrix|<tformat|<table|<row|<cell|0>>|<row|<cell|2*tan
+ \<theta\>*\<phi\>*u<rsup|1>*u<rsup|2>>>|<row|<cell|-sin \<theta\>*cos
+ \<theta\>*\<phi\>*u<rsup|1>*u<rsup|1>-<frac|g*\<phi\><rsup|2>|2*r<rsup|2>*>*tan
+ \<theta\>>>>>>
+ </equation*>
+
+ Equation (<reference|manifold>) can be developed as
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|\<partial\><rsub|t>\<phi\>+<frac|1|cos
+ \<theta\>>*<left|[>\<partial\><rsub|\<lambda\>>(*cos
+ \<theta\>*\<phi\>*u<rsup|1>)+\<partial\><rsub|\<theta\>>(cos
+ \<theta\>*\<phi\>*u<rsup|2>)<right|]>>|<cell|=>|<cell|0,>>|<row|<cell|\<partial\><rsub|t>(\<phi\>*u<rsup|1>)+<frac|1|cos
+ \<theta\>>*<left|[>\<partial\><rsub|\<lambda\>>(cos
+ \<theta\>*(\<phi\>*u<rsup|1>*u<rsup|1>+<frac|g*\<phi\><rsup|2>|2*r<rsup|2>*cos<rsup|2>\<theta\>>))+\<partial\><rsub|\<theta\>>(cos
+ \<theta\>*\<phi\>*u<rsup|1>*u<rsup|2>)<right|]>>|<cell|=>|<cell|-\<Gamma\><rsub|m
+ n><rsup|1>*T<rsup|n m>,>>|<row|<cell|\<partial\><rsub|t>(\<phi\>*u<rsup|2>)+<frac|1|cos
+ \<theta\>>*<left|[>\<partial\><rsub|\<lambda\>>(cos
+ \<theta\>**\<phi\>*u<rsup|1>*u<rsup|2>)+\<partial\><rsub|\<theta\>>(cos
+ \<theta\>*(\<phi\>*u<rsup|2>*u<rsup|2>+<frac|g*\<phi\><rsup|2>|2*r<rsup|2>*>))<right|]>>|<cell|=>|<cell|-\<Gamma\><rsub|m
+ n><rsup|2>*T<rsup|n m>.>>>>
+ </eqnarray*>
+
+ Using
+
+ <\equation*>
+ u\<equiv\>u<rsub|1>\<equiv\>r*cos \<theta\>
+ u<rsup|1>,<hspace|10pt>v\<equiv\>u<rsub|2>\<equiv\>r*u<rsup|2>,
+ </equation*>
+
+ (where exponents and indices indicate the covariant and contravariant
+ vector coordinates respectively) then gives
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|\<partial\><rsub|t>\<phi\>+<frac|1|r*cos
+ \<theta\>>*<left|[>\<partial\><rsub|\<lambda\>>(\<phi\>*u)+\<partial\><rsub|\<theta\>>(cos
+ \<theta\>*\<phi\>*v)<right|]>>|<cell|=>|<cell|0,>>|<row|<cell|\<partial\><rsub|t>u+<frac|1|r*cos
+ \<theta\>*>*(u*\<partial\><rsub|\<lambda\>>u+cos
+ \<theta\>*v*\<partial\><rsub|\<theta\>>u)+<frac|g|r*cos
+ \<theta\>>*\<partial\><rsub|\<lambda\>>\<phi\>-<frac|u*v|r>*tan
+ \<theta\>>|<cell|=>|<cell|*0,>>|<row|<cell|\<partial\><rsub|t>v+<frac|1|r*cos
+ \<theta\>>*(u*\<partial\><rsub|\<lambda\>>v+cos
+ \<theta\>*v**\<partial\><rsub|\<theta\>>*v)+<frac|g|r>*\<partial\><rsub|\<theta\>>\<phi\>>|<cell|=>|<cell|-<frac|u*u|r*>*tan
+ \<theta\>,>>>>
+ </eqnarray*>
+
+ which is the standard advective form.
+
+ The geometric source term can be re-expressed in the covariant coordinate
+ system as
+
+ <\equation*>
+ <hspace|10pt><with|math-font-series|bold|<wide|\<psi\>|^>>\<equiv\><matrix|<tformat|<table|<row|<cell|0>>|<row|<cell|2*tan
+ \<theta\>*\<phi\>*<frac|u<rsub|1>*u<rsub|2>|r*>>>|<row|<cell|-tan
+ \<theta\>*\<phi\>*<frac|u<rsub|1>*u<rsub|1>|r*>-<frac|g*\<phi\><rsup|2>|2*r*>*tan
+ \<theta\>>>>>>
+ </equation*>
+
+ <subsection|The Saint-Venant equations in general orthogonal coordinates>
+
+ \;
+
+ <\equation*>
+ m<rsub|x>*m<rsub|y>*d*x*d*y*\<partial\><rsub|t>\<b-U\>+<big|sum><rsub|f<rsub|x>>m<rsub|y>*d*y*\<b-F\><rsub|x>(\<b-U\>)+<big|sum><rsub|f<rsub|y>>m<rsub|x>*d*x*\<b-F\><rsub|y>(\<b-U\>)=0
+ </equation*>
+
+ can be rewritten
+
+ <\equation*>
+ m<rsub|x>*m<rsub|y>*\<partial\><rsub|t>\<b-U\>+<frac|1|d*x>*<big|sum><rsub|f<rsub|x>>m<rsub|y>*\<b-F\><rsub|x>(\<b-U\>)+<frac|1|d*y>*<big|sum><rsub|f<rsub|y>>m<rsub|x>*\<b-F\><rsub|y>(\<b-U\>)=0.
+ </equation*>
+
+ The differential form can be recovered by making <math|d*x> and <math|d*y>
+ tend to zero
+
+ <\with|mode|math>
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|m<rsub|x>*m<rsub|y>*\<partial\><rsub|t>\<phi\>+\<partial\><rsub|x>(m<rsub|y>*\<phi\>*u)+\<partial\><rsub|y>(m<rsub|x>*\<phi\>*v)=0,>|<cell|>>|<row|<cell|>|<cell|m<rsub|x>*m<rsub|y>*\<partial\><rsub|t>(\<phi\>*u)+\<partial\><rsub|x><left|[>m<rsub|y>*<left|(>\<phi\>*u<rsup|2>+g*<frac|\<phi\><rsup|2>|2><right|)><right|]>+\<partial\><rsub|y><left|(>m<rsub|x>*\<phi\>*u*v<right|)>=0,>|<cell|>>|<row|<cell|>|<cell|m<rsub|x>*m<rsub|y>*\<partial\><rsub|t>(\<phi\>*v)+\<partial\><rsub|x>(m<rsub|y>*\<phi\>*u*v)+\<partial\><rsub|y><left|[>m<rsub|x>*<left|(>\<phi\>*v<rsup|2>+g*<frac|\<phi\><rsup|2>|2><right|)><right|]>=0.>|<cell|>>>>
+ </eqnarray*>
+ </with>
+
+ which can be expanded as
+
+ <\with|mode|math>
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|\<partial\><rsub|t>\<phi\>+<frac|u*|m<rsub|x>>*\<partial\><rsub|x>\<phi\>*+<frac|v*|m<rsub|y>>*\<partial\><rsub|y>\<phi\>+<frac|\<phi\>*|m<rsub|x>*m<rsub|y>>*<left|[>\<partial\><rsub|x>(m<rsub|y>*u)+\<partial\><rsub|y>(m<rsub|x>*v)<right|]>=0,>|<cell|>>|<row|<cell|>|<cell|\<partial\><rsub|t>u*+<frac|u*|m<rsub|x>>*\<partial\><rsub|x>u+<frac|v*|m<rsub|y>>*\<partial\><rsub|y>u*+<frac|g|m<rsub|x>>*\<partial\><rsub|x>\<phi\>+g*<frac|\<phi\>|2*m<rsub|x>*m<rsub|y>>*\<partial\><rsub|x>m<rsub|y>=0,>|<cell|>>|<row|<cell|>|<cell|\<partial\><rsub|t>v+<frac|u|m<rsub|x>>*\<partial\><rsub|x>v+<frac|v|m<rsub|y>>*\<partial\><rsub|y>v+<frac|g|m<rsub|y>>*\<partial\><rsub|y>\<phi\>+g*<frac|\<phi\>|2*m<rsub|x>*m<rsub|y>>*\<partial\><rsub|y>m<rsub|x>=0.>|<cell|>>>>
+ </eqnarray*>
+ </with>
+
+ Introducing the notation
+
+ <\equation*>
+ d<rsub|t>=\<partial\><rsub|t>+<frac|u*|m<rsub|x>>*\<partial\><rsub|x>+<frac|v*|m<rsub|y>>*\<partial\><rsub|y>,
+ </equation*>
+
+ this can be rewritten
+
+ <\with|mode|math>
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|d<rsub|t>\<phi\>+<frac|\<phi\>*|m<rsub|x>*m<rsub|y>>*<left|[>\<partial\><rsub|x>(m<rsub|y>*u)+\<partial\><rsub|y>(m<rsub|x>*v)<right|]>=0,>|<cell|>>|<row|<cell|>|<cell|d<rsub|t>u*<with|color|blue|<with|color|green|<with|color|green|-f<rsub|G>*v>>>+<frac|g|m<rsub|x>>*\<partial\><rsub|x>\<phi\><with|color|red|+<frac|g*\<phi\>|2*m<rsub|x>*m<rsub|y>>*\<partial\><rsub|x>m<rsub|y>>=0,>|<cell|>>|<row|<cell|>|<cell|d<rsub|t>v<with|color|green|<with|color|blue|<with|color|green|+f<rsub|G>*u>>>+<frac|g|m<rsub|y>>*\<partial\><rsub|y>\<phi\><with|color|red|+*<frac|g*\<phi\>|2*m<rsub|x>*m<rsub|y>>*\<partial\><rsub|y>m<rsub|x>>=0.>|<cell|>>>>
+ </eqnarray*>
+ </with>
+
+ where additional and missing terms (respective to equations 43, 44 and 45
+ of Williamson et al, 1991) are indicated in red and green respectively and
+
+ <\equation*>
+ f<rsub|G>\<equiv\><frac|v*\<partial\><rsub|x>m<rsub|y>-u*\<partial\><rsub|y>m<rsub|x>|m<rsub|x>*m<rsub|y>>*.*
+ </equation*>
+
+ <subsubsection|Application to spherical coordinates>
+
+ For spherical coordinates
+
+ <\equation*>
+ m<rsub|x>=r*cos y,<hspace|10pt>m<rsub|y>=r,<hspace|10pt>\<partial\><rsub|x>m<rsub|y>=0,<hspace|10pt>\<partial\><rsub|y>m<rsub|x>=-r*sin
+ y,
+ </equation*>
+
+ which gives
+
+ <\with|mode|math>
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|d<rsub|t>\<phi\>+<frac|\<phi\>*|r*cos
+ y>*<left|[>\<partial\><rsub|x>u+\<partial\><rsub|y>(v*cos
+ y)<right|]>=0,>|<cell|>>|<row|<cell|>|<cell|d<rsub|t>u*<with|color|blue|<with|color|green|<with|color|green|-<frac|u*v|r>*tan
+ y>>>+<frac|g|r*cos y>*\<partial\><rsub|x>\<phi\>=0,>|<cell|>>|<row|<cell|>|<cell|d<rsub|t>v<with|color|green|<with|color|blue|<with|color|green|<with|color|green|+<frac|u<rsup|2>|r>*tan
+ y>>>>+<frac|g|m<rsub|y>>*\<partial\><rsub|y>\<phi\><with|color|red|-<frac|g*\<phi\>|2*r>*tan
+ y>=0.>|<cell|>>>>
+ </eqnarray*>
+ </with>
+
+ <subsubsection|Application to polar coordinates>
+
+ For polar coordinates
+
+ <\equation*>
+ m<rsub|x>=1,<hspace|10pt>m<rsub|y>=x,<hspace|10pt>\<partial\><rsub|x>m<rsub|y>=1,<hspace|10pt>\<partial\><rsub|y>m<rsub|x>=0,
+ </equation*>
+
+ which gives
+
+ <\with|mode|math>
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|>|<cell|d<rsub|t>\<phi\>+<frac|\<phi\>*|x>*<left|[>\<partial\><rsub|x>(x*u)+\<partial\><rsub|y>v<right|]>=0,>|<cell|>>|<row|<cell|>|<cell|d<rsub|t>u*<with|color|blue|<with|color|green|<with|color|green|-<frac|v<rsup|2>|x>>>>+g*\<partial\><rsub|x>\<phi\><with|color|red|+g*<frac|\<phi\>|2*x>>=0,>|<cell|>>|<row|<cell|>|<cell|d<rsub|t>v<with|color|green|<with|color|blue|<with|color|green|<with|color|green|+<frac|u*v|x>>>>>+<frac|g|x>*\<partial\><rsub|y>\<phi\>=0.>|<cell|>>>>
+ </eqnarray*>
+ </with>
+</body>
+
+<\references>
+ <\collection>
+ <associate|auto-1|<tuple|1|1>>
+ <associate|auto-2|<tuple|1|2>>
+ <associate|auto-3|<tuple|2|2>>
+ <associate|auto-4|<tuple|3|?>>
+ <associate|auto-5|<tuple|3.1|?>>
+ <associate|auto-6|<tuple|3.2|?>>
+ <associate|conservation|<tuple|1|1>>
+ <associate|gerris|<tuple|2|1>>
+ <associate|lonlat|<tuple|4|1>>
+ <associate|manifold|<tuple|5|2>>
+ </collection>
+</references>
+
+<\auxiliary>
+ <\collection>
+ <\associate|figure>
+ <tuple|normal|Elementary lengths for a surface
+ element.|<pageref|auto-1>>
+ </associate>
+ <\associate|toc>
+ <with|par-left|<quote|1.5fn>|1<space|2spc>Application to the
+ Saint-Venant equations <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
+ <no-break><pageref|auto-2>>
+
+ <with|par-left|<quote|1.5fn>|2<space|2spc>Derivation using the
+ ``manifold approach'' of Rossmanith et al., 2004
+ <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
+ <no-break><pageref|auto-3>>
+
+ <with|par-left|<quote|1.5fn>|3<space|2spc>The Saint-Venant equations in
+ orthogonal curvilinear coordinates <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
+ <no-break><pageref|auto-4>>
+
+ <with|par-left|<quote|3fn>|3.1<space|2spc>Application to spherical
+ coordinates <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
+ <no-break><pageref|auto-5>>
+ </associate>
+ </collection>
+</auxiliary>
\ No newline at end of file
diff --git a/src/domain.h b/src/domain.h
index 356a464..b4c0bb4 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -73,15 +73,18 @@ struct _GfsDomain {
gpointer array;
gboolean overlap; /* whether to overlap MPI communications with computation */
+
+ /* coordinate metrics */
+ gpointer metric_data;
+ gdouble (* face_metric) (const GfsDomain *, const FttCellFace *, const gpointer);
+ gdouble (* cell_metric) (const GfsDomain *, const FttCell *, const gpointer);
+ gdouble (* solid_metric) (const GfsDomain *, const FttCell *, const gpointer);
};
struct _GfsDomainClass {
GtsWGraphClass parent_class;
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,\
@@ -339,14 +342,14 @@ GSList * gfs_receive_boxes (GfsDomain * domain,
* @face: a #FttCellFace.
*
* Returns: the surface fraction of @face taking into account any
- * orthogonal mapping of @domain.
+ * orthogonal metric 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_map)
- f *= (* GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->face_map) (domain, face);
+ if (domain->face_metric)
+ f *= (* domain->face_metric) (domain, face, domain->metric_data);
return f;
}
@@ -356,17 +359,17 @@ gdouble gfs_domain_face_fraction (const GfsDomain * domain, const FttCellFace *
* @face: a #FttCellFace.
*
* Returns: the surface fraction "to the right" of @face taking into account any
- * orthogonal mapping of @domain.
+ * orthogonal metric 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_map) {
+ if (domain->face_metric) {
FttCellFace face1;
face1.cell = face->neighbor;
face1.d = FTT_OPPOSITE_DIRECTION (face->d);
- f *= (* GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->face_map) (domain, &face1);
+ f *= (* domain->face_metric) (domain, &face1, domain->metric_data);
}
return f;
}
@@ -377,30 +380,30 @@ gdouble gfs_domain_face_fraction_right (const GfsDomain * domain, const FttCellF
* @cell: a #FttCell.
*
* Returns: the volume fraction of @cell taking into account any
- * orthogonal mapping of @domain.
+ * orthogonal metric of @domain.
*/
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_map)
- a *= (* GFS_DOMAIN_CLASS (GTS_OBJECT (domain)->klass)->cell_map) (domain, cell);
+ if (domain->cell_metric)
+ a *= (* domain->cell_metric) (domain, cell, domain->metric_data);
return a;
}
/**
- * gfs_domain_solid_map:
+ * gfs_domain_solid_metric:
* @domain; a #GfsDomain.
* @cell: a mixed #FttCell.
*
- * Returns: the coordinate mapping at the center of area of the solid
+ * Returns: the coordinate metric 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)
+gdouble gfs_domain_solid_metric (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);
+ if (domain->solid_metric)
+ return (* domain->solid_metric) (domain, cell, domain->metric_data);
return 1.;
}
diff --git a/src/init.c b/src/init.c
index b3d98a3..c5f0ea7 100644
--- a/src/init.c
+++ b/src/init.c
@@ -164,6 +164,8 @@ GtsObjectClass ** gfs_classes (void)
#endif /* FTT_2D */
gfs_init_wave_class (),
+ gfs_metric_lon_lat_class (),
+
gfs_adapt_class (),
gfs_adapt_vorticity_class (),
gfs_adapt_streamline_curvature_class (),
diff --git a/src/poisson.c b/src/poisson.c
index 0115a21..8dc8d21 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -793,7 +793,7 @@ 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_domain_solid_map (c->domain, cell)*gfs_source_diffusion_cell (c->d, cell);
+ c->dt*gfs_domain_solid_metric (c->domain, cell)*gfs_source_diffusion_cell (c->d, cell);
if (c->rhoc) {
gdouble rho = c->alpha ? 1./gfs_function_value (c->alpha, cell) : 1.;
if (rho <= 0.) {
diff --git a/src/river.c b/src/river.c
index e2bb06e..e18d571 100644
--- a/src/river.c
+++ b/src/river.c
@@ -193,6 +193,24 @@ static void sources (FttCell * cell, GfsRiver * r)
zbR = GFS_VALUE (cell, r->v[3]) + GFS_VALUE (cell, r->dv[1][3]);
GFS_VALUE (cell, r->v[2]) += r->dt*r->g/2.*(etaL + etaR)*(zbL - zbR)/delta;
+
+ /* longitude-latitude geometric source terms */
+#if 0
+ FttVector p;
+ ftt_cell_pos (cell, &p);
+ gdouble radius = 3./(2.*M_PI), g = 1.;
+ gdouble theta = p.y/radius, tantheta = tan (theta);
+ gdouble
+ phiu = GFS_VALUE (cell, r->v[1]),
+ phiv = GFS_VALUE (cell, r->v[2]),
+ phi = GFS_VALUE (cell, r->v[0]);
+#if 1
+ GFS_VALUE (cell, r->v[1]) += r->dt*tantheta/radius*phiu*phiv/phi;
+ GFS_VALUE (cell, r->v[2]) += r->dt*tantheta/radius*(- phiu*phiu/phi - g*phi*phi/2.);
+#else
+ GFS_VALUE (cell, r->v[2]) += r->dt*tantheta/radius*(- g*phi*phi/2.);
+#endif
+#endif
}
static void advance (GfsRiver * r, gdouble dt)
diff --git a/src/simulation.c b/src/simulation.c
index bcd3b8f..9f453df 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -1968,21 +1968,24 @@ static void axi_read (GtsObject ** object, GtsFile * fp)
GFS_DOMAIN (*object)->refpos.y = 0.5;
}
-static gdouble axi_face_map (const GfsDomain * domain, const FttCellFace * face)
+static gdouble axi_face_metric (const GfsDomain * domain, const FttCellFace * face,
+ const gpointer data)
{
FttVector p;
ftt_face_pos (face, &p);
return p.y;
}
-static gdouble axi_cell_map (const GfsDomain * domain, const FttCell * cell)
+static gdouble axi_cell_metric (const GfsDomain * domain, const FttCell * cell,
+ const gpointer data)
{
FttVector p;
gfs_cell_cm (cell, &p);
return p.y;
}
-static gdouble axi_solid_map (const GfsDomain * domain, const FttCell * cell)
+static gdouble axi_solid_metric (const GfsDomain * domain, const FttCell * cell,
+ const gpointer data)
{
g_assert (GFS_IS_MIXED (cell));
return GFS_STATE (cell)->solid->ca.y;
@@ -1991,9 +1994,13 @@ static gdouble axi_solid_map (const GfsDomain * domain, const FttCell * cell)
static void axi_class_init (GfsSimulationClass * klass)
{
GTS_OBJECT_CLASS (klass)->read = axi_read;
- 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;
+}
+
+static void axi_init (GfsDomain * domain)
+{
+ domain->face_metric = axi_face_metric;
+ domain->cell_metric = axi_cell_metric;
+ domain->solid_metric = axi_solid_metric;
}
GfsSimulationClass * gfs_axi_class (void)
@@ -2006,7 +2013,7 @@ GfsSimulationClass * gfs_axi_class (void)
sizeof (GfsSimulation),
sizeof (GfsSimulationClass),
(GtsObjectClassInitFunc) axi_class_init,
- (GtsObjectInitFunc) NULL,
+ (GtsObjectInitFunc) axi_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
@@ -2015,3 +2022,96 @@ GfsSimulationClass * gfs_axi_class (void)
return klass;
}
+
+/* GfsMetricLonLat: Object */
+
+static void metric_lon_lat_write (GtsObject * o, FILE * fp)
+{
+ (* GTS_OBJECT_CLASS (gfs_metric_lon_lat_class ())->parent_class->write) (o, fp);
+ fprintf (fp, " %g", GFS_METRIC_LON_LAT (o)->r);
+}
+
+static gdouble metric_lon_lat_face_metric (const GfsDomain * domain, const FttCellFace * face,
+ const gpointer data)
+{
+ if (face->d/2 != FTT_Y)
+ return 1.;
+ FttVector p;
+ ftt_face_pos (face, &p);
+ return cos (p.y*GFS_SIMULATION (domain)->physical_params.L/GFS_METRIC_LON_LAT (data)->r);
+}
+
+static gdouble metric_lon_lat_cell_metric (const GfsDomain * domain, const FttCell * cell,
+ const gpointer data)
+{
+ FttVector p;
+ gfs_cell_cm (cell, &p);
+ return cos (p.y*GFS_SIMULATION (domain)->physical_params.L/GFS_METRIC_LON_LAT (data)->r);
+}
+
+static gdouble metric_lon_lat_solid_metric (const GfsDomain * domain, const FttCell * cell,
+ const gpointer data)
+{
+ g_assert (GFS_IS_MIXED (cell));
+ g_assert_not_implemented ();
+ return 1.;
+}
+
+static void metric_lon_lat_read (GtsObject ** o, GtsFile * fp)
+{
+ (* GTS_OBJECT_CLASS (gfs_metric_lon_lat_class ())->parent_class->read) (o, fp);
+ if (fp->type == GTS_ERROR)
+ return;
+
+ GfsDomain * domain = GFS_DOMAIN (gfs_object_simulation (*o));
+ if (domain->metric_data || domain->face_metric || domain->cell_metric || domain->solid_metric) {
+ gts_file_error (fp, "cannot use multiple metrics (yet)");
+ return;
+ }
+
+ GFS_METRIC_LON_LAT (*o)->r = gfs_read_constant (fp, gfs_object_simulation (*o));
+ if (fp->type == GTS_ERROR)
+ return;
+ if (GFS_METRIC_LON_LAT (*o)->r <= 0.) {
+ gts_file_error (fp, "radius must be strictly positive");
+ return;
+ }
+
+ domain->metric_data = *o;
+ domain->face_metric = metric_lon_lat_face_metric;
+ domain->cell_metric = metric_lon_lat_cell_metric;
+ domain->solid_metric = metric_lon_lat_solid_metric;
+}
+
+static void metric_lon_lat_class_init (GtsObjectClass * klass)
+{
+ klass->read = metric_lon_lat_read;
+ klass->write = metric_lon_lat_write;
+}
+
+static void metric_lon_lat_init (GfsMetricLonLat * m)
+{
+ GFS_EVENT (m)->istep = G_MAXINT/2;
+ m->r = 1.;
+}
+
+GfsEventClass * gfs_metric_lon_lat_class (void)
+{
+ static GfsEventClass * klass = NULL;
+
+ if (klass == NULL) {
+ GtsObjectClassInfo gfs_metric_lon_lat_info = {
+ "GfsMetricLonLat",
+ sizeof (GfsMetricLonLat),
+ sizeof (GfsEventClass),
+ (GtsObjectClassInitFunc) metric_lon_lat_class_init,
+ (GtsObjectInitFunc) metric_lon_lat_init,
+ (GtsArgSetFunc) NULL,
+ (GtsArgGetFunc) NULL
+ };
+ klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_event_class ()),
+ &gfs_metric_lon_lat_info);
+ }
+
+ return klass;
+}
diff --git a/src/simulation.h b/src/simulation.h
index f01c08d..544ac95 100644
--- a/src/simulation.h
+++ b/src/simulation.h
@@ -147,6 +147,26 @@ GfsSimulationClass * gfs_poisson_class (void);
GfsSimulationClass * gfs_axi_class (void);
+/* GfsMetricLonLat: Header */
+
+typedef struct _GfsMetricLonLat GfsMetricLonLat;
+
+struct _GfsMetricLonLat {
+ /*< private >*/
+ GfsEvent parent;
+
+ /*< public >*/
+ gdouble r;
+};
+
+#define GFS_METRIC_LON_LAT(obj) GTS_OBJECT_CAST (obj,\
+ GfsMetricLonLat,\
+ gfs_metric_lon_lat_class ())
+#define GFS_IS_METRIC_LON_LAT(obj) (gts_object_is_from_class (obj,\
+ gfs_metric_lon_lat_class ()))
+
+GfsEventClass * gfs_metric_lon_lat_class (void);
+
#ifdef __cplusplus
}
#endif /* __cplusplus */
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list