# Diff of /trunk/doc/user/Models.tex

revision 1966 by jfenwick, Wed Nov 5 04:59:22 2008 UTC revision 2108 by gross, Fri Nov 28 05:09:23 2008 UTC
# Line 17  Line 17
17  The following sections give a breif overview of the model classes and their corresponding methods.  The following sections give a breif overview of the model classes and their corresponding methods.
18
19  \section{Stokes Problem}  \section{Stokes Problem}
20  The velocity field $v$ and pressure $p$ of an incompressible fluid is given as the solution of their  The velocity \index{velocity} field $v$ and pressure $p$ of an incompressible fluid \index{incompressible fluid} is given as the solution of the Stokes problem\index{Stokes problem}
Stokes problem
21  \label{Stokes 1}  \label{Stokes 1}
22  -\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}  -\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=f\hackscore{i}-\sigma\hackscore{ij,j}
23
24  where $\eta$ is the viscosity and $F_i$ defines an internal force. We assume an incompressible media:  where $\eta$ is the viscosity, $F\hackscore{i}$ defines an internal force \index{force, internal} and $\sigma\hackscore{ij}$ is an intial stress \index{stress, initial}. We assume an incompressible media:
25  \label{Stokes 2}  \label{Stokes 2}
26  -v\hackscore{i,i}=0  -v\hackscore{i,i}=0
27
28  Natural boundary conditions are taken in the form  Natural boundary conditions are taken in the form
29  \label{Stokes Boundary}  \label{Stokes Boundary}
30  \left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)n\hackscore{j}-n\hackscore{i}p=f  \left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)n\hackscore{j}-n\hackscore{i}p=s\hackscore{i}+\sigma\hackscore{ij} n\hackscore{i}
31
32  which can be overwritten by a constraint  which can be overwritten by constraints of the form
33  \label{Stokes Boundary0}  \label{Stokes Boundary0}
34  v\hackscore{i}(x)=0  v\hackscore{i}(x)=v^D\hackscore{i}(x)
35
36  where the index $i$ may depend on the location $x$ on the bondary.  at some locations $x$ at the boundary of the domain. The index $i$ may depend on the location $x$ on the boundary.
37    $v^D$ is a given function on the domain.
38
39  \subsection{Solution Method \label{STOKES SOLVE}}  \subsection{Solution Method \label{STOKES SOLVE}}
40  In block form equation equations~\ref{Stokes 1} and~\ref{Stokes 2} takes the form of a saddle point problem  In block form equation equations~\ref{Stokes 1} and~\ref{Stokes 2} takes the form of a saddle point problem
41    \index{saddle point problem}
42
43  \left[ \begin{array}{cc}  \left[ \begin{array}{cc}
44  A     & B^{*} \\  A     & B^{*} \\
45  B & 0 \\  B & 0 \\
46  \end{array} \right]  \end{array} \right]
47  \left[ \begin{array}{c}  \left[ \begin{array}{c}
48  u \\  v \\
49  p \\  p \\
50  \end{array} \right]  \end{array} \right]
51  =\left[ \begin{array}{c}  =\left[ \begin{array}{c}
52  F \\  G \\
53  0 \\  0 \\
54  \end{array} \right]  \end{array} \right]
55  \label{SADDLEPOINT}  \label{SADDLEPOINT}
56
57  where $A$ is coercive, self-adjoint linear operator in $V$, $B$ is a divergence operator and $B^{*}$ is it adjoint operator (=gradient operator)). For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}.  where $A$ is coercive, self-adjoint linear operator in a suitable Hilbert space, $B$ is the $(-1) \cdot$ divergence operator and $B^{*}$ is it adjoint operator (=gradient operator). For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}.
58  We apply the preconditioner matrix  We use iterative techniques to solve this problem. To make sure that the incomressibilty condition holds
59    with sufficient accuracy we check for
60
61    \|v\hackscore{k,k}\| \hackscore \le  \epsilon
62    \|\sqrt{v\hackscore{j,k}v\hackscore{j,k}}\|
63
64    where $\epsilon$ is the desired relative accuracy and
65
66    \|p\|^2= \int\hackscore{\Omega} p^2 \; dx
67    \label{PRESSURE NORM}
68
69    defines the $L^2$-norm.
70    There are two approaches to solve this problem. The first approach, called the Uzawa scheme \index{Uzawa scheme}
71    eliminates the velocity $v$ from the problem. The second approach solves the equation in coupled form after the application of a preconditioner.
72
73    \subsubsection{Uzawa scheme}
74    The first eqution in~\ref{SADDLEPOINT} gives $v=A^{-1}(G-B^{*}p)$ assuming $p$ is known. This is inserted into the
75    second eqution which leads to
76
77    S p =  B A^{-1} G
78
79    with the Schur complement \index{Schur complement} $S=BA^{-1}B^{*}$. This problem can be solved iteratively using the reconditioned Conjugate Gradient Method (PCG)~\index{PCG!Preconditioned Conjugate Gradient Method}
80    with the preconditioner $\hat{S}$ defined as $q=\hat{S}^{-1}p$ by solving
81
82    \frac{1}{\eta}q = p
83
84    see~\cite{ELMAN} for more details. The evaluation of $w=Sp$ is done in the form
85
86    \begin{array}{rcl}
87    A v & = & B^{*}p \\
88    w & = & Bv \\
89    \end{array}
90    \label{EVAL PCG}
91
92    The residual \index{residual}  $r=B A^{-1} G - S p$ is given as
93
94    r=B A^{-1} (G - B^* p) = Bv \mbox{ with } v = A^{-1}(G-B^{*}p)
95
96    Therefore one uses the tuple $(v,Bv)$ to represent the residual of the current pressure $p$. Notice that before the iteration is started the right hand side $B A^{-1} G$ needs to be calculated. The bilinear form $(.,.)$ used is defined as
97
98    (p,(v,Bv))=\int\hackscore{\Omega} p \cdot Bv \; dx
99
100    where $p$ is the pressure increment and $(v,Bv)$ represents an increment in the residual.
101
102    \subsubsection{Coupled Solver}
103    An alternative approach to solve the saddle point problem~\ref{SADDLEPOINT} directly using an iterative such as
104    the generalized minimal residual method (GMRES) \index{generalized minimal residual method!GMRES} with a suitable
105    preconditioner. Here we use the operator
106
107  \left[ \begin{array}{cc}  \left[ \begin{array}{cc}
108  A^{-1}     & 0 \\  A^{-1}     & 0 \\
# Line 63  S^{-1} B A^{-1}  & -S^{-1} \\ Line 110  S^{-1} B A^{-1}  & -S^{-1} \\
110  \end{array} \right]  \end{array} \right]
111  \label{SADDLEPOINT PRECODITIONER}  \label{SADDLEPOINT PRECODITIONER}
112
113  with the Schur complement $S=BA^{-1}B^{*}$ to solve the problem iteratively. The updates $dv$ and $dp$ for  where again $S$ is the Schur complement~\cite{ELMAN}. In partice we will use an approximation $\hat{S}$ for $S$. The evaluation $(w,q)$ of the iteration operator for a given $(v,p)$ is done as
given velocity $v$ and pressure $p$ is given as
114
115  \begin{array}{rcl}  \begin{array}{rcl}
116  A dv & = & F-Av-B^{*}p \\  A w & = & Av+B^{*}p \\
117  S dp & = & B(v+dv) \\  \hat{S} q & = & B(w-v) \\
118  \end{array}  \end{array}
119  \label{SADDLEPOINT iteration}  \label{COUPLES SADDLEPOINT iteration}
120
121  This scheme is called the Uzawa scheme.  We use the inner product induced by the norm
122
123  In the case of the Stokes problem it can be shown that $S$ can be approximated by $\frac{1}{\eta}$. With this the iteration scheme can be implemented as  \|(v,p)\|^2= \int\hackscore{\Omega}  v\hackscore{i,j}  v\hackscore{i,j} + \left( \frac{p}{\eta}\right)^2\; dx
124    \label{COUPLES NORM}
125
126    In PDE form~\ref{COUPLES SADDLEPOINT iteration} takes the form
127
128  \begin{array}{rcl}  \begin{array}{rcl}
129  -\left(\eta(dv\hackscore{i,j}+ dv\hackscore{i,j})\right)\hackscore{,j} & = & F\hackscore{i}+\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}-p\hackscore{,i} \\  -\left(\eta(w\hackscore{i,j}+ w\hackscore{i,j})\right)\hackscore{,j} & = & -\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i} \\
130  \frac{1}{\eta} dp & = & - \left(v\hackscore{,i}+dv\hackscore{,i} \right)\hackscore{,i} \\  \frac{1}{\eta}  q & = & - (w-v)\hackscore{i,i} \\
131  \end{array}  \end{array}
132  \label{SADDLEPOINT iteration 2}  \label{SADDLEPOINT iteration 2}
133
To accelerate the convergence we are using the restarted $GMRES$ method using the norm

