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

revision 2252 by gross, Fri Feb 6 07:51:28 2009 UTC revision 2502 by gross, Tue Jun 30 05:49:22 2009 UTC
# Line 28  where $\eta$ is the viscosity, $F\hacksc Line 28 where$\eta$is the viscosity,$F\hacksc
28  \end{equation}  \end{equation}
29  Natural boundary conditions are taken in the form  Natural boundary conditions are taken in the form
30  \begin{equation}\label{Stokes Boundary}  \begin{equation}\label{Stokes Boundary}
31  \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}  \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{j}
32  \end{equation}  \end{equation}
33  which can be overwritten by constraints of the form  which can be overwritten by constraints of the form
34  \begin{equation}\label{Stokes Boundary0}  \begin{equation}\label{Stokes Boundary0}
# Line 81  S p =  B A^{-1} G Line 81  S p =  B A^{-1} G
81  \end{equation}  \end{equation}
82  with the Schur complement \index{Schur complement} $S=BA^{-1}B^{*}$. This problem can be solved iteratively  with the Schur complement \index{Schur complement} $S=BA^{-1}B^{*}$. This problem can be solved iteratively
83  with the preconditioner $\hat{S}$ defined as $q=\hat{S}^{-1}p$ by solving  with the preconditioner $\hat{S}$ defined as $q=\hat{S}^{-1}p$ by solving
84  \begin{equation}  \begin{equation} \label{P PREC}
85  \frac{1}{\eta}q = p  \frac{1}{\eta}q = p
86  \end{equation}  \end{equation}
87  see~\cite{ELMAN} for more details. Note that the residual for the current approximation $p$ is given as  see~\cite{ELMAN} for more details. Note that the residual for the current approximation $p$ is given as
# Line 128  ATOL = \epsilon \frac{\|\sqrt{v^{0}\hack Line 128  ATOL = \epsilon \frac{\|\sqrt{v^{0}\hack
128
129  \subsection{Functions}  \subsection{Functions}
130
132  opens the Stokes problem\index{Stokes problem} on the \Domain domain. The approximation  opens the Stokes problem\index{Stokes problem} on the \Domain domain. The approximation
133  order needs to be two.  order needs to be two.
135    the tolerances for all subproblems are set automatically.
136  \end{classdesc}  \end{classdesc}
137
# Line 142  the values of \var{fixed_u_mask}. The me Line 144  the values of \var{fixed_u_mask}. The me
144  \Data class objects.  \Data class objects.
145  \end{methoddesc}  \end{methoddesc}
146
147  \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p,  \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p
148  \optional{max_iter=20, \optional{verbose=False, \optional{usePCG=True}}}}  \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}}
149  solves the problem and return approximations for velocity and pressure.  solves the problem and return approximations for velocity and pressure.
150  The arguments \var{v} and \var{p} define initial guess. The values of \var{v} marked  The arguments \var{v} and \var{p} define initial guess. The values of \var{v} marked
# Line 151  If \var{usePCG} is set to \True Line 153  If \var{usePCG} is set to \True
153  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  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
154  the PCG scheme is more efficient.  the PCG scheme is more efficient.
155  \var{max_iter} defines the maximum number of iteration steps.  \var{max_iter} defines the maximum number of iteration steps.
156
157  If \var{verbose} is set to \True informations on the progress of of the solver are printed.  If \var{verbose} is set to \True informations on the progress of of the solver are printed.
158  \end{methoddesc}  \end{methoddesc}
159
# Line 165  returns the current relative tolerance. Line 168  returns the current relative tolerance.
168  sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the  sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the
169  absolute talerance is set to 0.  absolute talerance is set to 0.
170  \end{methoddesc}  \end{methoddesc}
171
172  \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{}  \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{}
173  sreturns the current absolute tolerance.  sreturns the current absolute tolerance.
174  \end{methoddesc}  \end{methoddesc}
175  \begin{methoddesc}[StokesProblemCartesian]{setSubProblemTolerance}{\optional{rtol=None}}
176  sets the tolerance to solve the involved PDEs. The subtolerance \var{rtol} should not be choosen to large  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{}
177  in order to avoid feed back of errors in the subproblem solution into the outer iteration.  returns the solver options used  solve the equations~(\ref{V CALC}) for velocity.
178  On the otherhand is choosen to small compute time is wasted.  \end{methoddesc}
179  If \var{rtol} is set to \var{None} the sub-tolerance is set automatically depending on the
180  tolerance choosen for the oter iteration.  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{}
181    returns the solver options used  solve the equation~(\label{P PREC}) for pressure.
182  \end{methoddesc}  \end{methoddesc}
183  \begin{methoddesc}[StokesProblemCartesian]{getSubProblemTolerance}{}
184  return the tolerance for the involved PDEs.  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{}
185    set the solver options for solving the equation to project the divergence of the velocity onto the function space of pressure.
186  \end{methoddesc}  \end{methoddesc}
187
188
189  \subsection{Example: Lit Driven Cavity}  \subsection{Example: Lit Driven Cavity}
190   The following script \file{lit\hackscore driven\hackscore cavity.py}   The following script \file{lit\hackscore driven\hackscore cavity.py}
191  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
192  illustrates the usage of the \class{StokesProblemCartesian} class to solve  illustrates the usage of the \class{StokesProblemCartesian} class to solve
193  the lit driven cavity problem~\cite{LITDRIVENCAVITY}:  the lit driven cavity problem:
194  \begin{python}  \begin{python}
195  from esys.escript import *  from esys.escript import *
196  from esys.finley import Rectangle  from esys.finley import Rectangle
# Line 281  Q^*u  + Q^*Q p & = & Q^*g \\ Line 288  Q^*u  + Q^*Q p & = & Q^*g \\
288  \end{array}  \end{array}
289  \end{equation}  \end{equation}
290  where $D^*$ and $Q^*$ denote the adjoint operators.  where $D^*$ and $Q^*$ denote the adjoint operators.
291  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
292  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$)
293
294  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
295  \begin{equation}  \begin{equation}\label{DARCY V FORM}
296  v= (I+D^*D)^{-1} (D^*f + g - Qp)  v= (I+D^*D)^{-1} (D^*f + g - Qp)
297  \end{equation}  \end{equation}
298  (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 307  We use the PCG \index{linear solver!PCG}
307  \begin{equation}  \begin{equation}
308  \begin{array}{rcl}  \begin{array}{rcl}
309  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)\\
310  & =&  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) \\
311  & =&  Q^* \left( g - Qp - v \right)  & =&  Q^* \left( g - Qp - v \right)
312  \end{array}  \end{array}
313  \end{equation}  \end{equation}
314  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.
315  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
316  $(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
317  \begin{equation}\label{UPDATE W}  \begin{equation}\label{UPDATE W}
318  (I+D^*D)w = Qp  (I+D^*D)v = Qp
319  \end{equation}  \end{equation}
320  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
321    implemented by solving
322  The iteration PCG \index{linear solver!PCG}\index{PCG} is terminated if  \begin{equation}\label{UPDATE P}
323    Q^*Q q  = Q^*\hat{r}
324    \end{equation}
325    The residual norm used in the PCG is given as
326    \begin{equation}\label{DARCY R NORM}
327    \|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
328    \|\hat{r}\|\hackscore{0}^2
329    \end{equation}
330    The iteration is terminated if
331  \begin{equation}\label{DARCY STOP}  \begin{equation}\label{DARCY STOP}
332  \int r \cdot (Q^*Q)^{-1} r \; dx \le \mbox{ATOL}^2  \|r\|\hackscore{PCG} \le \mbox{ATOL}
333  \end{equation}  \end{equation}
334  where ATOL is a given absolute tolerance.  where we set
335  The initial residual $r\hackscore{0}$ is  \begin{equation}\label{DARCY ATOL DEF}
336  \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}
337  r\hackscore{0}=Q^* \left( g - v\hackscore{ref} \right) \mbox{ with } v\hackscore{ref} = (I+D^*D)^{-1} (D^*f + g)  \end{equation}
338  \end{equation}  where rtol is a given relative tolerance and $\mbox{atol}$ is a given absolute tolerance (typically $=0$).
339  so the  Notice that if $Qp$ and $v$ both are zero, the pair $(0,p)$ is a solution.
340  \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$
341  \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).
342
343  \subsection{Functions}  \subsection{Functions}
345  opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.  opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.
346    If \var{adaptSubTolerance} is set to \True,
347    the relative tolerances for solving~(\ref{DARCY V FORM}),~(\ref{UPDATE W})
348    and~(\ref{UPDATE P}) are set automatically.
349  \end{classdesc}  \end{classdesc}
350
351  \begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None, \optional{permeability=None}}}}}}  \begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None, \optional{permeability=None}}}}}}
352  assigns values to the model parameters. Values can be assigned using various calls - in particular  assigns values to the model parameters. Values can be assigned using various calls - in particular
353  in a time dependend problem only values that change over time needs to be reset. The permability can be defined as scalar (isotropic), a vector (orthotropic) or a matrix (anisotropic).  in a time dependend problem only values that change over time needs to be reset. The permability can be defined as scalar (isotropic), a vector (orthotropic) or a matrix (anisotropic).
# Line 353  The method will try to cast the given va Line 363  The method will try to cast the given va
363  \Data class objects.  \Data class objects.
364  \end{methoddesc}  \end{methoddesc}
365
366  \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}}
367  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.
368  \end{methoddesc}  \end{methoddesc}
369
370    \begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}}
371    sets the absolute tolerance \mbox{atol} in \ref{DARCY ATOL DEF}.
372    \end{methoddesc}
373
374    \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{}
375    Returns the solver options used to solve the flux problems~(\ref{DARCY V FORM}) and~(\ref{UPDATE W}). Use the returned \SolverOptions object to control the solution algorithms. If the adaption of subtolerance is choosen, the tolerance will
376    be overwritten before the solver is called.
377    \end{methoddesc}
378
379    \begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{}
380    Returns the solver options used to solve the pressure problems~(\ref{UPDATE P}).
381    Use the returned \SolverOptions object to control the solution algorithms. If the adaption of subtolerance is choosen, the tolerance will
382    be overwritten before the solver is called.
383    \end{methoddesc}
384
385  \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}}}
386  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$.
387  \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
388  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{max_iter} sets the maximum number of iterations steps allowed for solving the coupled problem.
\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.
389  \end{methoddesc}  \end{methoddesc}
390
391
392  \subsection{Example: Gravity Flow}  \subsection{Example: Gravity Flow}
393  later  later
394
395  %================================================  %\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}

