/[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 2432 by gross, Wed May 20 06:06:20 2009 UTC
# Line 183  return the tolerance for the involved PD Line 183  return the tolerance for the involved PD
183   The following script \file{lit\hackscore driven\hackscore cavity.py}   The following script \file{lit\hackscore driven\hackscore cavity.py}
184  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
185  illustrates the usage of the \class{StokesProblemCartesian} class to solve  illustrates the usage of the \class{StokesProblemCartesian} class to solve
186  the lit driven cavity problem~\cite{LITDRIVENCAVITY}:  the lit driven cavity problem:
187  \begin{python}  \begin{python}
188  from esys.escript import *  from esys.escript import *
189  from esys.finley import Rectangle  from esys.finley import Rectangle
# Line 281  Q^*u  + Q^*Q p & = & Q^*g \\ Line 281  Q^*u  + Q^*Q p & = & Q^*g \\
281  \end{array}  \end{array}
282  \end{equation}  \end{equation}
283  where $D^*$ and $Q^*$ denote the adjoint operators.  where $D^*$ and $Q^*$ denote the adjoint operators.
284  In~\cite{XXX} it has been shown that this problem is continuous and coercive in $W \times V$ and therefore has a unique solution. Also standart FEM methods can be used for discretization. It is also possible  In~\cite{LEASTSQUARESFEM1994} it has been shown that this problem is continuous and coercive in $W \times V$ and therefore has a unique solution. Also standart FEM methods can be used for discretization. It is also possible
285  to solve the problem is coupled form, however this approach leads in some cases to a very ill-conditioned stiffness matrix in particular in the case of a very small or large permability ($\alpha\hackscore{1} \ll 1$ or $\alpha\hackscore{0} \gg 1$)    to solve the problem is coupled form, however this approach leads in some cases to a very ill-conditioned stiffness matrix in particular in the case of a very small or large permability ($\alpha\hackscore{1} \ll 1$ or $\alpha\hackscore{0} \gg 1$)  
286    
287  The approach we are taking is to eliminate the velocity $u$ from the problem. Assuming that $p$ is known we have  The approach we are taking is to eliminate the velocity $u$ from the problem. Assuming that $p$ is known we have
288  \begin{equation}  \begin{equation}\label{DARCY V FORM}
289  v= (I+D^*D)^{-1} (D^*f + g - Qp)  v= (I+D^*D)^{-1} (D^*f + g - Qp)
290  \end{equation}  \end{equation}
291  (notice that $(I+D^*D)$ is coercive in $W$) which is inserted into the second equation  (notice that $(I+D^*D)$ is coercive in $W$) which is inserted into the second equation
# Line 300  We use the PCG \index{linear solver!PCG} Line 300  We use the PCG \index{linear solver!PCG}
300  \begin{equation}  \begin{equation}
301  \begin{array}{rcl}  \begin{array}{rcl}
302  r & = & Q^*  \left( g -(I+D^*D)^{-1} (D^*f + g) - Qp + (I+D^*D)^{-1}Q p \right)\\  r & = & Q^*  \left( g -(I+D^*D)^{-1} (D^*f + g) - Qp + (I+D^*D)^{-1}Q p \right)\\
303  & =&  Q^* \left( - Qp - (I+D^*D)^{-1} (D^*f + g - Qp) \right) \\  & =&  Q^* \left( g- Qp - (I+D^*D)^{-1} (D^*f + g - Qp) \right) \\
304  & =&  Q^* \left( g - Qp - v \right)  & =&  Q^* \left( g - Qp - v \right)
305  \end{array}  \end{array}
306  \end{equation}  \end{equation}
307  So in a partical implementation we use the pair $(Qp,v)$ to represent the residual. This will save the  So in a partical implementation we use $\hat{r}=g-Qp-v$ to represent the residual.
308  reconstruction of the velocity $v$. In this notation the right hand side is given as  The evaluation of the iteration operator for a given $p$ is then
309  $(0,(I+D^*D)^{-1} (D^*f + g))$. The evaluation of the iteration operator for a given $p$ is then  returning $Qp+v$ where $v$ is the solution of
 returning $(Qp,w)$ where $w$ is the solution of  
310  \begin{equation}\label{UPDATE W}  \begin{equation}\label{UPDATE W}
311  (I+D^*D)w = Qp  (I+D^*D)v = Qp
312  \end{equation}  \end{equation}
313  We use $Q^*Q$ as a a preconditioner for the iteration operator $Q^* ( I - (I+D^*D)^{-1} ) Q$.  We use $(Q^*Q)^{-1}$ as a preconditioner for the iteration operator $Q^* ( I - (I+D^*D)^{-1} ) Q$. So the application of the preconditioner to $\hat{r}$ representing the residual is given by solving
314    implemented by solving
315  The iteration PCG \index{linear solver!PCG}\index{PCG} is terminated if  \begin{equation}\label{UPDATE P}
316    Q^*Q q  = Q^*\hat{r}
317    \end{equation}
318    The residual norm used in the PCG is given as
319    \begin{equation}\label{DARCY R NORM}
320    \|r\|\hackscore{PCG}^2 = \int r \cdot (Q^*Q)^{-1} r \; dx =\int \hat{r} \cdot Q (Q^*Q)^{-1} Q^* \hat{r} \; dx \approx
321    \|\hat{r}\|\hackscore{0}^2
322    \end{equation}
323    The iteration is terminated if
324  \begin{equation}\label{DARCY STOP}  \begin{equation}\label{DARCY STOP}
325  \int r \cdot (Q^*Q)^{-1} r \; dx \le \mbox{ATOL}^2  \|r\|\hackscore{PCG} \le \mbox{ATOL}
326  \end{equation}  \end{equation}
327  where ATOL is a given absolute tolerance.  where we set
328  The initial residual $r\hackscore{0}$ is  \begin{equation}\label{DARCY ATOL DEF}
329  \begin{equation}\label{DARCY STOP 2}  \mbox{ATOL} = \mbox{atol} + \mbox{rtol} \cdot \left(\frac{1}{\|v\|\hackscore{0}} + \frac{1}{\|Qp\|\hackscore{0}} \right)^{-1}
330  r\hackscore{0}=Q^* \left( g - v\hackscore{ref} \right) \mbox{ with } v\hackscore{ref} = (I+D^*D)^{-1} (D^*f + g)  \end{equation}
331  \end{equation}  where rtol is a given relative tolerance and $\mbox{atol}$ is a given absolute tolerance (typically $=0$).  
332  so the  Notice that if $Qp$ and $v$ both are zero, the pair $(0,p)$ is a solution.
333  \begin{equation}\label{DARCY NORM 0}  The problem is that ATOL is depending on the solution $p$ (and $v$ calculated form~\ref{DARCY V FORM}). In partcice one use the initial guess for $p$
334  \int r\hackscore0 \cdot (Q^*Q)^{-1} r\hackscore0 \; dx = \int \left( g - v\hackscore{ref} \right)  \cdot  Q  p\hackscore{ref} \; dx \mbox{ with }p\hackscore{ref} = (Q^*Q)^{-1} Q^* \left( g - v\hackscore{ref} \right)  to get a first value for ATOL. If the stopping crierion is met in the PCG iteration, a new $v$ is calculated from the current pressure approximation and ATOL is recalculated. If \ref{DARCY STOP} is still fullfilled the calculation is terminated and $(v,p)$ is returned. Otherwise PCG is restarted with a new ATOL.
 \end{equation}  
 So we set  
 \begin{equation}\label{DARCY NORM 1}  
 ATOL = atol + rtol \cdot \max(|g - v\hackscore{ref}|\hackscore{0}, |Q p\hackscore{ref} |\hackscore{0} )  
 \end{equation}  
 where atol and rtol a given absolute and relative tolerances, respectively. The reference flux $v\hackscore{ref}$  
 and reference pressure $p\hackscore{ref}$ may be calcualated from their definition which would require to solve to  
 PDEs but in a practical application estimates can be used for instance solutions from previous time steps or for simplified scenarious (e.g. constant permability).  
335    
336  \subsection{Functions}  \subsection{Functions}
337  \begin{classdesc}{DarcyFlow}{domain}  \begin{classdesc}{DarcyFlow}{domain}
# Line 353  The method will try to cast the given va Line 352  The method will try to cast the given va
352  \Data class objects.  \Data class objects.
353  \end{methoddesc}  \end{methoddesc}
354    
355  \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{atol=0,\optional{rtol=1e-8,\optional{p_ref=None,\optional{v_ref=None}}}}}  \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e-4}}
356  sets the absolute tolerance ATOL according to~\ref{DARCY NORM 1}. If \var{p_ref} is not present $0$ is used.  sets the relative tolerance \mbox{rtol} in \ref{DARCY ATOL DEF}.
 If \var{v_ref} is not present $0$ is used. If the final result ATOL is not positive an exception is thrown.  
