[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8
Stephane Popinet
popinet at users.sf.net
Tue Nov 24 12:24:54 UTC 2009
The following commit has been merged in the upstream branch:
commit 66a9b31fe5ff607276464bc15b98a45d2b2c85fb
Author: Stephane Popinet <popinet at users.sf.net>
Date: Wed Aug 5 19:47:56 2009 +1000
GSE alleviation for GfsWave wind wave model
darcs-hash:20090805094756-d4795-b37c0318daee3473542fcd162e60221904e5c339.gz
diff --git a/doc/figures/gse.tm b/doc/figures/gse.tm
new file mode 100644
index 0000000..aea0ccd
--- /dev/null
+++ b/doc/figures/gse.tm
@@ -0,0 +1,252 @@
+<TeXmacs|1.0.6.11>
+
+<style|<tuple|article|maxima>>
+
+<\body>
+ <section|GSE alleviation using spatial filtering (or is it diffusion?)>
+
+ Following Tolman, 2002, Appendix A, the filtered value is defined as
+
+ <\equation>
+ F<rsub|avg>(\<b-x\>)=<frac|1|3>*F(\<b-x\>)+<frac|1|6>*<big|sum><rsup|4><rsub|n=1>F(\<b-x\>+\<b-r\><rsub|n>),<label|avg>
+ </equation>
+
+ with
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|\<b-r\><rsub|1>>|<cell|=>|<cell|\<b-s\>+\<b-n\>,>>|<row|<cell|\<b-r\><rsub|2>>|<cell|=>|<cell|-\<b-s\>+\<b-n\>,>>|<row|<cell|\<b-r\><rsub|3>>|<cell|=>|<cell|-\<b-s\>-\<b-n\>,>>|<row|<cell|\<b-r\><rsub|4>>|<cell|=>|<cell|\<b-s\>-\<b-n\>>>>>
+ </eqnarray*>
+
+ and
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|<with|math-font-series|bold|s>>|<cell|=>|<cell|\<alpha\>*<matrix|<tformat|<table|<row|<cell|cos
+ \<theta\>>>|<row|<cell|sin \<theta\>>>>>>,>>|<row|<cell|\<b-n\>>|<cell|=>|<cell|\<beta\>*<matrix|<tformat|<table|<row|<cell|-sin
+ \<theta\>>>|<row|<cell|cos \<theta\>>>>>>.>>>>
+ </eqnarray*>
+
+ To third-order the Taylor expansion of <math|F> can be written
+
+ <\equation>
+ F(\<b-x\>+\<b-r\>)=F(\<b-x\>)+\<b-r\>*\<b-F\><rprime|'>(\<b-x\>)+<frac|1|2>*\<b-r\><rsup|T>*\<b-F\><rprime|''>(\<b-x\>)*\<b-r\>,<label|taylor>
+ </equation>
+
+ with the gradient <math|\<b-F\><rprime|'>> and Hessian
+ <math|\<b-F\><rprime|''>> given by
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|<with|mode|text|<math|\<b-F\><rprime|'>>>>|<cell|\<equiv\>>|<cell|<matrix|<tformat|<table|<row|<cell|<frac|\<partial\>F|\<partial\>x>>>|<row|<cell|<frac|\<partial\>F|\<partial\>y>>>>>>,>>|<row|<cell|<with|mode|text|<math|\<b-F\><rprime|''>>>>|<cell|\<equiv\>>|<cell|<matrix|<tformat|<table|<row|<cell|<frac|\<partial\><rsup|2>F|\<partial\>x<rsup|2>>>|<cell|<frac|\<partial\><rsup|2>F|\<partial\>x\<partial\>y>>>|<row|<cell|<frac|\<partial\><rsup|2>F|\<partial\>x\<partial\>y>>|<cell|<frac|\<partial\><rsup|2>F|\<partial\>y<rsup|2>>>>>>>.>>>>
+ </eqnarray*>
+
+ Using (<reference|taylor>) in (<reference|avg>) gives
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|F<rsub|avg>(\<b-x\>)>|<cell|=>|<cell|<frac|1|3>*F(\<b-x\>)+<frac|1|6>*<big|sum><rsup|4><rsub|n=1>F(\<b-x\>)+\<b-r\><rsub|n>*\<b-F\><rprime|'>(\<b-x\>)+<frac|1|2>*\<b-r\><rsub|n><rsup|T>*\<b-F\><rprime|''>(\<b-x\>)*\<b-r\><rsub|n>>>|<row|<cell|>|<cell|=>|<cell|F(\<b-x\>)+<frac|1|6>*\<b-F\><rprime|'>(\<b-x\>)*<big|sum><rsup|4><rsub|n=1>\<b-r\><rsub|n>+<frac|1|12>*<big|sum><rsup|4><rsub|n=1>\<b-r\><rsub|n><rsup|T>*\<b-F\><rprime|''>(\<b-x\>)*\<b-r\><rsub|n>>>|<row|<cell|>|<cell|=>|<cell|F(\<b-x\>)+<frac|1|12>*<big|sum><rsup|4><rsub|n=1>\<b-r\><rsub|n><rsup|T>*\<b-F\><rprime|''>(\<b-x\>)*\<b-r\><rsub|n>>>>>
+ </eqnarray*>
+
+ I am lazy so I use Maxima to obtain the simplified form
+
+ <with|prog-language|maxima|prog-session|default|<\session>
+ <\input|<with|color|red|(<with|math-font-family|rm|%i>1)
+ <with|color|black|>>>
+ s:matrix([alpha*cos(theta)],[alpha*sin(theta)]);
+ </input>
+
+ <\output>
+ <with|mode|math|math-display|true|<with|mode|text|font-family|tt|color|red|(<with|math-font-family|rm|%o1>)
+ <with|color|black|>><left|(><tabular*|<tformat|<table|<row|<cell|\<alpha\>*cos
+ <left|(>\<vartheta\><right|)>>>|<row|<cell|\<alpha\>*sin
+ <left|(>\<vartheta\><right|)>>>>>><right|)>>
+ </output>
+
+ <\input|<with|color|red|(<with|math-font-family|rm|%i>2)
+ <with|color|black|>>>
+ n:matrix([-beta*sin(theta)],[beta*cos(theta)]);
+ </input>
+
+ <\output>
+ <with|mode|math|math-display|true|<with|mode|text|font-family|tt|color|red|(<with|math-font-family|rm|%o2>)
+ <with|color|black|>><left|(><tabular*|<tformat|<table|<row|<cell|-\<beta\>*sin
+ <left|(>\<vartheta\><right|)>>>|<row|<cell|\<beta\>*cos
+ <left|(>\<vartheta\><right|)>>>>>><right|)>>
+ </output>
+
+ <\input|<with|color|red|(<with|math-font-family|rm|%i>3)
+ <with|color|black|>>>
+ F:matrix( [Fxx,Fxy], [Fxy,Fyy] );
+ </input>
+
+ <\output>
+ <with|mode|math|math-display|true|<with|mode|text|font-family|tt|color|red|(<with|math-font-family|rm|%o3>)
+ <with|color|black|>><left|(><tabular*|<tformat|<table|<row|<cell|<with|math-font-family|rm|Fxx>>|<cell|<with|math-font-family|rm|Fxy>>>|<row|<cell|<with|math-font-family|rm|Fxy>>|<cell|<with|math-font-family|rm|Fyy>>>>>><right|)>>
+ </output>
+
+ <\input|<with|color|red|(<with|math-font-family|rm|%i>4)
+ <with|color|black|>>>
+ r1: s+n;
+ </input>
+
+ <\output>
+ <with|mode|math|math-display|true|<with|mode|text|font-family|tt|color|red|(<with|math-font-family|rm|%o4>)
+ <with|color|black|>><left|(><tabular*|<tformat|<table|<row|<cell|\<alpha\>*cos
+ <left|(>\<vartheta\><right|)>-\<beta\>*sin
+ <left|(>\<vartheta\><right|)>>>|<row|<cell|\<alpha\>*sin
+ <left|(>\<vartheta\><right|)>+\<beta\>*cos
+ <left|(>\<vartheta\><right|)>>>>>><right|)>>
+ </output>
+
+ <\input|<with|color|red|(<with|math-font-family|rm|%i>5)
+ <with|color|black|>>>
+ r2: -s+n;
+ </input>
+
+ <\output>
+ <with|mode|math|math-display|true|<with|mode|text|font-family|tt|color|red|(<with|math-font-family|rm|%o5>)
+ <with|color|black|>><left|(><tabular*|<tformat|<table|<row|<cell|-\<beta\>*sin
+ <left|(>\<vartheta\><right|)>-\<alpha\>*cos
+ <left|(>\<vartheta\><right|)>>>|<row|<cell|\<beta\>*cos
+ <left|(>\<vartheta\><right|)>-\<alpha\>*sin
+ <left|(>\<vartheta\><right|)>>>>>><right|)>>
+ </output>
+
+ <\input|<with|color|red|(<with|math-font-family|rm|%i>6)
+ <with|color|black|>>>
+ r3: -s-n;
+ </input>
+
+ <\output>
+ <with|mode|math|math-display|true|<with|mode|text|font-family|tt|color|red|(<with|math-font-family|rm|%o6>)
+ <with|color|black|>><left|(><tabular*|<tformat|<table|<row|<cell|\<beta\>*sin
+ <left|(>\<vartheta\><right|)>-\<alpha\>*cos
+ <left|(>\<vartheta\><right|)>>>|<row|<cell|-\<alpha\>*sin
+ <left|(>\<vartheta\><right|)>-\<beta\>*cos
+ <left|(>\<vartheta\><right|)>>>>>><right|)>>
+ </output>
+
+ <\input|<with|color|red|(<with|math-font-family|rm|%i>7)
+ <with|color|black|>>>
+ r4: s-n;
+ </input>
+
+ <\output>
+ <with|mode|math|math-display|true|<with|mode|text|font-family|tt|color|red|(<with|math-font-family|rm|%o10>)
+ <with|color|black|>><left|(><tabular*|<tformat|<table|<row|<cell|\<beta\>*sin
+ <left|(>\<vartheta\><right|)>+\<alpha\>*cos
+ <left|(>\<vartheta\><right|)>>>|<row|<cell|\<alpha\>*sin
+ <left|(>\<vartheta\><right|)>-\<beta\>*cos
+ <left|(>\<vartheta\><right|)>>>>>><right|)>>
+ </output>
+
+ <\input|<with|color|red|(<with|math-font-family|rm|%i>11)
+ <with|color|black|>>>
+ ratsimp(transpose(r1).F.r1+transpose(r2).F.r2+transpose(r3).F.r3+transpose(r4).F.r4);
+ </input>
+
+ <\output>
+ <with|mode|math|math-display|true|<with|mode|text|font-family|tt|color|red|(<with|math-font-family|rm|%o12>)
+ <with|color|black|>><left|(>4*\<alpha\><rsup|2>*<with|math-font-family|rm|Fyy>+4*\<beta\><rsup|2>*<with|math-font-family|rm|Fxx><right|)>*sin
+ <left|(>\<vartheta\><right|)><rsup|2>+<left|(>8*\<alpha\><rsup|2>-8*\<beta\><rsup|2><right|)>*<with|math-font-family|rm|Fxy>*cos
+ <left|(>\<vartheta\><right|)>*sin <left|(>\<vartheta\><right|)>+<left|(>4*\<beta\><rsup|2>*<with|math-font-family|rm|Fyy>+4*\<alpha\><rsup|2>*<with|math-font-family|rm|Fxx><right|)>*cos
+ <left|(>\<vartheta\><right|)><rsup|2>>
+ </output>
+ </session>>
+
+ We get
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|F<rsub|avg>(\<b-x\>)>|<cell|=>|<cell|F(\<b-x\>)+>>|<row|<cell|>|<cell|>|<cell|<frac|1|3>**(\<alpha\><rsup|2>*cos<rsup|2>
+ \<theta\>+\<beta\><rsup|2>*sin<rsup|2> \<theta\>)*F<rsub|x
+ x>*+>>|<row|<cell|>|<cell|>|<cell|<frac|1|3>**(\<alpha\><rsup|2>*sin<rsup|2>
+ \<theta\>+\<beta\><rsup|2>*cos<rsup|2> \<theta\>)*F<rsub|y
+ y>+>>|<row|<cell|>|<cell|>|<cell|<frac|2|3>
+ (\<alpha\><rsup|2>-\<beta\><rsup|2>)*cos \<theta\>*sin \<theta\>*F<rsub|x
+ y>,>>>>
+ </eqnarray*>
+
+ which can be rewritten
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|<frac|F<rsub|avg>(\<b-x\>)-F(\<b-x\>)|\<Delta\>t>>|<cell|=>|<cell|D<rsub|x
+ x>*F<rsub|x x>+D<rsub|y y>*F<rsub|y y>+2*D<rsub|x y>*F<rsub|x y>,>>>>
+ </eqnarray*>
+
+ with
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|D<rsub|x*x>>|<cell|\<equiv\>>|<cell|D<rsub|s*s>*cos<rsup|2>
+ \<theta\>+D<rsub|n*n>*sin<rsup|2> \<theta\>,>>|<row|<cell|D<rsub|y*y>>|<cell|\<equiv\>>|<cell|D<rsub|s*s>*sin<rsup|2>
+ \<theta\>+D<rsub|n*n>*cos<rsup|2> \<theta\>,>>|<row|<cell|D<rsub|x*y>>|<cell|\<equiv\>>|<cell|(D<rsub|s*s>-D<rsub|n*n>)*cos
+ \<theta\>*sin \<theta\>,>>>>
+ </eqnarray*>
+
+ and
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|D<rsub|s*s>>|<cell|\<equiv\>>|<cell|<frac|\<alpha\><rsup|2>|3*\<Delta\>t>,>>|<row|<cell|D<rsub|n*n>>|<cell|\<equiv\>>|<cell|<frac|\<beta\><rsup|2>|3*\<Delta\>t>.>>>>
+ </eqnarray*>
+
+ This means that to third-order accuracy the spatial-filtering scheme of
+ Tolman is formally equivalent to a first-order time discretisation of the
+ diffusion equation of Booij and Holthuijsen, 1987. Furthermore Tolman takes
+ <math|\<alpha\>> and <math|\<beta\>> as
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|\<alpha\>>|<cell|\<equiv\>>|<cell|\<alpha\><rsub|s>*\<Delta\>c<rsub|g>*\<Delta\>t,>>|<row|<cell|\<beta\>>|<cell|\<equiv\>>|<cell|\<alpha\><rsub|n>*c<rsub|g>*\<Delta\>\<theta\>*\<Delta\>t,>>>>
+ </eqnarray*>
+
+ with <math|\<Delta\>c<rsub|g>\<equiv\>><with|mode|math|(\<gamma\>-\<gamma\><rsup|-1>)*c<rsub|g>/2>,
+ whereas Booij and Holthuijsen take
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|<wide|D|^><rsub|s*s>>|<cell|\<equiv\>>|<cell|<frac|1|12>*(\<Delta\>c<rsub|g>)<rsup|2>*T<rsub|s>,>>|<row|<cell|<wide|D|^><rsub|n*n>>|<cell|\<equiv\>>|<cell|<frac|1|12>*(c<rsub|g>*\<Delta\>\<theta\>)<rsup|2>*T<rsub|s>.>>>>
+ </eqnarray*>
+
+ The two schemes are actually equivalent when <math|D<rsub|s*s>>,
+ <math|D<rsub|n*n>> are identical to <math|<wide|D|^><rsub|s*s>>,
+ <math|<wide|D|^><rsub|n*n>> which gives
+
+ <\eqnarray*>
+ <tformat|<table|<row|<cell|T<rsub|s>>|<cell|=>|<cell|4*\<alpha\><rsup|2><rsub|s>*\<Delta\>t,>>|<row|<cell|\<alpha\><rsub|n>>|<cell|=>|<cell|\<alpha\><rsub|s>.>>>>
+ </eqnarray*>
+
+ This is obviously a much smaller value for <math|T<rsub|s>> than the 4 days
+ used in Tolman 2002 and I am not sure I understand where this discrepancy
+ comes from... The stability criterion for the explicit diffusion scheme
+ then becomes
+
+ <\equation*>
+ \<Delta\>t\<less\><frac|(\<Delta\>x)<rsup|2>|(c<rsub|g>*\<Delta\>\<theta\>)<rsup|2>*4*\<alpha\><rsup|2><rsub|s>*\<Delta\>t>
+ </equation*>
+
+ which simplifies as
+
+ <\equation*>
+ \<Delta\>t\<less\><frac|\<Delta\>x|2*\<alpha\><rsub|s>*c<rsub|g>*\<Delta\>\<theta\>>
+ </equation*>
+
+ which is the standard CFL stability criterion with a CFL number of
+ <math|(2*\<alpha\><rsub|s>*\<Delta\>\<theta\>)<rsup|-1>>. For 24 wave
+ directions this number is larger than one for
+ <math|\<alpha\><rsub|s>\<lesssim\>2>, which means that this scheme does not
+ restrict the timestep if wave advection is resolved using an explicit
+ scheme (as is usually the case).
+</body>
+
+<\references>
+ <\collection>
+ <associate|auto-1|<tuple|1|1>>
+ <associate|avg|<tuple|1|1>>
+ <associate|taylor|<tuple|2|1>>
+ </collection>
+</references>
+
+<\auxiliary>
+ <\collection>
+ <\associate|toc>
+ <vspace*|1fn><with|font-series|<quote|bold>|math-font-series|<quote|bold>|1<space|2spc>GSE
+ alleviation using spatial filtering (or is it diffusion?)>
+ <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-1><vspace|0.5fn>
+ </associate>
+ </collection>
+</auxiliary>
\ No newline at end of file
diff --git a/src/adaptive.c b/src/adaptive.c
index a22415c..5bb81cc 100644
--- a/src/adaptive.c
+++ b/src/adaptive.c
@@ -542,7 +542,7 @@ GfsEventClass * gfs_adapt_gradient_class (void)
static void gfs_adapt_error_destroy (GtsObject * o)
{
- if (GFS_ADAPT_ERROR (o)->v)
+ if (GFS_ADAPT_ERROR (o)->v != GFS_ADAPT (o)->c)
gts_object_destroy (GTS_OBJECT (GFS_ADAPT_ERROR (o)->v));
(* GTS_OBJECT_CLASS (gfs_adapt_error_class ())->parent_class->destroy) (o);
@@ -554,60 +554,16 @@ static void gfs_adapt_error_read (GtsObject ** o, GtsFile * fp)
if (fp->type == GTS_ERROR)
return;
- GFS_ADAPT_ERROR (*o)->v = gfs_temporary_variable (GFS_DOMAIN (gfs_object_simulation (*o)));
+ GFS_ADAPT_ERROR (*o)->v = GFS_ADAPT (*o)->c ? GFS_ADAPT (*o)->c :
+ gfs_temporary_variable (GFS_DOMAIN (gfs_object_simulation (*o)));
GFS_ADAPT_ERROR (*o)->v->coarse_fine = none;
GFS_ADAPT_ERROR (*o)->v->fine_coarse = none;
}
-static gdouble center_regular_gradient (FttCell * cell, FttComponent c, GfsVariable * v)
-{
- FttCell * n1 = ftt_cell_neighbor (cell, 2*c);
- guint level = ftt_cell_level (cell);
- if (n1) {
- if (ftt_cell_level (n1) < level)
- return center_regular_gradient (ftt_cell_parent (cell), c, v)/2.;
- FttCell * n2 = ftt_cell_neighbor (cell, 2*c + 1);
- if (n2) {
- if (ftt_cell_level (n2) < level)
- return center_regular_gradient (ftt_cell_parent (cell), c, v)/2.;
- /* two neighbors: second-order differencing (parabola) */
- return (GFS_VALUE (n1, v) - GFS_VALUE (n2, v))/2.;
- }
- else
- /* one neighbor: first-order differencing */
- return GFS_VALUE (n1, v) - GFS_VALUE (cell, v);
- }
- else {
- FttCell * n2 = ftt_cell_neighbor (cell, 2*c + 1);
- if (n2) {
- if (ftt_cell_level (n2) < level)
- return center_regular_gradient (ftt_cell_parent (cell), c, v)/2.;
- /* one neighbor: first-order differencing */
- return GFS_VALUE (cell, v) - GFS_VALUE (n2, v);
- }
- }
- /* no neighbors */
- return 0.;
-}
-
-static gdouble center_regular_2nd_derivative (FttCell * cell, FttComponent c, GfsVariable * v)
-{
- FttCell * n1 = ftt_cell_neighbor (cell, 2*c);
- FttCell * n2 = ftt_cell_neighbor (cell, 2*c + 1);
- if (n1 && n2) {
- guint level = ftt_cell_level (cell);
- if (ftt_cell_level (n1) < level || ftt_cell_level (n2) < level)
- return center_regular_2nd_derivative (ftt_cell_parent (cell), c, v)/4.;
- return GFS_VALUE (n1, v) - 2.*GFS_VALUE (cell, v) + GFS_VALUE (n2, v);
- }
- /* one or no neighbors */
- return 0.;
-}
-
static void compute_gradient (FttCell * cell, GfsAdaptError * a)
{
- GFS_VALUE (cell, a->dv) = center_regular_gradient (cell, a->dv->component,
- GFS_ADAPT_GRADIENT (a)->v);
+ GFS_VALUE (cell, a->dv) = gfs_center_regular_gradient (cell, a->dv->component,
+ GFS_ADAPT_GRADIENT (a)->v);
}
static void add_hessian_norm (FttCell * cell, GfsAdaptError * a)
@@ -616,14 +572,19 @@ static void add_hessian_norm (FttCell * cell, GfsAdaptError * a)
FttComponent j;
for (j = 0; j < FTT_DIMENSION; j++)
if (j != a->dv->component) {
- gdouble g = center_regular_gradient (cell, j, a->dv);
+ gdouble g = gfs_center_regular_gradient (cell, j, a->dv);
GFS_VALUE (cell, a->v) += g*g;
}
/* diagonal */
- gdouble g = center_regular_2nd_derivative (cell, a->dv->component, GFS_ADAPT_GRADIENT (a)->v);
+ gdouble g = gfs_center_regular_2nd_derivative (cell, a->dv->component, GFS_ADAPT_GRADIENT (a)->v);
GFS_VALUE (cell, a->v) += g*g;
}
+static void scale (FttCell * cell, GfsAdaptError * a)
+{
+ GFS_VALUE (cell, a->v) = sqrt (GFS_VALUE (cell, a->v))/8.*GFS_ADAPT_GRADIENT (a)->dimension;
+}
+
static gboolean gfs_adapt_error_event (GfsEvent * event,
GfsSimulation * sim)
{
@@ -632,19 +593,18 @@ static gboolean gfs_adapt_error_event (GfsEvent * event,
GfsAdaptError * a = GFS_ADAPT_ERROR (event);
GfsDomain * domain = GFS_DOMAIN (sim);
- gfs_domain_cell_traverse (domain,
- FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) gfs_cell_reset, a->v);
a->dv = gfs_temporary_variable (domain);
for (a->dv->component = 0; a->dv->component < FTT_DIMENSION; a->dv->component++) {
- gfs_domain_cell_traverse (domain,
- FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) compute_gradient, a);
gfs_domain_bc (domain, FTT_TRAVERSE_ALL, -1, a->dv);
- gfs_domain_cell_traverse (domain,
- FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) add_hessian_norm, a);
}
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) scale, a);
gts_object_destroy (GTS_OBJECT (a->dv));
return TRUE;
}
@@ -660,7 +620,7 @@ static void gfs_adapt_error_class_init (GfsEventClass * klass)
static gdouble cost_error (FttCell * cell, GfsAdaptError * a)
{
- return sqrt (GFS_VALUE (cell, a->v))/8.*GFS_ADAPT_GRADIENT (a)->dimension;
+ return GFS_VALUE (cell, a->v);
}
static void gfs_adapt_error_init (GfsAdapt * object)
@@ -997,8 +957,8 @@ typedef struct {
gboolean changed;
} AdaptLocalParams;
-#define REFINABLE(cell, p) (GFS_VARIABLE (cell, (p)->r->i))
-#define COARSENABLE(cell, p) (GFS_VARIABLE (cell, (p)->c->i))
+#define REFINABLE(cell, p) (GFS_VALUE (cell, (p)->r))
+#define COARSENABLE(cell, p) (GFS_VALUE (cell, (p)->c))
static gboolean coarsen_cell (FttCell * cell, AdaptLocalParams * p)
{
diff --git a/src/fluid.c b/src/fluid.c
index 9186c3a..04153f5 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -592,6 +592,87 @@ gdouble gfs_center_minmod_gradient (FttCell * cell,
}
/**
+ * gfs_center_regular_gradient:
+ * @cell: a #FttCell.
+ * @c: a component.
+ * @v: a #GfsVariable.
+ *
+ * The gradient is normalized by the size of the cell. Only regular
+ * Cartesian stencils are used to compute the gradient.
+ *
+ * Returns: the value of the @c component of the gradient of variable @v
+ * at the center of the cell.
+ */
+gdouble gfs_center_regular_gradient (FttCell * cell,
+ FttComponent c,
+ GfsVariable * v)
+{
+ g_return_val_if_fail (cell != NULL, 0.);
+ g_return_val_if_fail (c < FTT_DIMENSION, 0.);
+ g_return_val_if_fail (v != NULL, 0.);
+
+ FttCell * n1 = ftt_cell_neighbor (cell, 2*c);
+ guint level = ftt_cell_level (cell);
+ if (n1) {
+ if (ftt_cell_level (n1) < level)
+ return gfs_center_regular_gradient (ftt_cell_parent (cell), c, v)/2.;
+ FttCell * n2 = ftt_cell_neighbor (cell, 2*c + 1);
+ if (n2) {
+ if (ftt_cell_level (n2) < level)
+ return gfs_center_regular_gradient (ftt_cell_parent (cell), c, v)/2.;
+ /* two neighbors: second-order differencing (parabola) */
+ return (GFS_VALUE (n1, v) - GFS_VALUE (n2, v))/2.;
+ }
+ else
+ /* one neighbor: first-order differencing */
+ return GFS_VALUE (n1, v) - GFS_VALUE (cell, v);
+ }
+ else {
+ FttCell * n2 = ftt_cell_neighbor (cell, 2*c + 1);
+ if (n2) {
+ if (ftt_cell_level (n2) < level)
+ return gfs_center_regular_gradient (ftt_cell_parent (cell), c, v)/2.;
+ /* one neighbor: first-order differencing */
+ return GFS_VALUE (cell, v) - GFS_VALUE (n2, v);
+ }
+ }
+ /* no neighbors */
+ return 0.;
+}
+
+/**
+ * gfs_center_regular_2nd_derivative:
+ * @cell: a #FttCell.
+ * @c: a component.
+ * @v: a #GfsVariable.
+ *
+ * The derivative is normalized by the size of the cell squared. Only
+ * regular Cartesian stencils are used to compute the derivative.
+ *
+ * Returns: the value of the @c component of the 2nd derivative of
+ * variable @v at the center of the cell.
+ */
+gdouble gfs_center_regular_2nd_derivative (FttCell * cell,
+ FttComponent c,
+ GfsVariable * v)
+{
+ g_return_val_if_fail (cell != NULL, 0.);
+ g_return_val_if_fail (c < FTT_DIMENSION, 0.);
+ g_return_val_if_fail (v != NULL, 0.);
+
+ FttCell * n1 = ftt_cell_neighbor (cell, 2*c);
+ FttCell * n2 = ftt_cell_neighbor (cell, 2*c + 1);
+ if (n1 && n2) {
+ guint level = ftt_cell_level (cell);
+ if (ftt_cell_level (n1) < level || ftt_cell_level (n2) < level)
+ return gfs_center_regular_2nd_derivative (ftt_cell_parent (cell), c, v)/4.;
+ return GFS_VALUE (n1, v) - 2.*GFS_VALUE (cell, v) + GFS_VALUE (n2, v);
+ }
+ /* one or no neighbors */
+ return 0.;
+}
+
+/**
* gfs_face_gradient:
* @face: a #FttCellFace.
* @g: the #GfsGradient.
diff --git a/src/fluid.h b/src/fluid.h
index 6af8975..c77990c 100644
--- a/src/fluid.h
+++ b/src/fluid.h
@@ -121,6 +121,12 @@ gdouble gfs_center_van_leer_gradient (FttCell * cell,
gdouble gfs_center_minmod_gradient (FttCell * cell,
FttComponent c,
guint v);
+gdouble gfs_center_regular_gradient (FttCell * cell,
+ FttComponent c,
+ GfsVariable * v);
+gdouble gfs_center_regular_2nd_derivative (FttCell * cell,
+ FttComponent c,
+ GfsVariable * v);
typedef struct _GfsGradient GfsGradient;
diff --git a/src/wave.c b/src/wave.c
index e0afa3f..9bba35a 100644
--- a/src/wave.c
+++ b/src/wave.c
@@ -35,7 +35,7 @@ static double theta (guint ith, guint ntheta)
return 2.*M_PI*ith/ntheta;
}
-static void cg (int ik, int ith, FttVector * u, guint ntheta, gdouble g)
+static void group_velocity (int ik, int ith, FttVector * u, guint ntheta, gdouble g)
{
double cg = g/(4.*M_PI*frequency (ik));
u->x = cg*cos (theta (ith, ntheta));
@@ -83,6 +83,76 @@ static void solid_flux (FttCell * cell, SolidFluxParams * par)
GFS_VALUE (cell, par->fv) = 0.;
}
+typedef struct {
+ GfsVariable * F, * Fn, * dF;
+ gdouble D[2][2];
+} GSEData;
+
+static void compute_gradient (FttCell * cell, GSEData * p)
+{
+ GFS_VALUE (cell, p->dF) = gfs_center_regular_gradient (cell, p->dF->component, p->Fn);
+}
+
+static void diffusion (FttCell * cell, GSEData * p)
+{
+ gdouble h2 = ftt_cell_size (cell);
+ h2 *= h2;
+ /* off-diagonal */
+ FttComponent j;
+ for (j = 0; j < 2; j++)
+ if (j != p->dF->component)
+ GFS_VALUE (cell, p->F) +=
+ p->D[j][p->dF->component]*gfs_center_regular_gradient (cell, j, p->dF)/h2;
+ /* diagonal */
+ GFS_VALUE (cell, p->F) +=
+ p->D[p->dF->component][p->dF->component]*
+ gfs_center_regular_2nd_derivative (cell, p->dF->component, p->Fn)/h2;
+}
+
+static void copy_F (FttCell * cell, GSEData * p)
+{
+ GFS_VALUE (cell, p->Fn) = GFS_VALUE (cell, p->F);
+}
+
+static void gse_alleviation_diffusion (GfsDomain * domain, GfsVariable * F,
+ FttVector * cg, gdouble dt)
+{
+ gfs_domain_timer_start (domain, "gse_alleviation");
+
+ gdouble ncg = sqrt (cg->x*cg->x + cg->y*cg->y);
+ gdouble dcg = (GFS_WAVE_GAMMA - 1./GFS_WAVE_GAMMA)*ncg/2.;
+ gdouble Ts = 4.*GFS_WAVE (domain)->alpha_s*GFS_WAVE (domain)->alpha_s*dt;
+ gdouble Dss = dt*dcg*dcg*Ts/12.;
+ gdouble dtheta = 2.*M_PI/GFS_WAVE (domain)->ntheta;
+ gdouble Dnn = dt*(ncg*dtheta)*(ncg*dtheta)*Ts/12.;
+ GSEData p;
+ gdouble cost = cg->x/ncg, sint = cg->y/ncg;
+ p.D[0][0] = Dss*cost*cost + Dnn*sint*sint;
+ p.D[1][1] = Dss*sint*sint + Dnn*cost*cost;
+ p.D[0][1] = p.D[1][0] = (Dss - Dnn)*cost*sint;
+ p.F = F;
+ p.Fn = gfs_temporary_variable (domain);
+ p.dF = gfs_temporary_variable (domain);
+ gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) copy_F, &p);
+ gfs_domain_cell_traverse (domain, FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
+ (FttCellTraverseFunc) p.Fn->fine_coarse, p.Fn);
+ for (p.dF->component = 0; p.dF->component < 2; p.dF->component++) {
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) compute_gradient, &p);
+ gfs_domain_bc (domain, FTT_TRAVERSE_ALL, -1, p.dF);
+ gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) diffusion, &p);
+ }
+ gts_object_destroy (GTS_OBJECT (p.Fn));
+ gts_object_destroy (GTS_OBJECT (p.dF));
+ gfs_domain_timer_stop (domain, "gse_alleviation");
+}
+
+static void redo_some_events (GfsEvent * event, GfsSimulation * sim)
+{
+ if (GFS_IS_ADAPT (event) || GFS_IS_INIT (event))
+ gfs_event_redo (event, sim);
+}
+
static void wave_run (GfsSimulation * sim)
{
GfsDomain * domain = GFS_DOMAIN (sim);
@@ -115,36 +185,45 @@ static void wave_run (GfsSimulation * sim)
/* spatial advection */
guint ik, ith;
for (ik = 0; ik < wave->nk; ik++) {
- FttVector u;
- cg (ik, 0, &u, wave->ntheta, g);
+ FttVector cg;
+ group_velocity (ik, 0, &cg, wave->ntheta, g);
gfs_domain_face_traverse (domain, FTT_XYZ,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttFaceTraverseFunc) set_group_velocity, &u);
- gfs_simulation_set_timestep (sim);
+ (FttFaceTraverseFunc) set_group_velocity, &cg);
+ if (wave->alpha_s > 0.) {
+ /* stability criterion for GSE diffusion */
+ gdouble cfl = sim->advection_params.cfl;
+ sim->advection_params.cfl = MIN (cfl, 2./(4.*wave->alpha_s*M_PI/wave->ntheta));
+ gfs_simulation_set_timestep (sim);
+ sim->advection_params.cfl = cfl;
+ }
+ else
+ gfs_simulation_set_timestep (sim);
/* subcycling */
guint n = rint (dt/sim->advection_params.dt);
g_assert (fabs (sim->time.t + sim->advection_params.dt*n - tnext) < 1e-12);
while (n--) {
for (ith = 0; ith < wave->ntheta; ith++) {
- FttVector u;
- cg (ik, ith, &u, wave->ntheta, g);
+ FttVector cg;
+ group_velocity (ik, ith, &cg, wave->ntheta, g);
gfs_domain_face_traverse (domain, FTT_XYZ,
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttFaceTraverseFunc) set_group_velocity, &u);
+ (FttFaceTraverseFunc) set_group_velocity, &cg);
GfsVariable * t = GFS_WAVE (sim)->F[ik][ith];
sim->advection_params.v = t;
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) solid_flux, &par);
+ gfs_domain_traverse_leaves (domain, (FttCellTraverseFunc) solid_flux, &par);
gfs_tracer_advection_diffusion (domain, &sim->advection_params);
sim->advection_params.fv = par.fv;
gfs_domain_traverse_merged (domain, (GfsMergedTraverseFunc) gfs_advection_update,
&sim->advection_params);
+ if (wave->alpha_s > 0.)
+ gse_alleviation_diffusion (domain, t, &cg, sim->advection_params.dt);
gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, t);
gfs_domain_cell_traverse (domain,
FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
(FttCellTraverseFunc) t->fine_coarse, t);
}
- gts_container_foreach (GTS_CONTAINER (sim->adapts), (GtsFunc) gfs_event_redo, sim);
+ gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) redo_some_events, sim);
gfs_simulation_adapt (sim);
}
}
@@ -184,12 +263,11 @@ static void wave_read (GtsObject ** o, GtsFile * fp)
GfsWave * wave = GFS_WAVE (*o);
if (fp->type == '{') {
GtsFileVariable var[] = {
- {GTS_UINT, "nk", TRUE},
- {GTS_UINT, "ntheta", TRUE},
+ {GTS_UINT, "nk", TRUE, &wave->nk},
+ {GTS_UINT, "ntheta", TRUE, &wave->ntheta},
+ {GTS_DOUBLE, "alpha_s", TRUE, &wave->alpha_s},
{GTS_NONE}
};
- var[0].data = &wave->nk;
- var[1].data = &wave->ntheta;
gts_file_assign_variables (fp, var);
if (fp->type == GTS_ERROR)
return;
@@ -218,8 +296,9 @@ static void wave_write (GtsObject * o, FILE * fp)
fprintf (fp, " {\n"
" nk = %d\n"
" ntheta = %d\n"
+ " alpha_s = %g\n"
"}",
- wave->nk, wave->ntheta);
+ wave->nk, wave->ntheta, wave->alpha_s);
}
static void gfs_wave_class_init (GfsSimulationClass * klass)
@@ -250,6 +329,7 @@ static void wave_init (GfsWave * wave)
{
wave->nk = 25;
wave->ntheta = 24;
+ wave->alpha_s = 0.;
/* default for g is acceleration of gravity on Earth with kilometres as
spatial units, hours as time units and Hz as frequency units */
GFS_SIMULATION (wave)->physical_params.g = 9.81/1000.*3600.;
diff --git a/src/wave.h b/src/wave.h
index fe685c9..545ebf7 100644
--- a/src/wave.h
+++ b/src/wave.h
@@ -41,6 +41,7 @@ struct _GfsWave {
/*< public >*/
guint nk, ntheta;
+ gdouble alpha_s;
GfsVariable *** F;
};
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list