/[escript]/trunk/doc/user/Models.tex
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2252 by gross, Fri Feb 6 07:51:28 2009 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
1    
2  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3  %  %
4  % Copyright (c) 2003-2008 by University of Queensland  % Copyright (c) 2003-2009 by University of Queensland
5  % Earth Systems Science Computational Center (ESSCC)  % Earth Systems Science Computational Center (ESSCC)
6  % http://www.uq.edu.au/esscc  % http://www.uq.edu.au/esscc
7  %  %
# 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    
131  \begin{classdesc}{StokesProblemCartesian}{domain}  \begin{classdesc}{StokesProblemCartesian}{domain\optional{, adaptSubTolerance=\True}}
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.
134    If \var{adaptSubTolerance} is \True
135    the tolerances for all subproblems are set automatically.
136  \end{classdesc}  \end{classdesc}
137    
138  \begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}}}}}}  \begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}}}}}}
# 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
151  by \var{fixed_u_mask} remain unchanged.  by \var{fixed_u_mask} remain unchanged.
# 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.
 On the otherhand is choosen to small compute time is wasted.  
 If \var{rtol} is set to \var{None} the sub-tolerance is set automatically depending on the  
 tolerance choosen for the oter iteration.  
178  \end{methoddesc}  \end{methoddesc}
179  \begin{methoddesc}[StokesProblemCartesian]{getSubProblemTolerance}{}  
180  return the tolerance for the involved PDEs.  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{}
181    returns the solver options used  solve the equation~(\ref{P PREC}) for pressure.
182  \end{methoddesc}  \end{methoddesc}
183    
184    \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}
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}
344  \begin{classdesc}{DarcyFlow}{domain}  \begin{classdesc}{DarcyFlow}{domain \optional{, adaptSubTolerance=\True}}
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}
 % \section{Temperature Advection Diffusion\label{TEMP ADV DIFF}}  
   
 %\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}
542    which leads to
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
561    \optional{, adaptSubTolerance=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    
606    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getGammaDot}{}
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.2548

  ViewVC Help
Powered by ViewVC 1.1.26