doc/user/Models.tex


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2003-2010 by University of Queensland
% Earth Systems Science Computational Center (ESSCC)
% http://www.uq.edu.au/esscc
%
% Primary Business: Queensland, Australia
% Licensed under the Open Software License version 3.0
% http://www.opensource.org/licenses/osl-3.0.php
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Models}
\label{MODELS CHAPTER}

The following sections give a brief overview of the model classes and their
corresponding methods.

\input{stokessolver} 
%\input{darcyfluxNew}
\input{darcyflux}
%\input{levelsetmodel}

\section{Isotropic Kelvin Material \label{IKM}}
As proposed by Kelvin~\cite{Muhlhaus2005} material strain $D_{ij}=\frac{1}{2}(v_{i,j}+v_{j,i})$ can be decomposed into
an elastic part $D_{ij}^{el}$ and visco-plastic part $D_{ij}^{vp}$:
\begin{equation}\label{IKM-EQU-2}
D_{ij}=D_{ij}^{el}+D_{ij}^{vp}
\end{equation}
with the elastic strain given as 
\begin{equation}\label{IKM-EQU-3}
D_{ij}^{el'}=\frac{1}{2 \mu} \dot{\sigma}_{ij}'
\end{equation}
where $\sigma'_{ij}$ is the deviatoric stress (Notice that $\sigma'_{ii}=0$).
If the material is composed by materials $q$ the visco-plastic strain can be decomposed as
\begin{equation}\label{IKM-EQU-4}
D_{ij}^{vp'}=\sum_{q} D_{ij}^{q'} 
\end{equation}
where $D_{ij}^{q}$ is the strain in material $q$ given as 
\begin{equation}\label{IKM-EQU-5}
D_{ij}^{q'}=\frac{1}{2 \eta^{q}} \sigma'_{ij} 
\end{equation}
where $\eta^{q}$ is the viscosity of material $q$. We assume the following 
betwee the the strain in material $q$ 
\begin{equation}\label{IKM-EQU-5b}
\eta^{q}=\eta^{q}_{N} \left(\frac{\tau}{\tau_{t}^q}\right)^{1-n^{q}}
\mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'_{ij} \sigma'_{ij}}
\end{equation}
for a given power law coefficients $n^{q}\ge1$ and transition stresses $\tau_{t}^q$, see~\cite{Muhlhaus2005}.
Notice that $n^{q}=1$ gives a constant viscosity.
After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:
\begin{equation}\label{IKM-EQU-6}
D_{ij}'^{vp}=\frac{1}{2 \eta^{vp}} \sigma'_{ij} \mbox{ with } \frac{1}{\eta^{vp}} = \sum_{q} \frac{1}{\eta^{q}} \;.
\end{equation}
and finally with~\ref{IKM-EQU-2} 
\begin{equation}\label{IKM-EQU-2bb}
D_{ij}'=\frac{1}{2 \eta^{vp}} \sigma'_{ij}+\frac{1}{2 \mu} \dot{\sigma}_{ij}'
\end{equation}
The total stress $\tau$ needs to fullfill the yield condition \index{yield condition} 
\begin{equation}\label{IKM-EQU-8c}
\tau \le \tau_{Y} + \beta \; p
\end{equation}
with the Drucker-Prager \index{Druck-Prager} cohesion factor \index{cohesion factor} $\tau_{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$.
The deviatoric stress needs to fullfill the equilibrion equation
\begin{equation}\label{IKM-EQU-1}
-\sigma'_{ij,j}+p_{,i}=F_{i}
\end{equation}
where $F_{j}$ is a given external fource. We assume an incompressible media:
\begin{equation}\label{IKM-EQU-2bbb}
-v_{i,i}=0
\end{equation}
Natural boundary conditions are taken in the form 
\begin{equation}\label{IKM-EQU-Boundary}
\sigma'_{ij}n_{j}-n_{i}p=f
\end{equation}
which can be overwritten by a constraint 
\begin{equation}\label{IKM-EQU-Boundary0}
v_{i}(x)=0
\end{equation}
where the index $i$ may depend on the location $x$ on the bondary.

\subsection{Solution Method \label{IKM-SOLVE}}
By using a first order finite difference approximation wit step size $dt>0$~\ref{IKM-EQU-3} get the form
\begin{equation}\label{IKM-EQU-3b}
\dot{\sigma}_{ij}=\frac{1}{dt } \left( \sigma_{ij} - \sigma_{ij}^{-} \right)
\end{equation}
and 
\begin{equation}\label{IKM-EQU-2b}
D_{ij}'=\left(\frac{1}{2 \eta^{vp}} + \frac{1}{2 \mu dt}\right) \sigma_{ij}'-\frac{1}{2 \mu dt } \sigma_{ij}^{-'}
\end{equation}
where $\sigma_{ij}^{-}$ is the stress at the precious time step. With 
\begin{equation}\label{IKM-EQU-2c}
\dot{\gamma} = \sqrt{ 2 \left( D_{ij}' +
\frac{1}{  2 \mu \; dt} \sigma_{ij}^{-'}\right)^2}
\end{equation} 
we have
\begin{equation}
\tau = \eta_{eff} \cdot \dot{\gamma}
\end{equation} 
where
\begin{equation}
\eta_{eff}= min( \left(\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}\right)^{-1} 
, \eta_{max}) \mbox{ with } 
\eta_{max} = \left\{ 
\begin{array}{rcl}
\frac{\tau_{Y} + \beta \; p}{\dot{\gamma}} & & \dot{\gamma}>0 \\
&\mbox{ if } \\ 
\infty & & \mbox{otherwise}
\end{array}
\right.
\end{equation}
The upper bound $\eta_{max}$ makes sure that yield condtion~\ref{IKM-EQU-8c} holds. With this setting the eqaution \ref{IKM-EQU-2b} takes the form 
\begin{equation}\label{IKM-EQU-10}
\sigma_{ij}' =  2 \eta_{eff}  \left( D_{ij}' +
\frac{1}{  2 \mu \; dt} \sigma_{ij}^{'-}\right)  
\end{equation}
After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get
\begin{equation}\label{IKM-EQU-1ib}
-\left(\eta_{eff} (v_{i,j}+ v_{i,j})
\right)_{,j}+p_{,i}=F_{i}+
 \left(\frac{\eta_{eff}}{\mu dt } \sigma_{ij}^{'-} \right)_{,j}
\end{equation}
Combining this with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a
Stokes problem as discussed in section~\ref{STOKES SOLVE} in each time step.

If we set 
\begin{equation}\label{IKM-EQU-44}
\frac{1}{\eta(\tau)}= \frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
\end{equation}
we need to solve the nonlinear problem
\begin{equation}
\eta_{eff} -  min(\eta( \dot{\gamma} \cdot \eta_{eff}) 
, \eta_{max}) =0 
\end{equation}
We use the Newton-Raphson Scheme \index{Newton-Raphson scheme} to solve this problem
\begin{equation}\label{IKM-EQU-45}
\eta_{eff}^{(n+1)} = \min(\eta_{max}, 
\eta_{eff}^{(n)} -
\frac{\eta_{eff}^{(n)} - \eta( \tau^{(n)}) }
{1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
=\min(\eta_{max},
\frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot  \eta'( \tau^{(n)} )  }
{1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
\end{equation} 
where $\eta'$ denotes the derivative of $\eta$ with respect of $\tau$
and $\tau^{(n)} =  \dot{\gamma} \cdot \eta_{eff}^{(n)}$. 

Looking at the evaluation of $\eta$ in~\ref{IKM-EQU-44} it makes sense formulate
the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$. 
In fact we have 
\begin{equation}
\eta' = - \frac{\Theta'}{\Theta^2} 
\mbox{ with } 
\Theta' = \sum_{q} \left(\frac{1}{\eta^{q}} \right)'
\label{IKM iteration 7}
\end{equation} 
As
\begin{equation}\label{IKM-EQU-47}
\left(\frac{1}{\eta^{q}} \right)'
= \frac{n^{q}-1}{\eta^{q}_{N}} \cdot \frac{\tau^{n^{q}-2}}{(\tau_{t}^q)^{n^{q}-1}}
= \frac{n^{q}-1}{\eta^{q}}\cdot\frac{1}{\tau} 
\end{equation}
we have
\begin{equation}\label{IKM-EQU-48}
\Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum_{q}\frac{n^{q}-1}{\eta^{q}} 
\end{equation}
which leads to 
\begin{equation}\label{IKM-EQU-49}
\eta_{eff}^{(n+1)} = \min(\eta_{max}, 
\eta_{eff}^{(n)}
\frac{\Theta^{(n)}  + \omega^{(n)}  }
{\eta_{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })
\end{equation} 


\subsection{Functions}

\begin{classdesc}{IncompressibleIsotropicFlowCartesian}{
domain
\optional{, stress=0
\optional{, v=0
\optional{, p=0
\optional{, t=0
\optional{, numMaterials=1
\optional{, verbose=True
\optional{, adaptSubTolerance=True
}}}}}}}}
opens an incompressible, isotropic flow problem in Cartesian cooridninates
on the domain \var{domain}.
\var{stress},
\var{v},
\var{p}, and
\var{t} set the initial deviatoric stress, velocity, pressure and time.
\var{numMaterials} specifies the number of materials used in the power law
model. Some progress information are printed if \var{verbose} is set to
\True. If \var{adaptSubTolerance} is equal to True the tolerances for subproblems are set automatically.

The domain
needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for detials~\index{LBB condition}.
For instance one can use second order polynomials for velocity and 
first order polynomials for the pressure on the same element. Alternativly, one can use 
macro elements\index{macro elements} using linear polynomial for both pressure and velocity bu with a subdivided
element for the velocity. Typically, the macro element is more cost effective. The fact that pressure and velocity are represented in different way is expressed by
\begin{python}
velocity=Vector(0.0, Solution(mesh))
pressure=Scalar(0.0, ReducedSolution(mesh))
\end{python}
\end{classdesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDomain}{}
returns the domain.
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTime}{}
Returns current time.
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getStress}{}
Returns current stress.
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStress}{}
Returns current deviatoric stress.
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getPressure}{}
Returns current pressure.
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getVelocity}{}
Returns current velocity.
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStrain}{}
Returns deviatoric strain of current velocity
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTau}{}
Returns current second invariant of deviatoric stress
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getGammaDot}{}
Returns current second invariant of deviatoric strain
\end{methoddesc}


\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setTolerance}{tol=1.e-4}
Sets the tolerance used to terminate the iteration on a time step.
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setFlowTolerance}{tol=1.e-4}
Sets the relative tolerance for the incompressible solver, see \class{StokesProblemCartesian} for details.
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setEtaTolerance=}{rtol=1.e-8}
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}.
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setDruckerPragerLaw}
{\optional{tau_Y=None, \optional{friction=None}}}
Sets the parameters $\tau_{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
condition is not applied.
\end{methoddesc}

\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.
\end{methoddesc}


\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setPowerLaws}{eta_N, tau_t, power}
Sets the parameters of the power-law for all materials as defined in 
equation~\ref{IKM-EQU-5b}.
\var{eta_N} is the list of viscosities $\eta^{q}_{N}$,
\var{tau_t} is the list of reference stresses  $\tau_{t}^q$,
and \var{power} is the list of  power law coefficients $n^{q}$.
\end{methoddesc}