357  \end{methoddesc}  \end{methoddesc}
358    
359    \begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}}
360    sets the absolute tolerance \mbox{atol} in \ref{DARCY ATOL DEF}.
361    \end{methoddesc}
362    
363    \begin{methoddesc}[DarcyFlow]{setSubProblemTolerance}{\optional{rtol=None}}
364    sets the relative tolerance used to solve the involved PDEs. If no argument is given,
365    the square of the current relative tolerance is used. The sub-problem tolerance should be choosen as large as possible to minimize the compute time. However, a too large value for the sub-problem tolerance may lead to slow convergence or even dibergence in the outer iteration.
366    \end{methoddesc}
367    
368  \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False \optional{sub_rtol=1.e-8}}}}  \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False \optional{sub_rtol=1.e-8}}}}
369  solves the problem. and returns approximations for the flux $v$ and the pressure $p$.  solves the problem. and returns approximations for the flux $v$ and the pressure $p$.
370  \var{u0} and \var{p0} define initial guess for flux and pressure. Values marked  \var{u0} and \var{p0} define initial guess for flux and pressure. Values marked
371  by positive values \var{location_of_fixed_flux} and \var{location_of_fixed_pressure}, respectively, are kept unchanged.  by positive values \var{location_of_fixed_flux} and \var{location_of_fixed_pressure}, respectively, are kept unchanged.
 \var{sub_rtol} defines the tolerance used to solve the involved PDEs. \var{sub_rtol} needs to be choosen sufficiently small to ensure convergence but users need to keep in mind that a very small value for \var{sub_rtol} will result in a long compute time. Typically  $\var{sub_rtol}=\var{rtol}^2$ is a good choice if $\var{rtol}$ is not choosen too small.  