\|(v,p)\|^2= \int\hackscore{\Omega} \eta^2 v\hackscore{i,j}^2 + p^2 \; dx
\label{SADDLEPOINT iteration 3}

or alternatively the $PCG$ method on the pressure only using the norm

\|p\|^2= \int\hackscore{\Omega} \frac{1}{\eta^2} p^2 \; dx
\label{SADDLEPOINT iteration 4}

134
135
136    \subsection{Functions}
137
138    \begin{classdesc}{StokesProblemCartesian}{domain}
139    opens the Stokes problem\index{Stokes problem} on the \Domain domain. The approximation
140    order needs to be two.
141    \end{classdesc}
142
143    \begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}}}}}}
144    assigns values to the model parameters. In any call all values must be set.
145    \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,
146    \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.
147    The locations and compontents where the velocity is fixed are set by
148    the values of \var{fixed_u_mask}. The method will try to cast the given values to appropriate
149    \Data class objects.
150    \end{methoddesc}
151
152  \begin{classdesc}{StokesProblemCartesian}{domain,debug}  \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p,
153  opens the stokes equations on the \Domain domain. Setting debug=True switches the debug mode to on.  \optional{max_iter=20, \optional{verbose=False, \optional{useUzawa=True}}}}
154  \end{classdesc}  solves the problem and return approximations for velocity and pressure.
155    The arguments \var{v} and \var{p} define initial guess. The values of \var{v} marked
156    by \var{fixed_u_mask} remain unchanged.
157    If \var{useUzawa} is set to \True
158    the Uzawa\index{Uszwa} scheme is used. Otherwise the problem is solved in coupled form. In most cases
159    the Uzawa scheme is more efficient.
160    \var{max_iter} defines the maximum number of iteration steps.
161    If \var{verbose} is set to \True informations on the progress of of the solver are printed.
162    \end{methoddesc}
163
164
165    \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-8}}
166    sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1.
167    \end{methoddesc}
168    \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{}
169    returns the current relative tolerance.
170    \end{methoddesc}
171    \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}}
172    sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the
173    absolute talerance is set to 0.
174    \end{methoddesc}
175    \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{}
176    sreturns the current absolute tolerance.
177    \end{methoddesc}
178    \begin{methoddesc}[StokesProblemCartesian]{setSubToleranceReductionFactor}{\optional{reduction=None}}
179    sets the reduction factor for the tolerance used to solve the PDEs. A reduction factor
180    in the order of one will minimize compute time per iteration step but my slow down convergence or even lead to
181    divergency. On the other hand a very small value for the PDE tolerance could result in a wast of compute time.
182    If \var{reduction} is set to \var{None} the sub-tolerance is solved adaptively but
183    in cases a very small tolerance is set ($<10^{-6}$) it is recommended to set the
184    reduction factor by hand. This may require some experiments.
185    \end{methoddesc}
186    \begin{methoddesc}[StokesProblemCartesian]{getSubToleranceReductionFactor}{}
187    return the current reduction factor for the sub-problem tolerance.
188    \end{methoddesc}
189
190    \subsection{Example: Lit Driven Cavity}
191     The following script \file{lit\hackscore driven\hackscore cavity.py}
192    \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
193    illustrates the usage of the \class{StokesProblemCartesian} class to solve
194    the lit driven cavity problem~\cite{LITDRIVENCAVITY}:
195    \begin{python}
196    from esys.escript import *
197    from esys.finley import Rectangle
198    from esys.escript.models import StokesProblemCartesian
199    NE=25
200    dom = Rectangle(NE,NE,order=2)
201    x = dom.getX()
202    sc=StokesProblemCartesian(dom)
203    mask= (whereZero(x[0])*[1.,0]+whereZero(x[0]-1))*[1.,0] + \
204          (whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1]
205    sc.initialize(eta=.1, fixed_u_mask= mask)
206    v=Vector(0.,Solution(dom))
207    v[0]+=whereZero(x[1]-1.)
208    p=Scalar(0.,ReducedSolution(dom))
209    v,p=sc.solve(v,p, verbose=True)
210    saveVTK("u.xml",velocity=v,pressure=p)
211    \end{python}
212
213    \section{Darcy Flux}
214    We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$ solving
215    the Darcy flux problem \index{Darcy flux}\index{Darcy flow}
216    \label{DARCY PROBLEM}
217    \begin{array}{rcl}
218    u\hackscore{i} + \kappa\hackscore{ij} p\hackscore{,j} & = & g\hackscore{i} \\
219    u\hackscore{k,k} & = & f
220    \end{array}
221
222    with the boundary conditions
223    \label{DARCY BOUNDARY}
224    \begin{array}{rcl}
225    u\hackscore{i} \; n\hackscore{i}  = u^{N}\hackscore{i}  \; n\hackscore{i} & \mbox{ on } & \Gamma\hackscore{N} \\
226    p = p^{D} &  \mbox{ on } & \Gamma\hackscore{D} \\
227    \end{array}
228
229    where $\Gamma\hackscore{N}$ and $\Gamma\hackscore{D}$ are a partition of the boundary of $\Omega$ with $\Gamma\hackscore{D}$ non empty, $n\hackscore{i}$ is the outer normal field of the boundary of $\Omega$, $u^{N}\hackscore{i}$ and $p^{D}$ are given functions on $\Omega$, $g\hackscore{i}$ and $f$ are given source terms and $\kappa\hackscore{ij}$ is the given permability. We assume that $\kappa\hackscore{ij}$ is symmetric (which is not really required) and positive definite, i.e there are positive constants $\alpha\hackscore{0}$ and $\alpha\hackscore{1}$ wich are independent from the location in $\Omega$ such that
230
231    \alpha\hackscore{0} \; x\hackscore{i} x\hackscore{i} \le \kappa\hackscore{ij} x\hackscore{i} x\hackscore{j} \le \alpha\hackscore{1} \; x\hackscore{i} x\hackscore{i}
232
233    for all $x\hackscore{i}$.
234
example usage:
235
236  solution=StokesProblemCartesian(mesh) \\  \subsection{Solution Method \label{DARCY SOLVE}}
237  solution.setTolerance(TOL) \\  Without loss of generality we can assume that $u^{N}\hackscore{i} \; n\hackscore{i}=0$ and
238  solution.initialize(fixed\_u\_mask=b\_c,eta=eta,f=Y) \\  $p^{D}$. Otherewise one solves for $u-u^{N}$ and $p-p^{D}$ and sets
239  velocity,pressure=solution.solve(velocity,pressure,max\_iter=max\_iter,solver=solver) \\
240    \begin{array}{rcl}
241    g\hackscore{i} & \leftarrow & g\hackscore{i} - u^{N}\hackscore{i} -  \kappa\hackscore{ij} p^{D}\hackscore{,j }\\
242    f & \leftarrow & f - u^{N}\hackscore{k,k}
243    \end{array}
244
245    We set
246
247    V = \{ q \in H^{1}(\Omega) : q=0 \mbox{ on } \Gamma\hackscore{D} \}
248
249    and
250
251    W = \{ v \in (L^2(\Omega))^{d} : v\hackscore{k,k} \in L^2(\Omega) \mbox{ and } u\hackscore{i} \; n\hackscore{i} =0  \mbox{ on } \Gamma\hackscore{N} \}
252
253    and define the operator $Q: V \rightarrow (L^2(\Omega))^{d}$ defined by
254
255    (Qp)\hackscore{i} = \kappa\hackscore{ij} p\hackscore{,j}
256
257    and the operator $D: W \rightarrow L^2(\Omega)$ defined by
258
259    Dv = v\hackscore{k,k}
260
261    In operator notation the Darcy problem~\ref{DARCY PROBLEM} is written in the form
262
263    \begin{array}{rcl}
264    u + Qp & = & g \\
265    Du & = & f
266    \end{array}
267
268    We solve this equation by minimising the functional
269
270    J(u,p):=\|u + Qp - g\|^2\hackscore{0} + \|Du-f\|\hackscore{0}^2
271
272    over $W \times V$ where $\|.\|\hackscore{0}$ denotes the norm in $L^2(\Omega)$. A simple calculation shows that
273    one has to solve
274
275    ( v + Qq , u + Qp - g) + (Dv,Du-f) =0
276
277    for all $v\in W$ and $q \in V$.which translates back into operator notation
278
279    \begin{array}{rcl}
280    (I+D^*D)u + Qp & = & D^*f + g \\
281    Q^*u  + Q^*Q p & = & Q^* g \\
282    \end{array}
283
284    where $D^*$ and $Q^*$ denote the adjoint operators.
285    In~\cite{XXX} it has been shown that this problem is continuous and coercive in $W \times V$ and therefore has a unique solution. Also standart FEM methods can be used for discretization. It is also possible
286    to solve the problem is coupled form, however this approach leads in some cases to a very ill-conditioned stiffness matrix in particular in the case of a very small or large permability ($\alpha\hackscore{1} \ll 1$ or $\alpha\hackscore{0} \gg 1$)
287
288  % \subsection{Benchmark Problem}  The approach we are taking is to eliminate the velocity $u$ from the problem. Assuming that $p$ is known we have
289  %
290  % Convection problem  v= (I+D^*D)^{-1} (D^*f + g - Qp)
291
292    (notice that $(I+D^*D)$ is coercive in $W$) which is inserted into the second equation
293
294    Q^* (I+D^*D)^{-1} (D^*f + g - Qp) + Q^* Q p = Q^* g
295
296    which is
297
298    Q^* ( I - (I+D^*D)^{-1} ) Q p = Q^* ( g -(I+D^*D)^{-1} (D^*f + g) )
299
300    We use the PCG method to solve this. The residual $r$ ($\in V^*$) is given as
301
302    \begin{array}{rcl}
303    r & = & Q^* ( g -(I+D^*D)^{-1} (D^*f + g) - Qp + (I+D^*D)^{-1}Q p \\
304    & =&  Q^* \left( (g-Qp) - (I+D^*D)^{-1} (D^*f + g - Qp) \right) \\
305    & =&  Q^* \left( (g-Qp) - v \right)
306    \end{array}
307
308    So in a partical implementation we use the pair $(g-Qp,v)$ to represent the residual. This will save the
309    reconstruction of the velocity $v$. In this notation the right hand side is given as
310    $(g,(I+D^*D)^{-1} (D^*f + g))$. The evaluation of the iteration operator for a given $p$ is then
311    returning $(Qp,w)$ where $w$ is the solution of
312    \label{UPDATE W}
313    (I+D^*D)w = Qp
314
315    We use $Q^*Q$ as a a preconditioner for the iteration operator $Q^* ( I - (I+D^*D)^{-1} ) Q$.
316
317    \subsection{Functions}
318    \begin{classdesc}{DarcyFlow}{domain}
319    opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.
320    \end{classdesc}
321
322    \begin{methoddesc}[DarcyFlow]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}}}}}}
323    assigns values to the model parameters. In any call all values must be set.
324    \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,
325    \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.
326    The locations and compontents where the velocity is fixed are set by
327    the values of \var{fixed_u_mask}. The method will try to cast the given values to appropriate
328    \Data class objects.
329    \end{methoddesc}
330
331  \section{Temperature Cartesian}
332    \subsection{Example: Gravity Flow}
333
334    %================================================
335    \section{Temperature Advection Diffusion\label{TEMP ADV DIFF}}
336
337
338  \rho c\hackscore{p} \left (\frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T \right ) = k \nabla^{2}T  \rho c\hackscore{p} \left (\frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T \right ) = k \nabla^{2}T
# Line 130  where $\vec{v}$ is the velocity vector, Line 349  where $\vec{v}$ is the velocity vector,
349  \end{classdesc}  \end{classdesc}
350
351  \subsection{Benchmark Problem}  \subsection{Benchmark Problem}
352    %===============================================================================================================
353
354    %=========================================================
355    % \section{Level Set Method}
356
357  \section{Level Set Method}  %\subsection{Description}
358
359  \subsection{Description}  %\subsection{Method}