\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{update}{dt
\optional{, iter_max=100
\optional{, inner_iter_max=20
}}}
Updates stress, velocity and pressure for time increment \var{dt}.
where \var{iter_max} is the maximum number of iteration steps on a time step to
update the effective viscosity and \var{inner_iter_max} is the maximum
number of itertion steps in the incompressible solver.
\end{methoddesc}

\subsection{Example}
later

\input{faultsystem}

% \section{Drucker Prager Model}
1	ksteube	1811
2			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3	lgraham	1702	%
4	jfenwick	2881	% Copyright (c) 2003-2010 by University of Queensland
5	ksteube	1811	% Earth Systems Science Computational Center (ESSCC)
6			% http://www.uq.edu.au/esscc
7	lgraham	1702	%
8	ksteube	1811	% Primary Business: Queensland, Australia
9			% Licensed under the Open Software License version 3.0
10			% http://www.opensource.org/licenses/osl-3.0.php
11	lgraham	1702	%
12	ksteube	1811	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13	lgraham	1702
14			\chapter{Models}
15	lgraham	2128	\label{MODELS CHAPTER}
16	lgraham	1702
17	caltinay	3325	The following sections give a brief overview of the model classes and their
18			corresponding methods.
19	lgraham	1702
20	gross	2719	\input{stokessolver}
21	gross	3072	%\input{darcyfluxNew}
22			\input{darcyflux}
23	gross	2432	%\input{levelsetmodel}
24	gross	2100
25	gross	1841	\section{Isotropic Kelvin Material \label{IKM}}
26	jfenwick	3295	As proposed by Kelvin~\cite{Muhlhaus2005} material strain $D_{ij}=\frac{1}{2}(v_{i,j}+v_{j,i})$ can be decomposed into
27			an elastic part $D_{ij}^{el}$ and visco-plastic part $D_{ij}^{vp}$:
28	gross	1841	\begin{equation}\label{IKM-EQU-2}
29	jfenwick	3295	D_{ij}=D_{ij}^{el}+D_{ij}^{vp}
30	gross	1841	\end{equation}
31	gross	1859	with the elastic strain given as
32	gross	1841	\begin{equation}\label{IKM-EQU-3}
33	jfenwick	3295	D_{ij}^{el'}=\frac{1}{2 \mu} \dot{\sigma}_{ij}'
34	gross	1841	\end{equation}
35	jfenwick	3295	where $\sigma'_{ij}$ is the deviatoric stress (Notice that $\sigma'_{ii}=0$).
36	gross	1859	If the material is composed by materials $q$ the visco-plastic strain can be decomposed as
37	gross	1841	\begin{equation}\label{IKM-EQU-4}
38	jfenwick	3295	D_{ij}^{vp'}=\sum_{q} D_{ij}^{q'}
39	gross	1841	\end{equation}
40	jfenwick	3295	where $D_{ij}^{q}$ is the strain in material $q$ given as
41	gross	1841	\begin{equation}\label{IKM-EQU-5}
42	jfenwick	3295	D_{ij}^{q'}=\frac{1}{2 \eta^{q}} \sigma'_{ij}
43	gross	1841	\end{equation}
44	gross	1859	where $\eta^{q}$ is the viscosity of material $q$. We assume the following
45			betwee the the strain in material $q$
46			\begin{equation}\label{IKM-EQU-5b}
47	jfenwick	3295	\eta^{q}=\eta^{q}_{N} \left(\frac{\tau}{\tau_{t}^q}\right)^{1-n^{q}}
48			\mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'_{ij} \sigma'_{ij}}
49	gross	1859	\end{equation}
50	jfenwick	3295	for a given power law coefficients $n^{q}\ge1$ and transition stresses $\tau_{t}^q$, see~\cite{Muhlhaus2005}.
51	gross	1859	Notice that $n^{q}=1$ gives a constant viscosity.
52	gross	1841	After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:
53	gross	1859	\begin{equation}\label{IKM-EQU-6}
54	jfenwick	3295	D_{ij}'^{vp}=\frac{1}{2 \eta^{vp}} \sigma'_{ij} \mbox{ with } \frac{1}{\eta^{vp}} = \sum_{q} \frac{1}{\eta^{q}} \;.
55	gross	1841	\end{equation}
56	gross	2300	and finally with~\ref{IKM-EQU-2}
57	gross	2371	\begin{equation}\label{IKM-EQU-2bb}
58	jfenwick	3295	D_{ij}'=\frac{1}{2 \eta^{vp}} \sigma'_{ij}+\frac{1}{2 \mu} \dot{\sigma}_{ij}'
59	gross	1859	\end{equation}
60	gross	2300	The total stress $\tau$ needs to fullfill the yield condition \index{yield condition}
61	gross	1859	\begin{equation}\label{IKM-EQU-8c}
62	jfenwick	3295	\tau \le \tau_{Y} + \beta \; p
63	gross	1859	\end{equation}
64	jfenwick	3295	with the Drucker-Prager \index{Druck-Prager} cohesion factor \index{cohesion factor} $\tau_{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$.
65	gross	1859	The deviatoric stress needs to fullfill the equilibrion equation
66	gross	1841	\begin{equation}\label{IKM-EQU-1}
67	jfenwick	3295	-\sigma'_{ij,j}+p_{,i}=F_{i}
68	gross	1841	\end{equation}
69	jfenwick	3295	where $F_{j}$ is a given external fource. We assume an incompressible media:
70	gross	2371	\begin{equation}\label{IKM-EQU-2bbb}
71	jfenwick	3295	-v_{i,i}=0
72	gross	1841	\end{equation}
73	gross	1878	Natural boundary conditions are taken in the form
74			\begin{equation}\label{IKM-EQU-Boundary}
75	jfenwick	3295	\sigma'_{ij}n_{j}-n_{i}p=f
76	gross	1878	\end{equation}
77			which can be overwritten by a constraint
78			\begin{equation}\label{IKM-EQU-Boundary0}
79	jfenwick	3295	v_{i}(x)=0
80	gross	1878	\end{equation}
81			where the index $i$ may depend on the location $x$ on the bondary.
82
83			\subsection{Solution Method \label{IKM-SOLVE}}
84			By using a first order finite difference approximation wit step size $dt>0$~\ref{IKM-EQU-3} get the form
85			\begin{equation}\label{IKM-EQU-3b}
86	jfenwick	3295	\dot{\sigma}_{ij}=\frac{1}{dt } \left( \sigma_{ij} - \sigma_{ij}^{-} \right)
87	gross	1878	\end{equation}
88	gross	2300	and
89			\begin{equation}\label{IKM-EQU-2b}
90	jfenwick	3295	D_{ij}'=\left(\frac{1}{2 \eta^{vp}} + \frac{1}{2 \mu dt}\right) \sigma_{ij}'-\frac{1}{2 \mu dt } \sigma_{ij}^{-'}
91	gross	2300	\end{equation}
92	jfenwick	3295	where $\sigma_{ij}^{-}$ is the stress at the precious time step. With
93	gross	2300	\begin{equation}\label{IKM-EQU-2c}
94	jfenwick	3295	\dot{\gamma} = \sqrt{ 2 \left( D_{ij}' +
95			\frac{1}{ 2 \mu \; dt} \sigma_{ij}^{-'}\right)^2}
96	gross	2300	\end{equation}
97			we have
98			\begin{equation}
99	jfenwick	3295	\tau = \eta_{eff} \cdot \dot{\gamma}
100	gross	2300	\end{equation}
101			where
102			\begin{equation}
103	jfenwick	3295	\eta_{eff}= min( \left(\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}\right)^{-1}
104			, \eta_{max}) \mbox{ with }
105			\eta_{max} = \left\{
106	gross	2300	\begin{array}{rcl}
107	jfenwick	3295	\frac{\tau_{Y} + \beta \; p}{\dot{\gamma}} & & \dot{\gamma}>0 \\
108	gross	2300	&\mbox{ if } \\
109			\infty & & \mbox{otherwise}
110			\end{array}
111			\right.
112			\end{equation}
113	jfenwick	3295	The upper bound $\eta_{max}$ makes sure that yield condtion~\ref{IKM-EQU-8c} holds. With this setting the eqaution \ref{IKM-EQU-2b} takes the form
114	gross	1878	\begin{equation}\label{IKM-EQU-10}
115	jfenwick	3295	\sigma_{ij}' = 2 \eta_{eff} \left( D_{ij}' +
116			\frac{1}{ 2 \mu \; dt} \sigma_{ij}^{'-}\right)
117	gross	1878	\end{equation}
118	gross	1859	After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get
119			\begin{equation}\label{IKM-EQU-1ib}
120	jfenwick	3295	-\left(\eta_{eff} (v_{i,j}+ v_{i,j})
121			\right)_{,j}+p_{,i}=F_{i}+
122			\left(\frac{\eta_{eff}}{\mu dt } \sigma_{ij}^{'-} \right)_{,j}
123	gross	1859	\end{equation}
124	gross	2252	Combining this with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a
125			Stokes problem as discussed in section~\ref{STOKES SOLVE} in each time step.
126	gross	2300
127			If we set
128			\begin{equation}\label{IKM-EQU-44}
129			\frac{1}{\eta(\tau)}= \frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
130			\end{equation}
131			we need to solve the nonlinear problem
132	gross	1878	\begin{equation}
133	jfenwick	3295	\eta_{eff} - min(\eta( \dot{\gamma} \cdot \eta_{eff})
134			, \eta_{max}) =0
135	gross	2300	\end{equation}
136			We use the Newton-Raphson Scheme \index{Newton-Raphson scheme} to solve this problem
137			\begin{equation}\label{IKM-EQU-45}
138	jfenwick	3295	\eta_{eff}^{(n+1)} = \min(\eta_{max},
139			\eta_{eff}^{(n)} -
140			\frac{\eta_{eff}^{(n)} - \eta( \tau^{(n)}) }
141	gross	2300	{1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
142	jfenwick	3295	=\min(\eta_{max},
143	gross	2300	\frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot \eta'( \tau^{(n)} ) }
144			{1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
145	gross	2100	\end{equation}
146	gross	2300	where $\eta'$ denotes the derivative of $\eta$ with respect of $\tau$
147	jfenwick	3295	and $\tau^{(n)} = \dot{\gamma} \cdot \eta_{eff}^{(n)}$.
148	gross	2300
149			Looking at the evaluation of $\eta$ in~\ref{IKM-EQU-44} it makes sense formulate
150			the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$.
151			In fact we have
152	gross	2100	\begin{equation}
153	gross	2300	\eta' = - \frac{\Theta'}{\Theta^2}
154	gross	2100	\mbox{ with }
155	jfenwick	3295	\Theta' = \sum_{q} \left(\frac{1}{\eta^{q}} \right)'
156	gross	2100	\label{IKM iteration 7}
157			\end{equation}
158	gross	2300	As
159			\begin{equation}\label{IKM-EQU-47}
160	gross	2100	\left(\frac{1}{\eta^{q}} \right)'
161	jfenwick	3295	= \frac{n^{q}-1}{\eta^{q}_{N}} \cdot \frac{\tau^{n^{q}-2}}{(\tau_{t}^q)^{n^{q}-1}}
162	gross	2438	= \frac{n^{q}-1}{\eta^{q}}\cdot\frac{1}{\tau}
163	gross	2100	\end{equation}
164	gross	2300	we have
165			\begin{equation}\label{IKM-EQU-48}
166	jfenwick	3295	\Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum_{q}\frac{n^{q}-1}{\eta^{q}}
167	gross	2300	\end{equation}
168			which leads to
169	gross	2371	\begin{equation}\label{IKM-EQU-49}
170	jfenwick	3295	\eta_{eff}^{(n+1)} = \min(\eta_{max},
171			\eta_{eff}^{(n)}
172	gross	2300	\frac{\Theta^{(n)} + \omega^{(n)} }
173	jfenwick	3295	{\eta_{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })
174	gross	2432	\end{equation}
175
176
177			\subsection{Functions}
178
179			\begin{classdesc}{IncompressibleIsotropicFlowCartesian}{
180			domain
181			\optional{, stress=0
182			\optional{, v=0
183			\optional{, p=0
184			\optional{, t=0
185			\optional{, numMaterials=1
186	gross	2474	\optional{, verbose=True
187			\optional{, adaptSubTolerance=True
188			}}}}}}}}
189	gross	2433	opens an incompressible, isotropic flow problem in Cartesian cooridninates
190			on the domain \var{domain}.
191	gross	2432	\var{stress},
192			\var{v},
193			\var{p}, and
194			\var{t} set the initial deviatoric stress, velocity, pressure and time.
195			\var{numMaterials} specifies the number of materials used in the power law
196			model. Some progress information are printed if \var{verbose} is set to
197	gross	2474	\True. If \var{adaptSubTolerance} is equal to True the tolerances for subproblems are set automatically.
198	gross	2793
199			The domain
200			needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for detials~\index{LBB condition}.
201			For instance one can use second order polynomials for velocity and
202			first order polynomials for the pressure on the same element. Alternativly, one can use
203			macro elements\index{macro elements} using linear polynomial for both pressure and velocity bu with a subdivided
204			element for the velocity. Typically, the macro element is more cost effective. The fact that pressure and velocity are represented in different way is expressed by
205			\begin{python}
206			velocity=Vector(0.0, Solution(mesh))
207			pressure=Scalar(0.0, ReducedSolution(mesh))
208			\end{python}
209	gross	2432	\end{classdesc}
210
211	gross	2433	\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDomain}{}
212			returns the domain.
213			\end{methoddesc}
214	gross	2432
215	gross	2433	\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTime}{}
216			Returns current time.
217	gross	2432	\end{methoddesc}
218
219	gross	2433	\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getStress}{}
220			Returns current stress.
221	gross	2432	\end{methoddesc}
222
223	gross	2433	\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStress}{}
224			Returns current deviatoric stress.
225			\end{methoddesc}
226	gross	2432
227	gross	2433	\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getPressure}{}
228			Returns current pressure.
229	gross	2432	\end{methoddesc}
230	gross	2433
231			\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getVelocity}{}
232			Returns current velocity.
233	gross	2432	\end{methoddesc}
234	gross	2433
235			\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStrain}{}
236			Returns deviatoric strain of current velocity
237	gross	2432	\end{methoddesc}
238	gross	2433
239			\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTau}{}
240			Returns current second invariant of deviatoric stress
241	gross	2432	\end{methoddesc}
242	gross	2433
243			\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getGammaDot}{}
244			Returns current second invariant of deviatoric strain
245	gross	2432	\end{methoddesc}
246	gross	2433
247
248			\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setTolerance}{tol=1.e-4}
249			Sets the tolerance used to terminate the iteration on a time step.
250	gross	2432	\end{methoddesc}
251
252	gross	2433	\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setFlowTolerance}{tol=1.e-4}
253			Sets the relative tolerance for the incompressible solver, see \class{StokesProblemCartesian} for details.
254			\end{methoddesc}
255	gross	2432
256	gross	2474	\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
257	gross	2433	Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.
258			\end{methoddesc}
259
260			\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setEtaTolerance=}{rtol=1.e-8}
261			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}.
262			\end{methoddesc}
263
264			\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setDruckerPragerLaw}
265			{\optional{tau_Y=None, \optional{friction=None}}}
266	jfenwick	3295	Sets the parameters $\tau_{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
267	gross	2433	condition is not applied.
268			\end{methoddesc}
269
270			\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
271			Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.
272			\end{methoddesc}
273
274
275			\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setPowerLaws}{eta_N, tau_t, power}
276			Sets the parameters of the power-law for all materials as defined in
277			equation~\ref{IKM-EQU-5b}.
278	jfenwick	3295	\var{eta_N} is the list of viscosities $\eta^{q}_{N}$,
279			\var{tau_t} is the list of reference stresses $\tau_{t}^q$,
280	gross	2433	and \var{power} is the list of power law coefficients $n^{q}$.
281			\end{methoddesc}
282
283
284			\begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{update}{dt
285			\optional{, iter_max=100
286			\optional{, inner_iter_max=20
287			}}}
288			Updates stress, velocity and pressure for time increment \var{dt}.
289			where \var{iter_max} is the maximum number of iteration steps on a time step to
290			update the effective viscosity and \var{inner_iter_max} is the maximum
291			number of itertion steps in the incompressible solver.
292			\end{methoddesc}
293
294			\subsection{Example}
295			later
296
297	gross	2647	\input{faultsystem}
298
299	gross	3072	% \section{Drucker Prager Model}