372  \end{methoddesc}  \end{methoddesc}
373    
374    
375  \subsection{Example: Gravity Flow}  \subsection{Example: Gravity Flow}
376  later  later
377    
378  %================================================  %\input{levelsetmodel}
 % \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}  
379    
 % where $\vec{v}$ is the velocity vector, $T$ is the temperature, $\rho$ is the density, $\eta$ is the viscosity, % % $c\hackscore{p}$ is the specific heat at constant pressure and $k$ is the thermal conductivity.  
380    
 % \subsection{Description}  
   
 % \subsection{Method}  
 %  
 % \begin{classdesc}{TemperatureCartesian}{dom,theta=THETA,useSUPG=SUPG}  
 % \end{classdesc}  
   
 % \subsection{Benchmark Problem}  
 %===============================================================================================================  
   
 %=========================================================  
 \input{levelsetmodel}  
   
 % \section{Drucker Prager Model}  
381    
382  \section{Isotropic Kelvin Material \label{IKM}}  \section{Isotropic Kelvin Material \label{IKM}}
383  As proposed by Kelvin~\ref{KELVN} material strain $D\hackscore{ij}=v\hackscore{i,j}+v\hackscore{j,i}$ can be decomposed into  As proposed by Kelvin~\cite{Muhlhaus2005} material strain $D\hackscore{ij}=\frac{1}{2}(v\hackscore{i,j}+v\hackscore{j,i})$ can be decomposed into
384  an elastic part $D\hackscore{ij}^{el}$ and visco-plastic part $D\hackscore{ij}^{vp}$:  an elastic part $D\hackscore{ij}^{el}$ and visco-plastic part $D\hackscore{ij}^{vp}$:
385  \begin{equation}\label{IKM-EQU-2}  \begin{equation}\label{IKM-EQU-2}
386  D\hackscore{ij}=D\hackscore{ij}^{el}+D\hackscore{ij}^{vp}  D\hackscore{ij}=D\hackscore{ij}^{el}+D\hackscore{ij}^{vp}
# Line 421  betwee the the strain in material $q$ Line 404  betwee the the strain in material $q$
404  \eta^{q}=\eta^{q}\hackscore{N} \left(\frac{\tau}{\tau\hackscore{t}^q}\right)^{\frac{1}{n^{q}}-1}  \eta^{q}=\eta^{q}\hackscore{N} \left(\frac{\tau}{\tau\hackscore{t}^q}\right)^{\frac{1}{n^{q}}-1}
405  \mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'\hackscore{ij} \sigma'\hackscore{ij}}  \mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'\hackscore{ij} \sigma'\hackscore{ij}}
406  \end{equation}  \end{equation}
407  for a given power law coefficients $n^{q}\ge1$ and transition stresses $\tau\hackscore{t}^q$, see~\ref{POERLAW}.  for a given power law coefficients $n^{q}\ge1$ and transition stresses $\tau\hackscore{t}^q$, see~\cite{Muhlhaus2005}.
408  Notice that $n^{q}=1$ gives a constant viscosity.  Notice that $n^{q}=1$ gives a constant viscosity.
409  After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:  After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:
410  \begin{equation}\label{IKM-EQU-6}  \begin{equation}\label{IKM-EQU-6}
411  D\hackscore{ij}'^{vp}=\frac{1}{2 \eta^{vp}} \sigma'\hackscore{ij} \mbox{ with } \frac{1}{\eta^{vp}} = \sum\hackscore{q} \frac{1}{\eta^{q}} \;.  D\hackscore{ij}'^{vp}=\frac{1}{2 \eta^{vp}} \sigma'\hackscore{ij} \mbox{ with } \frac{1}{\eta^{vp}} = \sum\hackscore{q} \frac{1}{\eta^{q}} \;.
412  \end{equation}  \end{equation}
413  With  and finally with~\ref{IKM-EQU-2}
414  \begin{equation}\label{IKM-EQU-8}  \begin{equation}\label{IKM-EQU-2bb}
415  \dot{\gamma}=\sqrt{2 D\hackscore{ij} D\hackscore{ij}}  D\hackscore{ij}'=\frac{1}{2 \eta^{vp}} \sigma'\hackscore{ij}+\frac{1}{2 \mu} \dot{\sigma}\hackscore{ij}'
 \end{equation}  
 one gets  
 \begin{equation}\label{IKM-EQU-8b}  
 \tau = \eta^{vp} \dot{\gamma}^{vp} \;.  