% 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.

% \subsection{Description}

% \subsection{Method}
%
% \begin{classdesc}{TemperatureCartesian}{dom,theta=THETA,useSUPG=SUPG}
% \end{classdesc}
396
% \subsection{Benchmark Problem}
%===============================================================================================================
397
%=========================================================
\input{levelsetmodel}

% \section{Drucker Prager Model}
398
399  \section{Isotropic Kelvin Material \label{IKM}}  \section{Isotropic Kelvin Material \label{IKM}}
400  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
401  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}$:
402  \begin{equation}\label{IKM-EQU-2}  \begin{equation}\label{IKM-EQU-2}
403  D\hackscore{ij}=D\hackscore{ij}^{el}+D\hackscore{ij}^{vp}  D\hackscore{ij}=D\hackscore{ij}^{el}+D\hackscore{ij}^{vp}
# Line 418  D\hackscore{ij}^{q'}=\frac{1}{2 \eta^{q} Line 418  D\hackscore{ij}^{q'}=\frac{1}{2 \eta^{q}
418  where $\eta^{q}$ is the viscosity of material $q$. We assume the following  where $\eta^{q}$ is the viscosity of material $q$. We assume the following
419  betwee the the strain in material $q$  betwee the the strain in material $q$
420  \begin{equation}\label{IKM-EQU-5b}  \begin{equation}\label{IKM-EQU-5b}
421  \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)^{1-n^{q}}
422  \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}}
423  \end{equation}  \end{equation}
424  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}.
425  Notice that $n^{q}=1$ gives a constant viscosity.  Notice that $n^{q}=1$ gives a constant viscosity.
426  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:
427  \begin{equation}\label{IKM-EQU-6}  \begin{equation}\label{IKM-EQU-6}
428  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}} \;.
429  \end{equation}  \end{equation}
430  With  and finally with~\ref{IKM-EQU-2}
431  \begin{equation}\label{IKM-EQU-8}  \begin{equation}\label{IKM-EQU-2bb}
432  \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} \;.
433  \end{equation}  \end{equation}
434  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}
435  \begin{equation}\label{IKM-EQU-8c}  \begin{equation}\label{IKM-EQU-8c}
436  \tau \le \tau\hackscore{Y} + \beta \; p  \tau \le \tau\hackscore{Y} + \beta \; p
437  \end{equation}  \end{equation}
438  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}
439  The deviatoric stress needs to fullfill the equilibrion equation  The deviatoric stress needs to fullfill the equilibrion equation
440  \begin{equation}\label{IKM-EQU-1}  \begin{equation}\label{IKM-EQU-1}
441  -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}  -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}
442  \end{equation}  \end{equation}
443  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:
444  \begin{equation}\label{IKM-EQU-2}  \begin{equation}\label{IKM-EQU-2bbb}
445  -v\hackscore{i,i}=0  -v\hackscore{i,i}=0
446  \end{equation}  \end{equation}
447  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 457  where the index $i$ may depend on the lo
457  \subsection{Solution Method \label{IKM-SOLVE}}  \subsection{Solution Method \label{IKM-SOLVE}}
458  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
459  \begin{equation}\label{IKM-EQU-3b}  \begin{equation}\label{IKM-EQU-3b}
460  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)
461    \end{equation}
462    and
463    \begin{equation}\label{IKM-EQU-2b}
464    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}^{-'}
465    \end{equation}
466    where $\sigma\hackscore{ij}^{-}$ is the stress at the precious time step. With
467    \begin{equation}\label{IKM-EQU-2c}
468    \dot{\gamma} = \sqrt{ 2 \left( D\hackscore{ij}' +
469    \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{-'}\right)^2}
470    \end{equation}
471    we have
472    \begin{equation}
473    \tau = \eta\hackscore{eff} \cdot \dot{\gamma}
474    \end{equation}
475    where
476    \begin{equation}
477    \eta\hackscore{eff}= min( \left(\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}\right)^{-1}
478    , \eta\hackscore{max}) \mbox{ with }
479    \eta\hackscore{max} = \left\{
480    \begin{array}{rcl}
481    \frac{\tau\hackscore{Y} + \beta \; p}{\dot{\gamma}} & & \dot{\gamma}>0 \\
482    &\mbox{ if } \\
483    \infty & & \mbox{otherwise}
484    \end{array}
485    \right.
486  \end{equation}  \end{equation}
487  where $\sigma\hackscore{ij}^{'-}$ is the deviatoric stress at the precious time step.  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
Now we can combine equations~\ref{IKM-EQU-2}, \ref{IKM-EQU-3b} and~\ref{IKM-EQU-6b} to get
488  \begin{equation}\label{IKM-EQU-10}  \begin{equation}\label{IKM-EQU-10}
489  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff}  \left( D\hackscore{ij}' +  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff}  \left( D\hackscore{ij}' +
490  \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}}
491  \end{equation}  \end{equation}
Notice that $\eta\hackscore{eff}$ is a function of diatoric stress $\sigma\hackscore{ij}'$.
492  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
493  \begin{equation}\label{IKM-EQU-1ib}  \begin{equation}\label{IKM-EQU-1ib}
494  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})
495  \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+  \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+
496  \frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij,j}^{'-}   \left(\frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij}^{'-} \right)\hackscore{,j}
497  \end{equation}  \end{equation}
498  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
499  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.
500  In oder to perform step~\ref{IKM iteration 2} we need to calculate the $\eta\hackscore{eff}$ which
501  is a function of $\sigma\hackscore{ij}$ via $\tau$.  To get $\tau$ and $\eta\hackscore{eff}$ we need to solve the  If we set
502  non-linear equation  \begin{equation}\label{IKM-EQU-44}
503  \begin{equation}  \frac{1}{\eta(\tau)}= \frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
504  \tau = \eta\hackscore{eff} \cdot \dot{\gamma}\hackscore{total} \mbox{ with }  \end{equation}
505  \dot{\gamma}\hackscore{total} = \sqrt{ 2 \left( D\hackscore{ij}' +  we need to solve the nonlinear problem
506  \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)^2}  \begin{equation}
507  \label{IKM iteration 5}  \eta\hackscore{eff} -  min(\eta( \dot{\gamma} \cdot \eta\hackscore{eff})
508  \end{equation}  , \eta\hackscore{max}) =0
509  The Newton scheme takes the form  \end{equation}
510  \begin{equation}  We use the Newton-Raphson Scheme \index{Newton-Raphson scheme} to solve this problem
511  \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}
512  \label{IKM iteration 6}  \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
513  \end{equation}  \eta\hackscore{eff}^{(n)} -
514  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)}) }
515    {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
516    =\min(\eta\hackscore{max},
517    \frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot  \eta'( \tau^{(n)} )  }
518    {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
519    \end{equation}
520    where $\eta'$ denotes the derivative of $\eta$ with respect of $\tau$
521    and $\tau^{(n)} = \dot{\gamma} \cdot \eta\hackscore{eff}^{(n)}$.
522
523    Looking at the evaluation of $\eta$ in~\ref{IKM-EQU-44} it makes sense formulate
524    the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$.
525    In fact we have
526  \begin{equation}  \begin{equation}
527  \eta\hackscore{eff}' = - \eta\hackscore{eff}^2 \left(\frac{1}{\eta\hackscore{eff}}\right)'  \eta' = - \frac{\Theta'}{\Theta^2}
528  \mbox{ with }  \mbox{ with }
529  \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)'
530  \label{IKM iteration 7}  \label{IKM iteration 7}
531  \end{equation}  \end{equation}
532  \begin{equation}\label{IKM-EQU-5XX}  As
533    \begin{equation}\label{IKM-EQU-47}
534  \left(\frac{1}{\eta^{q}} \right)'  \left(\frac{1}{\eta^{q}} \right)'
535  = \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{n^{q}-1}{\eta^{q}\hackscore{N}} \cdot \frac{\tau^{n^{q}-2}}{(\tau\hackscore{t}^q)^{n^{q}-1}}
536  = \frac{1-\frac{1}{n^{q}}}{\tau}\frac{1}{\eta^{q}}  = \frac{n^{q}-1}{\eta^{q}}\cdot\frac{1}{\tau}
537  \end{equation}  \end{equation}
538  Notice that allways $\eta\hackscore{eff}'\le 0$ which makes the denomionator in~\ref{IKM iteration 6}  we have
539  positive.  \begin{equation}\label{IKM-EQU-48}
540    \Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum\hackscore{q}\frac{n^{q}-1}{\eta^{q}}
541    \end{equation}
543    \begin{equation}\label{IKM-EQU-49}
544    \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
545    \eta\hackscore{eff}^{(n)}
546    \frac{\Theta^{(n)}  + \omega^{(n)}  }
547    {\eta\hackscore{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })
548    \end{equation}
549
550
551    \subsection{Functions}
552
553    \begin{classdesc}{IncompressibleIsotropicFlowCartesian}{
554    domain
555    \optional{, stress=0
556    \optional{, v=0
557    \optional{, p=0
558    \optional{, t=0
559    \optional{, numMaterials=1
560    \optional{, verbose=True
562    }}}}}}}}
563    opens an incompressible, isotropic flow problem in Cartesian cooridninates
564    on the domain \var{domain}.
565    \var{stress},
566    \var{v},
567    \var{p}, and
568    \var{t} set the initial deviatoric stress, velocity, pressure and time.
569    \var{numMaterials} specifies the number of materials used in the power law
570    model. Some progress information are printed if \var{verbose} is set to
571    \True. If \var{adaptSubTolerance} is equal to True the tolerances for subproblems are set automatically.
572    \end{classdesc}
573
574    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDomain}{}
575    returns the domain.
576    \end{methoddesc}
577
578    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTime}{}
579    Returns current time.
580    \end{methoddesc}
581
582    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getStress}{}
583    Returns current stress.
584    \end{methoddesc}
585
586    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStress}{}
587    Returns current deviatoric stress.
588    \end{methoddesc}
589
590    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getPressure}{}
591    Returns current pressure.
592    \end{methoddesc}
593
594    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getVelocity}{}
595    Returns current velocity.
596    \end{methoddesc}
597
598    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStrain}{}
599    Returns deviatoric strain of current velocity
600    \end{methoddesc}
601
602    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTau}{}
603    Returns current second invariant of deviatoric stress
604    \end{methoddesc}
605
607    Returns current second invariant of deviatoric strain
608    \end{methoddesc}
609
610
611    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setTolerance}{tol=1.e-4}
612    Sets the tolerance used to terminate the iteration on a time step.
613    \end{methoddesc}
614
615    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setFlowTolerance}{tol=1.e-4}
616    Sets the relative tolerance for the incompressible solver, see \class{StokesProblemCartesian} for details.
617    \end{methoddesc}
618
619    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
620    Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.
621    \end{methoddesc}
622
623    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setEtaTolerance=}{rtol=1.e-8}
624    sets the relative tolerance for the effectice viscosity. Iteration on a time step is completed if the realtive of the effective viscosity is less than \var{rtol}.
625    \end{methoddesc}
626
627    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setDruckerPragerLaw}
628    {\optional{tau_Y=None, \optional{friction=None}}}
629    Sets the parameters $\tau\hackscore{Y}$ and $\beta$ for the Drucker-Prager model in condition~\ref{IKM-EQU-8c}. If \var{tau_Y} is set to None (default) Drucker-Prager
630    condition is not applied.
631    \end{methoddesc}
632
633    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
634    Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.
635    \end{methoddesc}
636
637
638    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setPowerLaws}{eta_N, tau_t, power}
639    Sets the parameters of the power-law for all materials as defined in
640    equation~\ref{IKM-EQU-5b}.
641    \var{eta_N} is the list of viscosities $\eta^{q}\hackscore{N}$,
642    \var{tau_t} is the list of reference stresses  $\tau\hackscore{t}^q$,
643    and \var{power} is the list of  power law coefficients $n^{q}$.
644    \end{methoddesc}
645
646
647    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{update}{dt
648    \optional{, iter_max=100
649    \optional{, inner_iter_max=20
650    }}}
651    Updates stress, velocity and pressure for time increment \var{dt}.
652    where \var{iter_max} is the maximum number of iteration steps on a time step to
653    update the effective viscosity and \var{inner_iter_max} is the maximum
654    number of itertion steps in the incompressible solver.
655    \end{methoddesc}
656
657    \subsection{Example}
658    later
659
660    % \section{Drucker Prager Model}

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