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

revision 2252 by gross, Fri Feb 6 07:51:28 2009 UTC revision 2432 by gross, Wed May 20 06:06:20 2009 UTC
# Line 183  return the tolerance for the involved PD Line 183  return the tolerance for the involved PD
183   The following script \file{lit\hackscore driven\hackscore cavity.py}   The following script \file{lit\hackscore driven\hackscore cavity.py}
184  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
185  illustrates the usage of the \class{StokesProblemCartesian} class to solve  illustrates the usage of the \class{StokesProblemCartesian} class to solve
186  the lit driven cavity problem~\cite{LITDRIVENCAVITY}:  the lit driven cavity problem:
187  \begin{python}  \begin{python}
188  from esys.escript import *  from esys.escript import *
189  from esys.finley import Rectangle  from esys.finley import Rectangle
# Line 281  Q^*u  + Q^*Q p & = & Q^*g \\ Line 281  Q^*u  + Q^*Q p & = & Q^*g \\
281  \end{array}  \end{array}
282  \end{equation}  \end{equation}
283  where $D^*$ and $Q^*$ denote the adjoint operators.  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  In~\cite{LEASTSQUARESFEM1994} 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$)    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  The approach we are taking is to eliminate the velocity $u$ from the problem. Assuming that $p$ is known we have
288  \begin{equation}  \begin{equation}\label{DARCY V FORM}
289  v= (I+D^*D)^{-1} (D^*f + g - Qp)  v= (I+D^*D)^{-1} (D^*f + g - Qp)
290  \end{equation}  \end{equation}
291  (notice that $(I+D^*D)$ is coercive in $W$) which is inserted into the second equation  (notice that $(I+D^*D)$ is coercive in $W$) which is inserted into the second equation
# Line 300  We use the PCG \index{linear solver!PCG} Line 300  We use the PCG \index{linear solver!PCG}
300  \begin{equation}  \begin{equation}
301  \begin{array}{rcl}  \begin{array}{rcl}
302  r & = & Q^*  \left( g -(I+D^*D)^{-1} (D^*f + g) - Qp + (I+D^*D)^{-1}Q p \right)\\  r & = & Q^*  \left( g -(I+D^*D)^{-1} (D^*f + g) - Qp + (I+D^*D)^{-1}Q p \right)\\
303  & =&  Q^* \left( - Qp - (I+D^*D)^{-1} (D^*f + g - Qp) \right) \\  & =&  Q^* \left( g- Qp - (I+D^*D)^{-1} (D^*f + g - Qp) \right) \\
304  & =&  Q^* \left( g - Qp - v \right)  & =&  Q^* \left( g - Qp - v \right)
305  \end{array}  \end{array}
306  \end{equation}  \end{equation}
307  So in a partical implementation we use the pair $(Qp,v)$ to represent the residual. This will save the  So in a partical implementation we use $\hat{r}=g-Qp-v$ to represent the residual.
308  reconstruction of the velocity $v$. In this notation the right hand side is given as  The evaluation of the iteration operator for a given $p$ is then
309  $(0,(I+D^*D)^{-1} (D^*f + g))$. The evaluation of the iteration operator for a given $p$ is then  returning $Qp+v$ where $v$ is the solution of
returning $(Qp,w)$ where $w$ is the solution of
310  \begin{equation}\label{UPDATE W}  \begin{equation}\label{UPDATE W}
311  (I+D^*D)w = Qp  (I+D^*D)v = Qp
312  \end{equation}  \end{equation}
313  We use $Q^*Q$ as a a preconditioner for the iteration operator $Q^* ( I - (I+D^*D)^{-1} ) Q$.  We use $(Q^*Q)^{-1}$ as a preconditioner for the iteration operator $Q^* ( I - (I+D^*D)^{-1} ) Q$. So the application of the preconditioner to $\hat{r}$ representing the residual is given by solving
314    implemented by solving
315  The iteration PCG \index{linear solver!PCG}\index{PCG} is terminated if  \begin{equation}\label{UPDATE P}
316    Q^*Q q  = Q^*\hat{r}
317    \end{equation}
318    The residual norm used in the PCG is given as
319    \begin{equation}\label{DARCY R NORM}
320    \|r\|\hackscore{PCG}^2 = \int r \cdot (Q^*Q)^{-1} r \; dx =\int \hat{r} \cdot Q (Q^*Q)^{-1} Q^* \hat{r} \; dx \approx
321    \|\hat{r}\|\hackscore{0}^2
322    \end{equation}
323    The iteration is terminated if
324  \begin{equation}\label{DARCY STOP}  \begin{equation}\label{DARCY STOP}
325  \int r \cdot (Q^*Q)^{-1} r \; dx \le \mbox{ATOL}^2  \|r\|\hackscore{PCG} \le \mbox{ATOL}
326  \end{equation}  \end{equation}
327  where ATOL is a given absolute tolerance.  where we set
328  The initial residual $r\hackscore{0}$ is  \begin{equation}\label{DARCY ATOL DEF}
329  \begin{equation}\label{DARCY STOP 2}  \mbox{ATOL} = \mbox{atol} + \mbox{rtol} \cdot \left(\frac{1}{\|v\|\hackscore{0}} + \frac{1}{\|Qp\|\hackscore{0}} \right)^{-1}
330  r\hackscore{0}=Q^* \left( g - v\hackscore{ref} \right) \mbox{ with } v\hackscore{ref} = (I+D^*D)^{-1} (D^*f + g)  \end{equation}
331  \end{equation}  where rtol is a given relative tolerance and $\mbox{atol}$ is a given absolute tolerance (typically $=0$).
332  so the  Notice that if $Qp$ and $v$ both are zero, the pair $(0,p)$ is a solution.
333  \begin{equation}\label{DARCY NORM 0}  The problem is that ATOL is depending on the solution $p$ (and $v$ calculated form~\ref{DARCY V FORM}). In partcice one use the initial guess for $p$
334  \int r\hackscore0 \cdot (Q^*Q)^{-1} r\hackscore0 \; dx = \int \left( g - v\hackscore{ref} \right)  \cdot  Q  p\hackscore{ref} \; dx \mbox{ with }p\hackscore{ref} = (Q^*Q)^{-1} Q^* \left( g - v\hackscore{ref} \right)  to get a first value for ATOL. If the stopping crierion is met in the PCG iteration, a new $v$ is calculated from the current pressure approximation and ATOL is recalculated. If \ref{DARCY STOP} is still fullfilled the calculation is terminated and $(v,p)$ is returned. Otherwise PCG is restarted with a new ATOL.
\end{equation}
So we set
\begin{equation}\label{DARCY NORM 1}
ATOL = atol + rtol \cdot \max(|g - v\hackscore{ref}|\hackscore{0}, |Q p\hackscore{ref} |\hackscore{0} )
\end{equation}
where atol and rtol a given absolute and relative tolerances, respectively. The reference flux $v\hackscore{ref}$
and reference pressure $p\hackscore{ref}$ may be calcualated from their definition which would require to solve to
PDEs but in a practical application estimates can be used for instance solutions from previous time steps or for simplified scenarious (e.g. constant permability).
335
336  \subsection{Functions}  \subsection{Functions}
337  \begin{classdesc}{DarcyFlow}{domain}  \begin{classdesc}{DarcyFlow}{domain}
# Line 353  The method will try to cast the given va Line 352  The method will try to cast the given va
352  \Data class objects.  \Data class objects.
353  \end{methoddesc}  \end{methoddesc}
354
355  \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{atol=0,\optional{rtol=1e-8,\optional{p_ref=None,\optional{v_ref=None}}}}}  \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e-4}}
356  sets the absolute tolerance ATOL according to~\ref{DARCY NORM 1}. If \var{p_ref} is not present $0$ is used.  sets the relative tolerance \mbox{rtol} in \ref{DARCY ATOL DEF}.
If \var{v_ref} is not present $0$ is used. If the final result ATOL is not positive an exception is thrown.
357  \end{methoddesc}  \end{methoddesc}
358
359    \begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}}
360    sets the absolute tolerance \mbox{atol} in \ref{DARCY ATOL DEF}.
361    \end{methoddesc}
362
363    \begin{methoddesc}[DarcyFlow]{setSubProblemTolerance}{\optional{rtol=None}}
364    sets the relative tolerance used to solve the involved PDEs. If no argument is given,
365    the square of the current relative tolerance is used. The sub-problem tolerance should be choosen as large as possible to minimize the compute time. However, a too large value for the sub-problem tolerance may lead to slow convergence or even dibergence in the outer iteration.
366    \end{methoddesc}
367
368  \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False \optional{sub_rtol=1.e-8}}}}  \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False \optional{sub_rtol=1.e-8}}}}
369  solves the problem. and returns approximations for the flux $v$ and the pressure $p$.  solves the problem. and returns approximations for the flux $v$ and the pressure $p$.
370  \var{u0} and \var{p0} define initial guess for flux and pressure. Values marked  \var{u0} and \var{p0} define initial guess for flux and pressure. Values marked
371  by positive values \var{location_of_fixed_flux} and \var{location_of_fixed_pressure}, respectively, are kept unchanged.  by positive values \var{location_of_fixed_flux} and \var{location_of_fixed_pressure}, respectively, are kept unchanged.
\var{sub_rtol} defines the tolerance used to solve the involved PDEs. \var{sub_rtol} needs to be choosen sufficiently small to ensure convergence but users need to keep in mind that a very small value for \var{sub_rtol} will result in a long compute time. Typically  $\var{sub_rtol}=\var{rtol}^2$ is a good choice if $\var{rtol}$ is not choosen too small.
372  \end{methoddesc}  \end{methoddesc}
373
374
375  \subsection{Example: Gravity Flow}  \subsection{Example: Gravity Flow}
376  later  later
377
378  %================================================  %\input{levelsetmodel}

