/[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 2415 by gross, Wed May 13 02:48:39 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    
# Line 397  later Line 401  later
401  % \section{Drucker Prager Model}  % \section{Drucker Prager Model}
402    
403  \section{Isotropic Kelvin Material \label{IKM}}  \section{Isotropic Kelvin Material \label{IKM}}
404  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
405  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}$:
406  \begin{equation}\label{IKM-EQU-2}  \begin{equation}\label{IKM-EQU-2}
407  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 425  betwee the the strain in material $q$
425  \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}
426  \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}}
427  \end{equation}  \end{equation}
428  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}.
429  Notice that $n^{q}=1$ gives a constant viscosity.  Notice that $n^{q}=1$ gives a constant viscosity.
430  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:
431  \begin{equation}\label{IKM-EQU-6}  \begin{equation}\label{IKM-EQU-6}
432  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}} \;.
433  \end{equation}  \end{equation}
434  With  and finally with~\ref{IKM-EQU-2}
435  \begin{equation}\label{IKM-EQU-8}  \begin{equation}\label{IKM-EQU-2bb}
436  \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} \;.  
437  \end{equation}  \end{equation}
438  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}
439  \begin{equation}\label{IKM-EQU-8c}  \begin{equation}\label{IKM-EQU-8c}
440  \tau \le \tau\hackscore{Y} + \beta \; p  \tau \le \tau\hackscore{Y} + \beta \; p
441  \end{equation}  \end{equation}
442  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}  
443  The deviatoric stress needs to fullfill the equilibrion equation  The deviatoric stress needs to fullfill the equilibrion equation
444  \begin{equation}\label{IKM-EQU-1}  \begin{equation}\label{IKM-EQU-1}
445  -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}  -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}
446  \end{equation}  \end{equation}
447  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:
448  \begin{equation}\label{IKM-EQU-2}  \begin{equation}\label{IKM-EQU-2bbb}
449  -v\hackscore{i,i}=0  -v\hackscore{i,i}=0
450  \end{equation}  \end{equation}
451  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 461  where the index $i$ may depend on the lo
461  \subsection{Solution Method \label{IKM-SOLVE}}  \subsection{Solution Method \label{IKM-SOLVE}}
462  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
463  \begin{equation}\label{IKM-EQU-3b}  \begin{equation}\label{IKM-EQU-3b}
464  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)
465    \end{equation}
466    and
467    \begin{equation}\label{IKM-EQU-2b}
468    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}^{-'}
469    \end{equation}
470    where $\sigma\hackscore{ij}^{-}$ is the stress at the precious time step. With
471    \begin{equation}\label{IKM-EQU-2c}
472    \dot{\gamma} = \sqrt{ 2 \left( D\hackscore{ij}' +
473    \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{-'}\right)^2}
474    \end{equation}
475    we have
476    \begin{equation}
477    \tau = \eta\hackscore{eff} \cdot \dot{\gamma}
478    \end{equation}
479    where
480    \begin{equation}
481    \eta\hackscore{eff}= min( \left(\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}\right)^{-1}
482    , \eta\hackscore{max}) \mbox{ with }
483    \eta\hackscore{max} = \left\{
484    \begin{array}{rcl}
485    \frac{\tau\hackscore{Y} + \beta \; p}{\dot{\gamma}} & & \dot{\gamma}>0 \\
486    &\mbox{ if } \\
487    \infty & & \mbox{otherwise}
488    \end{array}
489    \right.
490  \end{equation}  \end{equation}
491  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  
492  \begin{equation}\label{IKM-EQU-10}  \begin{equation}\label{IKM-EQU-10}
493  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff}  \left( D\hackscore{ij}' +  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff}  \left( D\hackscore{ij}' +
494  \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}}  
495  \end{equation}  \end{equation}
 Notice that $\eta\hackscore{eff}$ is a function of diatoric stress $\sigma\hackscore{ij}'$.  
496  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
497  \begin{equation}\label{IKM-EQU-1ib}  \begin{equation}\label{IKM-EQU-1ib}
498  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})
499  \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+  \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+
500  \frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij,j}^{'-}   \left(\frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij}^{'-} \right)\hackscore{,j}
501  \end{equation}  \end{equation}
502  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
503  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.
504  In oder to perform step~\ref{IKM iteration 2} we need to calculate the $\eta\hackscore{eff}$ which  
505  is a function of $\sigma\hackscore{ij}$ via $\tau$.  To get $\tau$ and $\eta\hackscore{eff}$ we need to solve the  If we set
506  non-linear equation  \begin{equation}\label{IKM-EQU-44}
507  \begin{equation}  \frac{1}{\eta(\tau)}= \frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
508  \tau = \eta\hackscore{eff} \cdot \dot{\gamma}\hackscore{total} \mbox{ with }  \end{equation}
509  \dot{\gamma}\hackscore{total} = \sqrt{ 2 \left( D\hackscore{ij}' +  we need to solve the nonlinear problem
510  \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)^2}  \begin{equation}
511  \label{IKM iteration 5}  \eta\hackscore{eff} -  min(\eta( \dot{\gamma} \cdot \eta\hackscore{eff})
512  \end{equation}  , \eta\hackscore{max}) =0
513  The Newton scheme takes the form  \end{equation}
514  \begin{equation}  We use the Newton-Raphson Scheme \index{Newton-Raphson scheme} to solve this problem
515  \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}
516  \label{IKM iteration 6}  \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
517  \end{equation}  \eta\hackscore{eff}^{(n)} -
518  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)}) }
519    {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
520    =\min(\eta\hackscore{max},
521    \frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot  \eta'( \tau^{(n)} )  }
522    {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
523    \end{equation}
524    where $\eta'$ denotes the derivative of $\eta$ with respect of $\tau$
525    and $\tau^{(n)} =  \dot{\gamma} \cdot \eta\hackscore{eff}^{(n)}$.
526    
527    Looking at the evaluation of $\eta$ in~\ref{IKM-EQU-44} it makes sense formulate
528    the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$.
529    In fact we have
530  \begin{equation}  \begin{equation}
531  \eta\hackscore{eff}' = - \eta\hackscore{eff}^2 \left(\frac{1}{\eta\hackscore{eff}}\right)'  \eta' = - \frac{\Theta'}{\Theta^2}
532  \mbox{ with }  \mbox{ with }
533  \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)'
534  \label{IKM iteration 7}  \label{IKM iteration 7}
535  \end{equation}  \end{equation}
536  \begin{equation}\label{IKM-EQU-5XX}  As
537    \begin{equation}\label{IKM-EQU-47}
538  \left(\frac{1}{\eta^{q}} \right)'  \left(\frac{1}{\eta^{q}} \right)'
539  = \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}}}}
540  = \frac{1-\frac{1}{n^{q}}}{\tau}\frac{1}{\eta^{q}}  = \frac{1-\frac{1}{n^{q}}}{\eta^{q}}\frac{1}{\tau}
541  \end{equation}  \end{equation}
542  Notice that allways $\eta\hackscore{eff}'\le 0$ which makes the denomionator in~\ref{IKM iteration 6}  we have
543  positive.  \begin{equation}\label{IKM-EQU-48}
544    \Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum\hackscore{q}\frac{1-\frac{1}{n^{q}}}{\eta^{q}}
545    \end{equation}
546    which leads to
547    \begin{equation}\label{IKM-EQU-49}
548    \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
549    \eta\hackscore{eff}^{(n)}
550    \frac{\Theta^{(n)}  + \omega^{(n)}  }
551    {\eta\hackscore{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })
552    \end{equation}

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

  ViewVC Help
Powered by ViewVC 1.1.26