\subsection{Method}
360
361  Advection and Reinitialisation  %Advection and Reinitialisation
362
363  \begin{classdesc}{LevelSet}{mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth}  %\begin{classdesc}{LevelSet}{mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth}
364  \end{classdesc}  %\end{classdesc}
365
366  %example usage:  %example usage:
367
368  %levelset = LevelSet(mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth)  %levelset = LevelSet(mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth)
369
370  \begin{methoddesc}[LevelSet]{update\_parameter}{parameter}  %\begin{methoddesc}[LevelSet]{update\_parameter}{parameter}
371  Update the parameter.  %Update the parameter.
372  \end{methoddesc}  %\end{methoddesc}
373
374  \begin{methoddesc}[LevelSet]{update\_phi}{paramter}{velocity}{dt}{t\_step}  %\begin{methoddesc}[LevelSet]{update\_phi}{paramter}{velocity}{dt}{t\_step}
375  Update level set function; advection and reinitialization  %Update level set function; advection and reinitialization
376  \end{methoddesc}  %\end{methoddesc}
377
378  \subsection{Benchmark Problem}  %\subsection{Benchmark Problem}
379
380  Rayleigh-Taylor instability problem  %Rayleigh-Taylor instability problem
381
382
383  % \section{Drucker Prager Model}  % \section{Drucker Prager Model}
# Line 191  for a given power law coefficients $n^{q Line 411 for a given power law coefficients$n^{q
411  Notice that $n^{q}=1$ gives a constant viscosity.  Notice that $n^{q}=1$ gives a constant viscosity.
412  After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:  After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:
413  \label{IKM-EQU-6}  \label{IKM-EQU-6}
414  D\hackscore{ij}'^{vp}=\frac{1}{2 \eta^{vp}} \sigma'\hackscore{ij} \mbox{ with } \frac{1}{\eta^{vp}} \sum\hackscore{q} \frac{1}{\eta^{q}} \;.  D\hackscore{ij}'^{vp}=\frac{1}{2 \eta^{vp}} \sigma'\hackscore{ij} \mbox{ with } \frac{1}{\eta^{vp}} = \sum\hackscore{q} \frac{1}{\eta^{q}} \;.
415
416  With  With
417  \label{IKM-EQU-8}  \label{IKM-EQU-8}
# Line 239  D\hackscore{ij}'^{el}=\frac{1}{2 \mu dt Line 459  D\hackscore{ij}'^{el}=\frac{1}{2 \mu dt
459  where $\sigma\hackscore{ij}^{'-}$ is the deviatoric stress at the precious time step.  where $\sigma\hackscore{ij}^{'-}$ is the deviatoric stress at the precious time step.
460  Now we can combine equations~\ref{IKM-EQU-2}, \ref{IKM-EQU-3b} and~\ref{IKM-EQU-6b} to get  Now we can combine equations~\ref{IKM-EQU-2}, \ref{IKM-EQU-3b} and~\ref{IKM-EQU-6b} to get
461  \label{IKM-EQU-10}  \label{IKM-EQU-10}
462  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff} D\hackscore{ij}' +  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff}  \left( D\hackscore{ij}' +
463  \frac{\eta\hackscore{eff}}{\mu \; dt}  \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)  \mbox{ with }
\sigma\hackscore{ij}^{'-} \mbox{ with }
464  \frac{1}{\eta\hackscore{eff}}=\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}  \frac{1}{\eta\hackscore{eff}}=\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
465
466  Notice that $\eta\hackscore{eff}$ is a function of stress.  Notice that $\eta\hackscore{eff}$ is a function of diatoric stress $\sigma\hackscore{ij}'$.
467  After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get  After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get
468  \label{IKM-EQU-1ib}  \label{IKM-EQU-1ib}
469  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})
470    \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+
471  \frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij,j}^{'-}  \frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij,j}^{'-}
472
473  Together with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a problem with a form almost identical  Together with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a problem with a form almost identical
474  to the Stokes problem discussed in section~\ref{STOKES SOLVE} but with the difference that $\eta\hackscore{eff}$ is depending on the solution. Analog to the iteration scheme~\ref{SADDLEPOINT iteration 2} we can run  to the Stokes problem discussed in section~\ref{STOKES SOLVE} but with the difference that $\eta\hackscore{eff}$ is depending on the solution. Analog to the iteration scheme~\ref{SADDLEPOINT iteration 2} we can run
475
476  \begin{array}{rcl}  \begin{array}{rcl}
477  -\left(\eta\hackscore{eff}(dv\hackscore{i,j}+ dv\hackscore{i,j})\right)\hackscore{,j} & = & F\hackscore{i}+ \sigma\hackscore{ij,j}' \\  -\left(\eta\hackscore{eff}(dv\hackscore{i,j}+ dv\hackscore{i,j}
478  \frac{1}{\eta\hackscore{eff}} dp & = & - v\hackscore{i,i}^+ \\  )\right)\hackscore{,j} & = & F\hackscore{i}+ \sigma\hackscore{ij,j}'-p\hackscore{,i} \\
479  d\sigma\hackscore{ij}' & = & 2 \eta\hackscore{eff} D\hackscore{ij}^{+'} + \frac{\eta\hackscore{eff}}{\mu \; dt} \sigma\hackscore{ij}' -  \sigma\hackscore{ij}\\  \frac{1}{\eta\hackscore{eff}} dp & = & - v\hackscore{i,i}^+
480  \end{array}  \end{array}
481  \label{IKM iteration 2}  \label{IKM iteration 2}
482
483  where $v^+=v+dv$. As this problem is non-linear the Jacobi-free Newton-GMRES method is used with the norm  where $v^+=v+dv$. As this problem is non-linear the Jacobi-free Newton-GMRES method is used with the norm
484
485  \|(\sigma, p)\|^2= \int\hackscore{\Omega} \sigma\hackscore{i,j}^2 + p^2 \; dx  \|(v, p)\|^2= \int\hackscore{\Omega} v\hackscore{i,j}^2 + \frac{1}{\bar{\eta}^2\hackscore{eff}} p^2 \; dx
486  \label{IKM iteration 3}  \label{IKM iteration 3}
487
488    where  $\bar{\eta}\hackscore{eff}$ is the caracteristic viscosity, for instance:
489
490    \frac{1}{\bar{\eta}\hackscore{eff}} = \frac{1}{\tau^{-}}+\sum\hackscore{q}  \frac{1}{\eta^{q}\hackscore{N}}
491    \label{IKM iteration 4}
492
493    In oder to perform step~\ref{IKM iteration 2} we need to calculate the $\eta\hackscore{eff}$ as well as $\sigma\hackscore{ij}'$ while via $\tau$ the first is a function of the latter. The priority is the
494    calculation of $\eta\hackscore{eff}$ with the Newton-Raphson scheme. This value can then be used to calculate
495    $\sigma\hackscore{ij}'$ via~\ref{IKM-EQU-10}. We need to solve
496
497    \tau = \eta\hackscore{eff} \cdot \epsilon \mbox{ with }
498    \epsilon = \sqrt{ 2 \left( D\hackscore{ij}' +
499    \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)^2}
500    \label{IKM iteration 5}
501
502    The Newton scheme takes the form
503
504    \tau\hackscore{n+1} = \min(\tau\hackscore{n} - \frac{\tau\hackscore{n} - \eta\hackscore{eff}  \cdot \epsilon}{1 - \eta\hackscore{eff}'  \cdot  \epsilon}, \tau\hackscore{Y} + \beta \; p)
505    = \min(\frac{\eta\hackscore{eff} - \tau\hackscore{n}  \eta\hackscore{eff}'}
506    {1 - \eta\hackscore{eff}'  \cdot  \epsilon}, \frac{\tau\hackscore{Y} + \beta \; p}{\epsilon}) \epsilon
507    \label{IKM iteration 6}
508
509    where $\eta\hackscore{eff}'$ denotes the derivative of $\eta\hackscore{eff}$ with respect of $\tau$. The second term in $\min$ is droped of $\tau\hackscore{Y} + \beta \; p<0$ or $\epsilon=0$. In fact we have
510
511    \eta\hackscore{eff}' = - \eta\hackscore{eff}^2 \left(\frac{1}{\eta\hackscore{eff}}\right)'
512    \mbox{ with }
513    \left(\frac{1}{\eta\hackscore{eff}}\right)' = \sum\hackscore{q} \left(\frac{1}{\eta^{q}} \right)'
514    \label{IKM iteration 7}
515
516    \label{IKM-EQU-5XX}
517    \left(\frac{1}{\eta^{q}} \right)'
518    = \frac{1-\frac{1}{n^{q}}}{\eta^{q}\hackscore{N}} \frac{\tau^{-\frac{1}{n^{q}}}}{(\tau\hackscore{t}^q)^{1-\frac{1}{n^{q}}}}
519    = \frac{1-\frac{1}{n^{q}}}{ \tau \eta^{q}}
520
521    Notice that allways $\eta\hackscore{eff}'\le 0$ which makes the denomionator in~\ref{IKM iteration 6}
522    positive.
523
524
525

Legend:
 Removed from v.1966 changed lines Added in v.2108

 ViewVC Help Powered by ViewVC 1.1.26