%\begin{equation}
% \rho c\hackscore{p} \left (\frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T \right ) = k \nabla^{2}T
% \label{HEAT EQUATION}
% \end{equation}
379
% where $\vec{v}$ is the velocity vector, $T$ is the temperature, $\rho$ is the density, $\eta$ is the viscosity, % % $c\hackscore{p}$ is the specific heat at constant pressure and $k$ is the thermal conductivity.
380
% \subsection{Description}

% \subsection{Method}
%
% \begin{classdesc}{TemperatureCartesian}{dom,theta=THETA,useSUPG=SUPG}
% \end{classdesc}

% \subsection{Benchmark Problem}
%===============================================================================================================

%=========================================================
\input{levelsetmodel}

% \section{Drucker Prager Model}
381
382  \section{Isotropic Kelvin Material \label{IKM}}  \section{Isotropic Kelvin Material \label{IKM}}
383  As proposed by Kelvin~\ref{KELVN} material strain $D\hackscore{ij}=v\hackscore{i,j}+v\hackscore{j,i}$ can be decomposed into  As proposed by Kelvin~\cite{Muhlhaus2005} material strain $D\hackscore{ij}=\frac{1}{2}(v\hackscore{i,j}+v\hackscore{j,i})$ can be decomposed into
384  an elastic part $D\hackscore{ij}^{el}$ and visco-plastic part $D\hackscore{ij}^{vp}$:  an elastic part $D\hackscore{ij}^{el}$ and visco-plastic part $D\hackscore{ij}^{vp}$:
385  \begin{equation}\label{IKM-EQU-2}  \begin{equation}\label{IKM-EQU-2}
386  D\hackscore{ij}=D\hackscore{ij}^{el}+D\hackscore{ij}^{vp}  D\hackscore{ij}=D\hackscore{ij}^{el}+D\hackscore{ij}^{vp}
# Line 421  betwee the the strain in material $q$ Line 404  betwee the the strain in material $q$
404  \eta^{q}=\eta^{q}\hackscore{N} \left(\frac{\tau}{\tau\hackscore{t}^q}\right)^{\frac{1}{n^{q}}-1}  \eta^{q}=\eta^{q}\hackscore{N} \left(\frac{\tau}{\tau\hackscore{t}^q}\right)^{\frac{1}{n^{q}}-1}
405  \mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'\hackscore{ij} \sigma'\hackscore{ij}}  \mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'\hackscore{ij} \sigma'\hackscore{ij}}
406  \end{equation}  \end{equation}
407  for a given power law coefficients $n^{q}\ge1$ and transition stresses $\tau\hackscore{t}^q$, see~\ref{POERLAW}.  for a given power law coefficients $n^{q}\ge1$ and transition stresses $\tau\hackscore{t}^q$, see~\cite{Muhlhaus2005}.
408  Notice that $n^{q}=1$ gives a constant viscosity.  Notice that $n^{q}=1$ gives a constant viscosity.
409  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:
410  \begin{equation}\label{IKM-EQU-6}  \begin{equation}\label{IKM-EQU-6}
411  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}} \;.
412  \end{equation}  \end{equation}
413  With  and finally with~\ref{IKM-EQU-2}
414  \begin{equation}\label{IKM-EQU-8}  \begin{equation}\label{IKM-EQU-2bb}
415  \dot{\gamma}=\sqrt{2 D\hackscore{ij} D\hackscore{ij}}  D\hackscore{ij}'=\frac{1}{2 \eta^{vp}} \sigma'\hackscore{ij}+\frac{1}{2 \mu} \dot{\sigma}\hackscore{ij}'
\end{equation}
one gets
\begin{equation}\label{IKM-EQU-8b}
\tau = \eta^{vp} \dot{\gamma}^{vp} \;.
416  \end{equation}  \end{equation}
417  With the Drucker-Prager cohesion factor $\tau\hackscore{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$ we want to achieve  The total stress $\tau$ needs to fullfill the yield condition \index{yield condition}
418  \begin{equation}\label{IKM-EQU-8c}  \begin{equation}\label{IKM-EQU-8c}
419  \tau \le \tau\hackscore{Y} + \beta \; p  \tau \le \tau\hackscore{Y} + \beta \; p
420  \end{equation}  \end{equation}
421  which leads to the condition  with the Drucker-Prager \index{Druck-Prager} cohesion factor \index{cohesion factor} $\tau\hackscore{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$.
\begin{equation}\label{IKM-EQU-8d}
\eta^{vp} \le \frac{\tau\hackscore{Y} + \beta \; p}{ \dot{\gamma}^{vp}} \; .
\end{equation}
Therefore we modify the definition of $\eta^{vp}$ to the form
\begin{equation}\label{IKM-EQU-6b}
\frac{1}{\eta^{vp}}=\max(\sum\hackscore{q} \frac{1}{\eta^{q}}, \frac{\dot{\gamma}^{vp}} {\tau\hackscore{Y} + \beta \; p})
\end{equation}
422  The deviatoric stress needs to fullfill the equilibrion equation  The deviatoric stress needs to fullfill the equilibrion equation
423  \begin{equation}\label{IKM-EQU-1}  \begin{equation}\label{IKM-EQU-1}
424  -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}  -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}
425  \end{equation}  \end{equation}
426  where $F\hackscore{j}$ is a given external fource. We assume an incompressible media:  where $F\hackscore{j}$ is a given external fource. We assume an incompressible media:
427  \begin{equation}\label{IKM-EQU-2}  \begin{equation}\label{IKM-EQU-2bbb}
428  -v\hackscore{i,i}=0  -v\hackscore{i,i}=0
429  \end{equation}  \end{equation}
430  Natural boundary conditions are taken in the form  Natural boundary conditions are taken in the form
# Line 468  where the index $i$ may depend on the lo Line 440  where the index $i$ may depend on the lo
440  \subsection{Solution Method \label{IKM-SOLVE}}  \subsection{Solution Method \label{IKM-SOLVE}}
441  By using a first order finite difference approximation wit step size $dt>0$~\ref{IKM-EQU-3} get the form  By using a first order finite difference approximation wit step size $dt>0$~\ref{IKM-EQU-3} get the form
442  \begin{equation}\label{IKM-EQU-3b}  \begin{equation}\label{IKM-EQU-3b}
443  D\hackscore{ij}'^{el}=\frac{1}{2 \mu dt } \left( \sigma\hackscore{ij}' - \sigma\hackscore{ij}^{'-} \right)  \dot{\sigma}\hackscore{ij}=\frac{1}{dt } \left( \sigma\hackscore{ij} - \sigma\hackscore{ij}^{-} \right)
444    \end{equation}
445    and
446    \begin{equation}\label{IKM-EQU-2b}
447    D\hackscore{ij}'=\left(\frac{1}{2 \eta^{vp}} + \frac{1}{2 \mu dt}\right) \sigma\hackscore{ij}'-\frac{1}{2 \mu dt } \sigma\hackscore{ij}^{-'}
448  \end{equation}  \end{equation}
449  where $\sigma\hackscore{ij}^{'-}$ is the deviatoric stress at the precious time step.  where $\sigma\hackscore{ij}^{-}$ is the stress at the precious time step. With
450  Now we can combine equations~\ref{IKM-EQU-2}, \ref{IKM-EQU-3b} and~\ref{IKM-EQU-6b} to get  \begin{equation}\label{IKM-EQU-2c}
451    \dot{\gamma} = \sqrt{ 2 \left( D\hackscore{ij}' +
452    \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{-'}\right)^2}
453    \end{equation}
454    we have
455    \begin{equation}
456    \tau = \eta\hackscore{eff} \cdot \dot{\gamma}
457    \end{equation}
458    where
459    \begin{equation}
460    \eta\hackscore{eff}= min( \left(\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}\right)^{-1}
461    , \eta\hackscore{max}) \mbox{ with }
462    \eta\hackscore{max} = \left\{
463    \begin{array}{rcl}
464    \frac{\tau\hackscore{Y} + \beta \; p}{\dot{\gamma}} & & \dot{\gamma}>0 \\
465    &\mbox{ if } \\
466    \infty & & \mbox{otherwise}
467    \end{array}
468    \right.
469    \end{equation}
470    The upper bound $\eta\hackscore{max}$ makes sure that yield condtion~\ref{IKM-EQU-8c} holds. With this setting the eqaution \ref{IKM-EQU-2b} takes the form
471  \begin{equation}\label{IKM-EQU-10}  \begin{equation}\label{IKM-EQU-10}
472  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff}  \left( D\hackscore{ij}' +  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff}  \left( D\hackscore{ij}' +
473  \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)  \mbox{ with }  \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)
\frac{1}{\eta\hackscore{eff}}=\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
474  \end{equation}  \end{equation}
Notice that $\eta\hackscore{eff}$ is a function of diatoric stress $\sigma\hackscore{ij}'$.
475  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
476  \begin{equation}\label{IKM-EQU-1ib}  \begin{equation}\label{IKM-EQU-1ib}
477  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})
478  \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+  \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+
479  \frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij,j}^{'-}   \left(\frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij}^{'-} \right)\hackscore{,j}
480  \end{equation}  \end{equation}
481  Combining this with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a  Combining this with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a
482  Stokes problem as discussed in section~\ref{STOKES SOLVE} in each time step.  Stokes problem as discussed in section~\ref{STOKES SOLVE} in each time step.
483  In oder to perform step~\ref{IKM iteration 2} we need to calculate the $\eta\hackscore{eff}$ which
484  is a function of $\sigma\hackscore{ij}$ via $\tau$.  To get $\tau$ and $\eta\hackscore{eff}$ we need to solve the  If we set
485  non-linear equation  \begin{equation}\label{IKM-EQU-44}
486  \begin{equation}  \frac{1}{\eta(\tau)}= \frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
487  \tau = \eta\hackscore{eff} \cdot \dot{\gamma}\hackscore{total} \mbox{ with }  \end{equation}
488  \dot{\gamma}\hackscore{total} = \sqrt{ 2 \left( D\hackscore{ij}' +  we need to solve the nonlinear problem
489  \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)^2}  \begin{equation}
490  \label{IKM iteration 5}  \eta\hackscore{eff} -  min(\eta( \dot{\gamma} \cdot \eta\hackscore{eff})
491  \end{equation}  , \eta\hackscore{max}) =0
492  The Newton scheme takes the form  \end{equation}
493  \begin{equation}  We use the Newton-Raphson Scheme \index{Newton-Raphson scheme} to solve this problem
494  \tau\hackscore{n+1} = \min(\tau\hackscore{n} - \frac{\tau\hackscore{n} - \eta\hackscore{eff}  \cdot \dot{\gamma}\hackscore{total}}{1 - \eta\hackscore{eff}'  \cdot  \dot{\gamma}\hackscore{total}}, \tau\hackscore{Y} + \beta \; p)  \begin{equation}\label{IKM-EQU-45}
495  \label{IKM iteration 6}  \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
496  \end{equation}  \eta\hackscore{eff}^{(n)} -
497  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$ (?)). We have  \frac{\eta\hackscore{eff}^{(n)} - \eta( \tau^{(n)}) }
498    {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
499    =\min(\eta\hackscore{max},
500    \frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot  \eta'( \tau^{(n)} )  }
501    {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
502    \end{equation}
503    where $\eta'$ denotes the derivative of $\eta$ with respect of $\tau$
504    and $\tau^{(n)} = \dot{\gamma} \cdot \eta\hackscore{eff}^{(n)}$.
505
506    Looking at the evaluation of $\eta$ in~\ref{IKM-EQU-44} it makes sense formulate
507    the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$.
508    In fact we have
509  \begin{equation}  \begin{equation}
510  \eta\hackscore{eff}' = - \eta\hackscore{eff}^2 \left(\frac{1}{\eta\hackscore{eff}}\right)'  \eta' = - \frac{\Theta'}{\Theta^2}
511  \mbox{ with }  \mbox{ with }
512  \left(\frac{1}{\eta\hackscore{eff}}\right)' = \sum\hackscore{q} \left(\frac{1}{\eta^{q}} \right)'  \Theta' = \sum\hackscore{q} \left(\frac{1}{\eta^{q}} \right)'
513  \label{IKM iteration 7}  \label{IKM iteration 7}
514  \end{equation}  \end{equation}
515  \begin{equation}\label{IKM-EQU-5XX}  As
516    \begin{equation}\label{IKM-EQU-47}
517  \left(\frac{1}{\eta^{q}} \right)'  \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}}}}  = \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}\frac{1}{\eta^{q}}  = \frac{1-\frac{1}{n^{q}}}{\eta^{q}}\frac{1}{\tau}
520  \end{equation}  \end{equation}
521  Notice that allways $\eta\hackscore{eff}'\le 0$ which makes the denomionator in~\ref{IKM iteration 6}  we have
522  positive.  \begin{equation}\label{IKM-EQU-48}
523    \Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum\hackscore{q}\frac{1-\frac{1}{n^{q}}}{\eta^{q}}
524    \end{equation}
526    \begin{equation}\label{IKM-EQU-49}
527    \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
528    \eta\hackscore{eff}^{(n)}
529    \frac{\Theta^{(n)}  + \omega^{(n)}  }
530    {\eta\hackscore{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })
531    \end{equation}
532
533
534    \subsection{Functions}
535
536    #=============================================================================
537    \begin{classdesc}{IncompressibleIsotropicFlowCartesian}{
538    domain
539    \optional{, stress=0
540    \optional{, v=0
541    \optional{, p=0
542    \optional{, t=0
543    \optional{, numMaterials=1
544    \optional{, verbose=True}}}}}}}
545    opens an incompressible, isotropic flow problem in Cartesian cooridninates.
546    \var{stress},
547    \var{v},
548    \var{p}, and
549    \var{t} set the initial deviatoric stress, velocity, pressure and time.
550    \var{numMaterials} specifies the number of materials used in the power law
551    model. Some progress information are printed if \var{verbose} is set to
552    \True.
553    \end{classdesc}
554
555     setExternals(self, F=None, f=None, fixed_v_mask=None, v_boundary=None):
557    assigns values to external forces and boundary conditions. Between two calls only variables with a new values need to be set.
558
559    In any call all values must be set.
560    \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,
561    \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.
562    The locations and compontents where the velocity is fixed are set by
563    the values of \var{fixed_u_mask}. The method will try to cast the given values to appropriate
564    \Data class objects.
565    \end{methoddesc}
566
567    \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p,
568    \optional{max_iter=20, \optional{verbose=False, \optional{usePCG=True}}}}
569    solves the problem and return approximations for velocity and pressure.
570    The arguments \var{v} and \var{p} define initial guess. The values of \var{v} marked
572    If \var{usePCG} is set to \True
573    reconditioned conjugate gradient method (PCG) \index{preconditioned conjugate gradient method!PCG}  scheme is used. Otherwise the problem is solved generalized minimal residual method (GMRES) \index{generalized minimal residual method!GMRES}. In most cases
574    the PCG scheme is more efficient.
575    \var{max_iter} defines the maximum number of iteration steps.
576    If \var{verbose} is set to \True informations on the progress of of the solver are printed.
577    \end{methoddesc}
578
579
580    \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}}
581    sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1.
582    \end{methoddesc}
583    \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{}
584    returns the current relative tolerance.
585    \end{methoddesc}
586    \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}}
587    sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the
588    absolute talerance is set to 0.
589    \end{methoddesc}
590    \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{}
591    sreturns the current absolute tolerance.
592    \end{methoddesc}
593    \begin{methoddesc}[StokesProblemCartesian]{setSubProblemTolerance}{\optional{rtol=None}}
594    sets the tolerance to solve the involved PDEs. The subtolerance \var{rtol} should not be choosen to large
595    in order to avoid feed back of errors in the subproblem solution into the outer iteration.
596    On the otherhand is choosen to small compute time is wasted.
597    If \var{rtol} is set to \var{None} the sub-tolerance is set automatically depending on the
598    tolerance choosen for the oter iteration.
599    \end{methoddesc}
600    \begin{methoddesc}[StokesProblemCartesian]{getSubProblemTolerance}{}
601    return the tolerance for the involved PDEs.
602    \end{methoddesc}
603
604
605    % \section{Drucker Prager Model}

Legend:
 Removed from v.2252 changed lines Added in v.2432