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

revision 1966 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  \begin{equation}\label{Stokes 1}  \begin{equation}\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  \end{equation}  \end{equation}
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  \begin{equation}\label{Stokes 2}  \begin{equation}\label{Stokes 2}
26  -v\hackscore{i,i}=0  -v\hackscore{i,i}=0
27  \end{equation}  \end{equation}
28  Natural boundary conditions are taken in the form  Natural boundary conditions are taken in the form
29  \begin{equation}\label{Stokes Boundary}  \begin{equation}\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  \end{equation}  \end{equation}
32  which can be overwritten by a constraint  which can be overwritten by constraints of the form
33  \begin{equation}\label{Stokes Boundary0}  \begin{equation}\label{Stokes Boundary0}
34  v\hackscore{i}(x)=0  v\hackscore{i}(x)=v^D\hackscore{i}(x)
35  \end{equation}  \end{equation}
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  \begin{equation}  \begin{equation}
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  \end{equation}  \end{equation}
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    \begin{equation}
61    \|v\hackscore{k,k}\| \hackscore \le  \epsilon
62    \|\sqrt{v\hackscore{j,k}v\hackscore{j,k}}\|
63    \end{equation}
64    where $\epsilon$ is the desired relative accuracy and
65    \begin{equation}
66    \|p\|^2= \int\hackscore{\Omega} p^2 \; dx
67    \label{PRESSURE NORM}
68    \end{equation}
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    \begin{equation}
77    S p =  B A^{-1} G
78    \end{equation}
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    \begin{equation}
82    \frac{1}{\eta}q = p
83    \end{equation}
84    see~\cite{ELMAN} for more details. The evaluation of $w=Sp$ is done in the form
85    \begin{equation}
86    \begin{array}{rcl}
87    A v & = & B^{*}p \\
88    w & = & Bv \\
89    \end{array}
90    \label{EVAL PCG}
91    \end{equation}
92    The residual \index{residual}  $r=B A^{-1} G - S p$ is given as
93    \begin{equation}
94    r=B A^{-1} (G - B^* p) = Bv \mbox{ with } v = A^{-1}(G-B^{*}p)
95    \end{equation}
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    \begin{equation}
98    (p,(v,Bv))=\int\hackscore{\Omega} p \cdot Bv \; dx
99    \end{equation}
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  \begin{equation}  \begin{equation}
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  \end{equation}  \end{equation}
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  \begin{equation}  \begin{equation}
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  \end{equation}  \end{equation}
121  This scheme is called the Uzawa scheme.  We use the inner product induced by the norm
122    \begin{equation}
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    \end{equation}
126    In PDE form~\ref{COUPLES SADDLEPOINT iteration} takes the form
127  \begin{equation}  \begin{equation}
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  \end{equation}  \end{equation}
To accelerate the convergence we are using the restarted $GMRES$ method using the norm
\begin{equation}
\|(v,p)\|^2= \int\hackscore{\Omega} \eta^2 v\hackscore{i,j}^2 + p^2 \; dx
\end{equation}
or alternatively the $PCG$ method on the pressure only using the norm
\begin{equation}
\|p\|^2= \int\hackscore{\Omega} \frac{1}{\eta^2} p^2 \; dx
\end{equation}
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)*[0.,1.]+whereZero(x-1))*[1.,1]
206    v=Vector(0.,Solution(dom))
207    v+=whereZero(x-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    \begin{equation}\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    \end{equation}
221    with the boundary conditions
222    \begin{equation}\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    \end{equation}
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    \begin{equation}
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    \end{equation}
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) \\  \begin{equation}
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    \end{equation}
244    We set
245    \begin{equation}
246    V = \{ q \in H^{1}(\Omega) : q=0 \mbox{ on } \Gamma\hackscore{D} \}
247    \end{equation}
248    and
249    \begin{equation}
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    \end{equation}
252    and define the operator $Q: V \rightarrow (L^2(\Omega))^{d}$ defined by
253    \begin{equation}
254    (Qp)\hackscore{i} = \kappa\hackscore{ij} p\hackscore{,j}
255    \end{equation}
256    and the operator $D: W \rightarrow L^2(\Omega)$ defined by
257    \begin{equation}
258    Dv = v\hackscore{k,k}
259    \end{equation}
260    In operator notation the Darcy problem~\ref{DARCY PROBLEM} is written in the form
261    \begin{equation}
262    \begin{array}{rcl}
263    u + Qp & = & g \\
264    Du & = & f
265    \end{array}
266    \end{equation}
267    We solve this equation by minimising the functional
268    \begin{equation}
269    J(u,p):=\|u + Qp - g\|^2\hackscore{0} + \|Du-f\|\hackscore{0}^2
270    \end{equation}
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    \begin{equation}
274    ( v + Qq , u + Qp - g) + (Dv,Du-f) =0
275    \end{equation}
276    for all $v\in W$ and $q \in V$.which translates back into operator notation
277    \begin{equation}
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    \end{equation}
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    \begin{equation}
289    v= (I+D^*D)^{-1} (D^*f + g - Qp)
290    \end{equation}
291    (notice that $(I+D^*D)$ is coercive in $W$) which is inserted into the second equation
292    \begin{equation}
293    Q^* (I+D^*D)^{-1} (D^*f + g - Qp) + Q^* Q p = Q^* g
294    \end{equation}
295    which is
296    \begin{equation}
297    Q^* ( I - (I+D^*D)^{-1} ) Q p = Q^* ( g -(I+D^*D)^{-1} (D^*f + g) )
298    \end{equation}
299    We use the PCG method to solve this. The residual $r$ ($\in V^*$) is given as
300    \begin{equation}
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    \end{equation}
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    \begin{equation}\label{UPDATE W}
312    (I+D^*D)w = Qp
313    \end{equation}
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  \begin{equation}  \begin{equation}
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  \begin{equation}\label{IKM-EQU-6}  \begin{equation}\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  \end{equation}  \end{equation}
402  With  With
403  \begin{equation}\label{IKM-EQU-8}  \begin{equation}\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  \begin{equation}\label{IKM-EQU-10}  \begin{equation}\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  \end{equation}  \end{equation}
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  \begin{equation}\label{IKM-EQU-1ib}  \begin{equation}\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  \end{equation}  \end{equation}
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  \begin{equation}  \begin{equation}
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  \end{equation}  \end{equation}
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  \begin{equation}  \begin{equation}
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  \end{equation}  \end{equation}
474    where  $\bar{\eta}\hackscore{eff}$ is the caracteristic viscosity, for instance:
475    \begin{equation}
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    \end{equation}
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    \begin{equation}
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    \end{equation}
488    The Newton scheme takes the form
489    \begin{equation}
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    \end{equation}
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    \begin{equation}
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    \end{equation}
502    \begin{equation}\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    \end{equation}
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.1966 changed lines Added in v.2100