416  \end{equation}  \end{equation}
417  With the Drucker-Prager cohesion factor $\tau\hackscore{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$ we want to achieve  The total stress $\tau$ needs to fullfill the yield condition \index{yield condition}
418  \begin{equation}\label{IKM-EQU-8c}  \begin{equation}\label{IKM-EQU-8c}
419  \tau \le \tau\hackscore{Y} + \beta \; p  \tau \le \tau\hackscore{Y} + \beta \; p
420  \end{equation}  \end{equation}
421  which leads to the condition  with the Drucker-Prager \index{Druck-Prager} cohesion factor \index{cohesion factor} $\tau\hackscore{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$.
 \begin{equation}\label{IKM-EQU-8d}  
 \eta^{vp} \le \frac{\tau\hackscore{Y} + \beta \; p}{ \dot{\gamma}^{vp}} \; .  
 \end{equation}  
 Therefore we modify the definition of $\eta^{vp}$ to the form  
 \begin{equation}\label{IKM-EQU-6b}  
 \frac{1}{\eta^{vp}}=\max(\sum\hackscore{q} \frac{1}{\eta^{q}}, \frac{\dot{\gamma}^{vp}} {\tau\hackscore{Y} + \beta \; p})  
 \end{equation}  
422  The deviatoric stress needs to fullfill the equilibrion equation  The deviatoric stress needs to fullfill the equilibrion equation
423  \begin{equation}\label{IKM-EQU-1}  \begin{equation}\label{IKM-EQU-1}
424  -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}  -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}
425  \end{equation}  \end{equation}
426  where $F\hackscore{j}$ is a given external fource. We assume an incompressible media:  where $F\hackscore{j}$ is a given external fource. We assume an incompressible media:
427  \begin{equation}\label{IKM-EQU-2}  \begin{equation}\label{IKM-EQU-2bbb}
428  -v\hackscore{i,i}=0  -v\hackscore{i,i}=0
429  \end{equation}  \end{equation}
430  Natural boundary conditions are taken in the form  Natural boundary conditions are taken in the form
# Line 468  where the index $i$ may depend on the lo Line 440  where the index $i$ may depend on the lo
440  \subsection{Solution Method \label{IKM-SOLVE}}  \subsection{Solution Method \label{IKM-SOLVE}}
441  By using a first order finite difference approximation wit step size $dt>0$~\ref{IKM-EQU-3} get the form  By using a first order finite difference approximation wit step size $dt>0$~\ref{IKM-EQU-3} get the form
442  \begin{equation}\label{IKM-EQU-3b}  \begin{equation}\label{IKM-EQU-3b}
443  D\hackscore{ij}'^{el}=\frac{1}{2 \mu dt } \left( \sigma\hackscore{ij}' - \sigma\hackscore{ij}^{'-} \right)  \dot{\sigma}\hackscore{ij}=\frac{1}{dt } \left( \sigma\hackscore{ij} - \sigma\hackscore{ij}^{-} \right)
444    \end{equation}
445    and
446    \begin{equation}\label{IKM-EQU-2b}
447    D\hackscore{ij}'=\left(\frac{1}{2 \eta^{vp}} + \frac{1}{2 \mu dt}\right) \sigma\hackscore{ij}'-\frac{1}{2 \mu dt } \sigma\hackscore{ij}^{-'}
448  \end{equation}  \end{equation}
449  where $\sigma\hackscore{ij}^{'-}$ is the deviatoric stress at the precious time step.  where $\sigma\hackscore{ij}^{-}$ is the stress at the precious time step. With
450  Now we can combine equations~\ref{IKM-EQU-2}, \ref{IKM-EQU-3b} and~\ref{IKM-EQU-6b} to get  \begin{equation}\label{IKM-EQU-2c}
451    \dot{\gamma} = \sqrt{ 2 \left( D\hackscore{ij}' +
452    \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{-'}\right)^2}
453    \end{equation}
454    we have
455    \begin{equation}
456    \tau = \eta\hackscore{eff} \cdot \dot{\gamma}
457    \end{equation}
458    where
459    \begin{equation}
460    \eta\hackscore{eff}= min( \left(\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}\right)^{-1}
461    , \eta\hackscore{max}) \mbox{ with }
462    \eta\hackscore{max} = \left\{
463    \begin{array}{rcl}
464    \frac{\tau\hackscore{Y} + \beta \; p}{\dot{\gamma}} & & \dot{\gamma}>0 \\
465    &\mbox{ if } \\
466    \infty & & \mbox{otherwise}
467    \end{array}
468    \right.
469    \end{equation}
470    The upper bound $\eta\hackscore{max}$ makes sure that yield condtion~\ref{IKM-EQU-8c} holds. With this setting the eqaution \ref{IKM-EQU-2b} takes the form
471  \begin{equation}\label{IKM-EQU-10}  \begin{equation}\label{IKM-EQU-10}
472  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff}  \left( D\hackscore{ij}' +  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff}  \left( D\hackscore{ij}' +
473  \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)  \mbox{ with }  \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)  
 \frac{1}{\eta\hackscore{eff}}=\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}  
474  \end{equation}  \end{equation}
 Notice that $\eta\hackscore{eff}$ is a function of diatoric stress $\sigma\hackscore{ij}'$.  
475  After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get  After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get
476  \begin{equation}\label{IKM-EQU-1ib}  \begin{equation}\label{IKM-EQU-1ib}
477  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})
478  \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+  \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+
479  \frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij,j}^{'-}   \left(\frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij}^{'-} \right)\hackscore{,j}
480  \end{equation}  \end{equation}
481  Combining this with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a  Combining this with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a
482  Stokes problem as discussed in section~\ref{STOKES SOLVE} in each time step.  Stokes problem as discussed in section~\ref{STOKES SOLVE} in each time step.
483  In oder to perform step~\ref{IKM iteration 2} we need to calculate the $\eta\hackscore{eff}$ which  
484  is a function of $\sigma\hackscore{ij}$ via $\tau$.  To get $\tau$ and $\eta\hackscore{eff}$ we need to solve the  If we set
485  non-linear equation  \begin{equation}\label{IKM-EQU-44}
486  \begin{equation}  \frac{1}{\eta(\tau)}= \frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
487  \tau = \eta\hackscore{eff} \cdot \dot{\gamma}\hackscore{total} \mbox{ with }  \end{equation}
488  \dot{\gamma}\hackscore{total} = \sqrt{ 2 \left( D\hackscore{ij}' +  we need to solve the nonlinear problem
489  \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)^2}  \begin{equation}
490  \label{IKM iteration 5}  \eta\hackscore{eff} -  min(\eta( \dot{\gamma} \cdot \eta\hackscore{eff})
491  \end{equation}  , \eta\hackscore{max}) =0
492  The Newton scheme takes the form  \end{equation}
493  \begin{equation}  We use the Newton-Raphson Scheme \index{Newton-Raphson scheme} to solve this problem
494  \tau\hackscore{n+1} = \min(\tau\hackscore{n} - \frac{\tau\hackscore{n} - \eta\hackscore{eff}  \cdot \dot{\gamma}\hackscore{total}}{1 - \eta\hackscore{eff}'  \cdot  \dot{\gamma}\hackscore{total}}, \tau\hackscore{Y} + \beta \; p)  \begin{equation}\label{IKM-EQU-45}
495  \label{IKM iteration 6}  \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
496  \end{equation}  \eta\hackscore{eff}^{(n)} -
497  where $\eta\hackscore{eff}'$ denotes the derivative of $\eta\hackscore{eff}$ with respect of $\tau$. The second term in $\min$ is droped of $\tau\hackscore{Y} + \beta \; p<0$ (?)). We have  \frac{\eta\hackscore{eff}^{(n)} - \eta( \tau^{(n)}) }
498    {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
499    =\min(\eta\hackscore{max},
500    \frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot  \eta'( \tau^{(n)} )  }
501    {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
502    \end{equation}
503    where $\eta'$ denotes the derivative of $\eta$ with respect of $\tau$
504    and $\tau^{(n)} =  \dot{\gamma} \cdot \eta\hackscore{eff}^{(n)}$.
505    
506    Looking at the evaluation of $\eta$ in~\ref{IKM-EQU-44} it makes sense formulate
507    the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$.
508    In fact we have
509  \begin{equation}  \begin{equation}
510  \eta\hackscore{eff}' = - \eta\hackscore{eff}^2 \left(\frac{1}{\eta\hackscore{eff}}\right)'  \eta' = - \frac{\Theta'}{\Theta^2}
511  \mbox{ with }  \mbox{ with }
512  \left(\frac{1}{\eta\hackscore{eff}}\right)' = \sum\hackscore{q} \left(\frac{1}{\eta^{q}} \right)'  \Theta' = \sum\hackscore{q} \left(\frac{1}{\eta^{q}} \right)'
513  \label{IKM iteration 7}  \label{IKM iteration 7}
514  \end{equation}  \end{equation}
515  \begin{equation}\label{IKM-EQU-5XX}  As
516    \begin{equation}\label{IKM-EQU-47}
517  \left(\frac{1}{\eta^{q}} \right)'  \left(\frac{1}{\eta^{q}} \right)'
518  = \frac{1-\frac{1}{n^{q}}}{\eta^{q}\hackscore{N}} \frac{\tau^{-\frac{1}{n^{q}}}}{(\tau\hackscore{t}^q)^{1-\frac{1}{n^{q}}}}  = \frac{1-\frac{1}{n^{q}}}{\eta^{q}\hackscore{N}} \frac{\tau^{-\frac{1}{n^{q}}}}{(\tau\hackscore{t}^q)^{1-\frac{1}{n^{q}}}}
519  = \frac{1-\frac{1}{n^{q}}}{\tau}\frac{1}{\eta^{q}}  = \frac{1-\frac{1}{n^{q}}}{\eta^{q}}\frac{1}{\tau}
520  \end{equation}  \end{equation}
521  Notice that allways $\eta\hackscore{eff}'\le 0$ which makes the denomionator in~\ref{IKM iteration 6}  we have
522  positive.  \begin{equation}\label{IKM-EQU-48}
523    \Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum\hackscore{q}\frac{1-\frac{1}{n^{q}}}{\eta^{q}}
524    \end{equation}
525    which leads to
526    \begin{equation}\label{IKM-EQU-49}
527    \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
528    \eta\hackscore{eff}^{(n)}
529    \frac{\Theta^{(n)}  + \omega^{(n)}  }
530    {\eta\hackscore{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })
531    \end{equation}
532    
533    
534    \subsection{Functions}
535    
536    #=============================================================================
537    \begin{classdesc}{IncompressibleIsotropicFlowCartesian}{
538    domain
539    \optional{, stress=0
540    \optional{, v=0
541    \optional{, p=0
542    \optional{, t=0
543    \optional{, numMaterials=1
544    \optional{, verbose=True}}}}}}}
545    opens an incompressible, isotropic flow problem in Cartesian cooridninates.
546    \var{stress},
547    \var{v},
548    \var{p}, and
549    \var{t} set the initial deviatoric stress, velocity, pressure and time.
550    \var{numMaterials} specifies the number of materials used in the power law
551    model. Some progress information are printed if \var{verbose} is set to
552    \True.
553    \end{classdesc}
554    
555     setExternals(self, F=None, f=None, fixed_v_mask=None, v_boundary=None):
556    \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setExternals}{\optional{F=None, \optional{f=None, \optional{fixed_v_mask=None, \optional{v_boundary=None}}}}}
557    assigns values to external forces and boundary conditions. Between two calls only variables with a new values need to be set.
558    
559    In any call all values must be set.
560    \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,
561    \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.
562    The locations and compontents where the velocity is fixed are set by
563    the values of \var{fixed_u_mask}. The method will try to cast the given values to appropriate
564    \Data class objects.
565    \end{methoddesc}
566    
567    \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p,
568    \optional{max_iter=20, \optional{verbose=False, \optional{usePCG=True}}}}
569    solves the problem and return approximations for velocity and pressure.
570    The arguments \var{v} and \var{p} define initial guess. The values of \var{v} marked
571    by \var{fixed_u_mask} remain unchanged.
572    If \var{usePCG} is set to \True
573    reconditioned conjugate gradient method (PCG) \index{preconditioned conjugate gradient method!PCG}  scheme is used. Otherwise the problem is solved generalized minimal residual method (GMRES) \index{generalized minimal residual method!GMRES}. In most cases
574    the PCG scheme is more efficient.
575    \var{max_iter} defines the maximum number of iteration steps.
576    If \var{verbose} is set to \True informations on the progress of of the solver are printed.
577    \end{methoddesc}
578    
579    
580    \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}}
581    sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1.
582    \end{methoddesc}
583    \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{}
584    returns the current relative tolerance.
585    \end{methoddesc}
586    \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}}
587    sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the
588    absolute talerance is set to 0.
589    \end{methoddesc}
590    \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{}
591    sreturns the current absolute tolerance.
592    \end{methoddesc}
593    \begin{methoddesc}[StokesProblemCartesian]{setSubProblemTolerance}{\optional{rtol=None}}
594    sets the tolerance to solve the involved PDEs. The subtolerance \var{rtol} should not be choosen to large
595    in order to avoid feed back of errors in the subproblem solution into the outer iteration.
596    On the otherhand is choosen to small compute time is wasted.
597    If \var{rtol} is set to \var{None} the sub-tolerance is set automatically depending on the
598    tolerance choosen for the oter iteration.
599    \end{methoddesc}
600    \begin{methoddesc}[StokesProblemCartesian]{getSubProblemTolerance}{}
601    return the tolerance for the involved PDEs.
602    \end{methoddesc}
603    
604    
605    % \section{Drucker Prager Model}

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

  ViewVC Help
Powered by ViewVC 1.1.26