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

revision 2099 by jfenwick, Wed Nov 5 04:59:22 2008 UTC revision 2100 by gross, Wed Nov 26 08:13:00 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
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]
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]
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}
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}
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

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

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

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
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{LIT DRIVEN CAVITY}:
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)
204          (whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1]
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$ such that
215    \label{DARCY PROBLEM}
216    \begin{array}{rcl}
217    u\hackscore{i} + \kappa\hackscore{ij} p\hackscore{,j} & = & g\hackscore{i} \\
218    u\hackscore{k,k} & = & f
219    \end{array}
220
221    with the boundary conditions
222    \label{DARCY BOUNDARY}
223    \begin{array}{rcl}
224    u\hackscore{i} \; n\hackscore{i}  = u^{N}\hackscore{i}  \; n\hackscore{i} & \mbox{ on } & \Gamma\hackscore{N} \\
225    p = p^{D} &  \mbox{ on } & \Gamma\hackscore{D} \\
226    \end{array}
227
228    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
229
230    \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}
231
232    for all $x\hackscore{i}$.
233
example usage:
234
235  solution=StokesProblemCartesian(mesh) \\  \subsection{Solution Method \label{DARCY SOLVE}}
236  solution.setTolerance(TOL) \\  Without loss of generality we can assume that $u^{N}\hackscore{i} \; n\hackscore{i}=0$ and
237  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
238  velocity,pressure=solution.solve(velocity,pressure,max\_iter=max\_iter,solver=solver) \\
239    \begin{array}{rcl}
240    g\hackscore{i} & \leftarrow & g\hackscore{i} - u^{N}\hackscore{i} -  \kappa\hackscore{ij} p^{D}\hackscore{,j }\\
241    f & \leftarrow & f - u^{N}\hackscore{k,k}
242    \end{array}
243
244    We set
245
246    V = \{ q \in H^{1}(\Omega) : q=0 \mbox{ on } \Gamma\hackscore{D} \}
247
248    and
249
250    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} \}
251
252    and define the operator $Q: V \rightarrow (L^2(\Omega))^{d}$ defined by
253
254    (Qp)\hackscore{i} = \kappa\hackscore{ij} p\hackscore{,j}
255
256    and the operator $D: W \rightarrow L^2(\Omega)$ defined by
257
258    Dv = v\hackscore{k,k}
259
260    In operator notation the Darcy problem~\ref{DARCY PROBLEM} is written in the form
261
262    \begin{array}{rcl}
263    u + Qp & = & g \\
264    Du & = & f
265    \end{array}
266
267    We solve this equation by minimising the functional
268
269    J(u,p):=\|u + Qp - g\|^2\hackscore{0} + \|Du-f\|\hackscore{0}^2
270
271    over $W \times V$ where $\|.\|\hackscore{0}$ denotes the norm in $L^2(\Omega)$. A simple calculation shows that
272    one has to solve
273
274    ( v + Qq , u + Qp - g) + (Dv,Du-f) =0
275
276    for all $v\in W$ and $q \in V$.which translates back into operator notation
277
278    \begin{array}{rcl}
279    (I+D^*D)u + Qp & = & D^*f + g \\
280    Q^*u  + Q^*Q p & = & Q^* g \\
281    \end{array}
282
283    where $D^*$ and $Q^*$ denote the adjoint operators.
284    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
285    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$)
286
287    The approach we are taking is to eliminate the velocity $u$ from the problem. Assuming that $p$ is known we have
288
289    v= (I+D^*D)^{-1} (D^*f + g - Qp)
290
291    (notice that $(I+D^*D)$ is coercive in $W$) which is inserted into the second equation
292
293    Q^* (I+D^*D)^{-1} (D^*f + g - Qp) + Q^* Q p = Q^* g
294
295    which is
296
297    Q^* ( I - (I+D^*D)^{-1} ) Q p = Q^* ( g -(I+D^*D)^{-1} (D^*f + g) )
298
299    We use the PCG method to solve this. The residual $r$ ($\in V^*$) is given as
300
301    \begin{array}{rcl}
302    r & = & Q^* ( g -(I+D^*D)^{-1} (D^*f + g) - Qp + (I+D^*D)^{-1}Q p \\
303    & =&  Q^* \left( (g-Qp) - (I+D^*D)^{-1} (D^*f + g - Qp) \right) \\
304    & =&  Q^* \left( (g-Qp) - v \right)
305    \end{array}
306
307    So in a partical implementation we use the pair $(g-Qp,v)$ to represent the residual. This will save the
308    reconstruction of the velocity $v$. In this notation the right hand side is given as
309    $(g,(I+D^*D)^{-1} (D^*f + g))$. The evaluation of the iteration operator for a given $p$ is then
310    returning $(Qp,w)$ where $w$ is the solution of
311    \label{UPDATE W}
312    (I+D^*D)w = Qp
313
314    We use $Q^*Q$ as a a preconditioner for the iteration operator $Q^* ( I - (I+D^*D)^{-1} ) Q$.
315
316  % \subsection{Benchmark Problem}  \subsection{Functions}
%
% Convection problem
317
318    \subsection{Example: Gravity Flow}
319
320  \section{Temperature Cartesian}  %================================================
322
323
324  \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 335  where $\vec{v}$ is the velocity vector,
335  \end{classdesc}  \end{classdesc}
336
337  \subsection{Benchmark Problem}  \subsection{Benchmark Problem}
338    %===============================================================================================================
339
340    %=========================================================
341    % \section{Level Set Method}
342
343  \section{Level Set Method}  %\subsection{Description}
344
345  \subsection{Description}  %\subsection{Method}

\subsection{Method}
346
348
349  \begin{classdesc}{LevelSet}{mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth}  %\begin{classdesc}{LevelSet}{mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth}
350  \end{classdesc}  %\end{classdesc}
351
352  %example usage:  %example usage:
353
354  %levelset = LevelSet(mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth)  %levelset = LevelSet(mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth)
355
356  \begin{methoddesc}[LevelSet]{update\_parameter}{parameter}  %\begin{methoddesc}[LevelSet]{update\_parameter}{parameter}
357  Update the parameter.  %Update the parameter.
358  \end{methoddesc}  %\end{methoddesc}
359
360  \begin{methoddesc}[LevelSet]{update\_phi}{paramter}{velocity}{dt}{t\_step}  %\begin{methoddesc}[LevelSet]{update\_phi}{paramter}{velocity}{dt}{t\_step}
361  Update level set function; advection and reinitialization  %Update level set function; advection and reinitialization
362  \end{methoddesc}  %\end{methoddesc}
363
364  \subsection{Benchmark Problem}  %\subsection{Benchmark Problem}
365
366  Rayleigh-Taylor instability problem  %Rayleigh-Taylor instability problem
367
368
369  % \section{Drucker Prager Model}  % \section{Drucker Prager Model}
# Line 191  for a given power law coefficients $n^{q Line 397 for a given power law coefficients$n^{q
397  Notice that $n^{q}=1$ gives a constant viscosity.  Notice that $n^{q}=1$ gives a constant viscosity.
398  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:
399  \label{IKM-EQU-6}  \label{IKM-EQU-6}
400  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}} \;.
401
402  With  With
403  \label{IKM-EQU-8}  \label{IKM-EQU-8}
# Line 239  D\hackscore{ij}'^{el}=\frac{1}{2 \mu dt Line 445  D\hackscore{ij}'^{el}=\frac{1}{2 \mu dt
445  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.
446  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
447  \label{IKM-EQU-10}  \label{IKM-EQU-10}
448  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff} D\hackscore{ij}' +  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff}  \left( D\hackscore{ij}' +
449  \frac{\eta\hackscore{eff}}{\mu \; dt}  \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)  \mbox{ with }
\sigma\hackscore{ij}^{'-} \mbox{ with }
450  \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}}
451
452  Notice that $\eta\hackscore{eff}$ is a function of stress.  Notice that $\eta\hackscore{eff}$ is a function of diatoric stress $\sigma\hackscore{ij}'$.
453  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
454  \label{IKM-EQU-1ib}  \label{IKM-EQU-1ib}
455  -\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})
456    \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+
457  \frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij,j}^{'-}  \frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij,j}^{'-}
458
459  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
460  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
461
462  \begin{array}{rcl}  \begin{array}{rcl}
463  -\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}
464  \frac{1}{\eta\hackscore{eff}} dp & = & - v\hackscore{i,i}^+ \\  )\right)\hackscore{,j} & = & F\hackscore{i}+ \sigma\hackscore{ij,j}'-p\hackscore{,i} \\
465  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}^+
466  \end{array}  \end{array}
467  \label{IKM iteration 2}  \label{IKM iteration 2}
468
469  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
470
471  \|(\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
472  \label{IKM iteration 3}  \label{IKM iteration 3}
473
474    where  $\bar{\eta}\hackscore{eff}$ is the caracteristic viscosity, for instance:
475
476    \frac{1}{\bar{\eta}\hackscore{eff}} = \frac{1}{\tau^{-}}+\sum\hackscore{q}  \frac{1}{\eta^{q}\hackscore{N}}
477    \label{IKM iteration 4}
478
479    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
480    calculation of $\eta\hackscore{eff}$ with the Newton-Raphson scheme. This value can then be used to calculate
481    $\sigma\hackscore{ij}'$ via~\ref{IKM-EQU-10}. We need to solve
482
483    \tau = \eta\hackscore{eff} \cdot \epsilon \mbox{ with }
484    \epsilon = \sqrt{ 2 \left( D\hackscore{ij}' +
485    \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)^2}
486    \label{IKM iteration 5}
487
488    The Newton scheme takes the form
489
490    \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)
491    = \min(\frac{\eta\hackscore{eff} - \tau\hackscore{n}  \eta\hackscore{eff}'}
492    {1 - \eta\hackscore{eff}'  \cdot  \epsilon}, \frac{\tau\hackscore{Y} + \beta \; p}{\epsilon}) \epsilon
493    \label{IKM iteration 6}
494
495    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
496
497    \eta\hackscore{eff}' = - \eta\hackscore{eff}^2 \left(\frac{1}{\eta\hackscore{eff}}\right)'
498    \mbox{ with }
499    \left(\frac{1}{\eta\hackscore{eff}}\right)' = \sum\hackscore{q} \left(\frac{1}{\eta^{q}} \right)'
500    \label{IKM iteration 7}
501
502    \label{IKM-EQU-5XX}
503    \left(\frac{1}{\eta^{q}} \right)'
504    = \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}}}}
505    = \frac{1-\frac{1}{n^{q}}}{ \tau \eta^{q}}
506
507    Notice that allways $\eta\hackscore{eff}'\le 0$ which makes the denomionator in~\ref{IKM iteration 6}
508    positive.
509
510
511

Legend:
 Removed from v.2099 changed lines Added in v.2100