/[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 1966 by jfenwick, Wed Nov 5 04:59:22 2008 UTC revision 2108 by gross, Fri Nov 28 05:09:23 2008 UTC
# Line 17  Line 17 
17  The following sections give a breif overview of the model classes and their corresponding methods.  The following sections give a breif overview of the model classes and their corresponding methods.
18    
19  \section{Stokes Problem}  \section{Stokes Problem}
20  The velocity field $v$ and pressure $p$ of an incompressible fluid is given as the solution of their  The velocity \index{velocity} field $v$ and pressure $p$ of an incompressible fluid \index{incompressible fluid} is given as the solution of the Stokes problem\index{Stokes problem}
 Stokes problem  
21  \begin{equation}\label{Stokes 1}  \begin{equation}\label{Stokes 1}
22  -\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}  -\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=f\hackscore{i}-\sigma\hackscore{ij,j}
23  \end{equation}  \end{equation}
24  where $\eta$ is the viscosity and $F_i$ defines an internal force. We assume an incompressible media:  where $\eta$ is the viscosity, $F\hackscore{i}$ defines an internal force \index{force, internal} and $\sigma\hackscore{ij}$ is an intial stress \index{stress, initial}. We assume an incompressible media:
25  \begin{equation}\label{Stokes 2}  \begin{equation}\label{Stokes 2}
26  -v\hackscore{i,i}=0  -v\hackscore{i,i}=0
27  \end{equation}  \end{equation}
28  Natural boundary conditions are taken in the form  Natural boundary conditions are taken in the form
29  \begin{equation}\label{Stokes Boundary}  \begin{equation}\label{Stokes Boundary}
30  \left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)n\hackscore{j}-n\hackscore{i}p=f  \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}
31  \end{equation}  \end{equation}
32  which can be overwritten by a constraint  which can be overwritten by constraints of the form
33  \begin{equation}\label{Stokes Boundary0}  \begin{equation}\label{Stokes Boundary0}
34  v\hackscore{i}(x)=0  v\hackscore{i}(x)=v^D\hackscore{i}(x)
35  \end{equation}  \end{equation}
36  where the index $i$ may depend on the location $x$ on the bondary.  at some locations $x$ at the boundary of the domain. The index $i$ may depend on the location $x$ on the boundary.
37    $v^D$ is a given function on the domain.
38    
39  \subsection{Solution Method \label{STOKES SOLVE}}  \subsection{Solution Method \label{STOKES SOLVE}}
40  In block form equation equations~\ref{Stokes 1} and~\ref{Stokes 2} takes the form of a saddle point problem  In block form equation equations~\ref{Stokes 1} and~\ref{Stokes 2} takes the form of a saddle point problem
41    \index{saddle point problem}
42  \begin{equation}  \begin{equation}
43  \left[ \begin{array}{cc}  \left[ \begin{array}{cc}
44  A     & B^{*} \\  A     & B^{*} \\
45  B & 0 \\  B & 0 \\
46  \end{array} \right]  \end{array} \right]
47  \left[ \begin{array}{c}  \left[ \begin{array}{c}
48  u \\  v \\
49  p \\  p \\
50  \end{array} \right]  \end{array} \right]
51  =\left[ \begin{array}{c}  =\left[ \begin{array}{c}
52  F \\  G \\
53  0 \\  0 \\
54  \end{array} \right]  \end{array} \right]
55  \label{SADDLEPOINT}  \label{SADDLEPOINT}
56  \end{equation}  \end{equation}
57  where $A$ is coercive, self-adjoint linear operator in $V$, $B$ is a divergence operator and $B^{*}$ is it adjoint operator (=gradient operator)). For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}.  where $A$ is coercive, self-adjoint linear operator in a suitable Hilbert space, $B$ is the $(-1) \cdot$ divergence operator and $B^{*}$ is it adjoint operator (=gradient operator). For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}.
58  We apply the preconditioner matrix  We use iterative techniques to solve this problem. To make sure that the incomressibilty condition holds
59    with sufficient accuracy we check for
60    \begin{equation}
61    \|v\hackscore{k,k}\| \hackscore \le  \epsilon
62    \|\sqrt{v\hackscore{j,k}v\hackscore{j,k}}\|
63    \end{equation}
64    where $\epsilon$ is the desired relative accuracy and
65    \begin{equation}
66    \|p\|^2= \int\hackscore{\Omega} p^2 \; dx
67    \label{PRESSURE NORM}
68    \end{equation}
69    defines the $L^2$-norm.
70    There are two approaches to solve this problem. The first approach, called the Uzawa scheme \index{Uzawa scheme}
71    eliminates the velocity $v$ from the problem. The second approach solves the equation in coupled form after the application of a preconditioner.
72    
73    \subsubsection{Uzawa scheme}
74    The first eqution in~\ref{SADDLEPOINT} gives $v=A^{-1}(G-B^{*}p)$ assuming $p$ is known. This is inserted into the
75    second eqution which leads to
76    \begin{equation}
77    S p =  B A^{-1} G
78    \end{equation}
79    with the Schur complement \index{Schur complement} $S=BA^{-1}B^{*}$. This problem can be solved iteratively using the reconditioned Conjugate Gradient Method (PCG)~\index{PCG!Preconditioned Conjugate Gradient Method}
80    with the preconditioner $\hat{S}$ defined as $q=\hat{S}^{-1}p$ by solving
81    \begin{equation}
82    \frac{1}{\eta}q = p
83    \end{equation}
84    see~\cite{ELMAN} for more details. The evaluation of $w=Sp$ is done in the form
85    \begin{equation}
86    \begin{array}{rcl}
87    A v & = & B^{*}p \\
88    w & = & Bv \\
89    \end{array}
90    \label{EVAL PCG}
91    \end{equation}
92    The residual \index{residual}  $r=B A^{-1} G - S p$ is given as
93    \begin{equation}
94    r=B A^{-1} (G - B^* p) = Bv \mbox{ with } v = A^{-1}(G-B^{*}p)
95    \end{equation}
96    Therefore one uses the tuple $(v,Bv)$ to represent the residual of the current pressure $p$. Notice that before the iteration is started the right hand side $B A^{-1} G$ needs to be calculated. The bilinear form $(.,.)$ used is defined as
97    \begin{equation}
98    (p,(v,Bv))=\int\hackscore{\Omega} p \cdot Bv \; dx
99    \end{equation}
100    where $p$ is the pressure increment and $(v,Bv)$ represents an increment in the residual.
101    
102    \subsubsection{Coupled Solver}
103    An alternative approach to solve the saddle point problem~\ref{SADDLEPOINT} directly using an iterative such as
104    the generalized minimal residual method (GMRES) \index{generalized minimal residual method!GMRES} with a suitable
105    preconditioner. Here we use the operator
106  \begin{equation}  \begin{equation}
107  \left[ \begin{array}{cc}  \left[ \begin{array}{cc}
108  A^{-1}     & 0 \\  A^{-1}     & 0 \\
# Line 63  S^{-1} B A^{-1}  & -S^{-1} \\ Line 110  S^{-1} B A^{-1}  & -S^{-1} \\
110  \end{array} \right]  \end{array} \right]
111  \label{SADDLEPOINT PRECODITIONER}  \label{SADDLEPOINT PRECODITIONER}
112  \end{equation}  \end{equation}
113  with the Schur complement $S=BA^{-1}B^{*}$ to solve the problem iteratively. The updates $dv$ and $dp$ for  where again $S$ is the Schur complement~\cite{ELMAN}. In partice we will use an approximation $\hat{S}$ for $S$. The evaluation $(w,q)$ of the iteration operator for a given $(v,p)$ is done as
 given velocity $v$ and pressure $p$ is given as  
114  \begin{equation}  \begin{equation}
115  \begin{array}{rcl}  \begin{array}{rcl}
116  A dv & = & F-Av-B^{*}p \\  A w & = & Av+B^{*}p \\
117  S dp & = & B(v+dv) \\  \hat{S} q & = & B(w-v) \\
118  \end{array}  \end{array}
119  \label{SADDLEPOINT iteration}  \label{COUPLES SADDLEPOINT iteration}
120  \end{equation}  \end{equation}
121  This scheme is called the Uzawa scheme.  We use the inner product induced by the norm
122    \begin{equation}
123  In the case of the Stokes problem it can be shown that $S$ can be approximated by $\frac{1}{\eta}$. With this the iteration scheme can be implemented as  \|(v,p)\|^2= \int\hackscore{\Omega}  v\hackscore{i,j}  v\hackscore{i,j} + \left( \frac{p}{\eta}\right)^2\; dx
124    \label{COUPLES NORM}
125    \end{equation}
126    In PDE form~\ref{COUPLES SADDLEPOINT iteration} takes the form
127  \begin{equation}  \begin{equation}
128  \begin{array}{rcl}  \begin{array}{rcl}
129  -\left(\eta(dv\hackscore{i,j}+ dv\hackscore{i,j})\right)\hackscore{,j} & = & F\hackscore{i}+\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}-p\hackscore{,i} \\  -\left(\eta(w\hackscore{i,j}+ w\hackscore{i,j})\right)\hackscore{,j} & = & -\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i} \\
130  \frac{1}{\eta} dp & = & - \left(v\hackscore{,i}+dv\hackscore{,i} \right)\hackscore{,i} \\  \frac{1}{\eta}  q & = & - (w-v)\hackscore{i,i} \\
131  \end{array}  \end{array}
132  \label{SADDLEPOINT iteration 2}  \label{SADDLEPOINT iteration 2}
133  \end{equation}  \end{equation}
 To accelerate the convergence we are using the restarted $GMRES$ method using the norm  
 \begin{equation}  
 \|(v,p)\|^2= \int\hackscore{\Omega} \eta^2 v\hackscore{i,j}^2 + p^2 \; dx  
 \label{SADDLEPOINT iteration 3}  
 \end{equation}  
 or alternatively the $PCG$ method on the pressure only using the norm  
 \begin{equation}  
 \|p\|^2= \int\hackscore{\Omega} \frac{1}{\eta^2} p^2 \; dx  
 \label{SADDLEPOINT iteration 4}  
 \end{equation}  
134    
135    
136    \subsection{Functions}
137    
138    \begin{classdesc}{StokesProblemCartesian}{domain}
139    opens the Stokes problem\index{Stokes problem} on the \Domain domain. The approximation
140    order needs to be two.
141    \end{classdesc}
142    
143    \begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}}}}}}
144    assigns values to the model parameters. In any call all values must be set.
145    \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,
146    \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.
147    The locations and compontents where the velocity is fixed are set by
148    the values of \var{fixed_u_mask}. The method will try to cast the given values to appropriate
149    \Data class objects.
150    \end{methoddesc}
151    
152  \begin{classdesc}{StokesProblemCartesian}{domain,debug}  \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p,
153  opens the stokes equations on the \Domain domain. Setting debug=True switches the debug mode to on.  \optional{max_iter=20, \optional{verbose=False, \optional{useUzawa=True}}}}
154  \end{classdesc}  solves the problem and return approximations for velocity and pressure.
155    The arguments \var{v} and \var{p} define initial guess. The values of \var{v} marked
156    by \var{fixed_u_mask} remain unchanged.
157    If \var{useUzawa} is set to \True
158    the Uzawa\index{Uszwa} scheme is used. Otherwise the problem is solved in coupled form. In most cases
159    the Uzawa scheme is more efficient.
160    \var{max_iter} defines the maximum number of iteration steps.
161    If \var{verbose} is set to \True informations on the progress of of the solver are printed.
162    \end{methoddesc}
163    
164    
165    \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-8}}
166    sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1.
167    \end{methoddesc}
168    \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{}
169    returns the current relative tolerance.
170    \end{methoddesc}
171    \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}}
172    sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the
173    absolute talerance is set to 0.
174    \end{methoddesc}
175    \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{}
176    sreturns the current absolute tolerance.
177    \end{methoddesc}
178    \begin{methoddesc}[StokesProblemCartesian]{setSubToleranceReductionFactor}{\optional{reduction=None}}
179    sets the reduction factor for the tolerance used to solve the PDEs. A reduction factor
180    in the order of one will minimize compute time per iteration step but my slow down convergence or even lead to
181    divergency. On the other hand a very small value for the PDE tolerance could result in a wast of compute time.
182    If \var{reduction} is set to \var{None} the sub-tolerance is solved adaptively but
183    in cases a very small tolerance is set ($<10^{-6}$) it is recommended to set the
184    reduction factor by hand. This may require some experiments.
185    \end{methoddesc}
186    \begin{methoddesc}[StokesProblemCartesian]{getSubToleranceReductionFactor}{}
187    return the current reduction factor for the sub-problem tolerance.
188    \end{methoddesc}
189    
190    \subsection{Example: Lit Driven Cavity}
191     The following script \file{lit\hackscore driven\hackscore cavity.py}
192    \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
193    illustrates the usage of the \class{StokesProblemCartesian} class to solve
194    the lit driven cavity problem~\cite{LITDRIVENCAVITY}:
195    \begin{python}
196    from esys.escript import *
197    from esys.finley import Rectangle
198    from esys.escript.models import StokesProblemCartesian
199    NE=25
200    dom = Rectangle(NE,NE,order=2)
201    x = dom.getX()
202    sc=StokesProblemCartesian(dom)
203    mask= (whereZero(x[0])*[1.,0]+whereZero(x[0]-1))*[1.,0] + \
204          (whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1]
205    sc.initialize(eta=.1, fixed_u_mask= mask)
206    v=Vector(0.,Solution(dom))
207    v[0]+=whereZero(x[1]-1.)
208    p=Scalar(0.,ReducedSolution(dom))
209    v,p=sc.solve(v,p, verbose=True)
210    saveVTK("u.xml",velocity=v,pressure=p)
211    \end{python}
212    
213    \section{Darcy Flux}
214    We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$ solving
215    the Darcy flux problem \index{Darcy flux}\index{Darcy flow}
216    \begin{equation}\label{DARCY PROBLEM}
217    \begin{array}{rcl}
218    u\hackscore{i} + \kappa\hackscore{ij} p\hackscore{,j} & = & g\hackscore{i} \\
219    u\hackscore{k,k} & = & f
220    \end{array}
221    \end{equation}
222    with the boundary conditions
223    \begin{equation}\label{DARCY BOUNDARY}
224    \begin{array}{rcl}
225    u\hackscore{i} \; n\hackscore{i}  = u^{N}\hackscore{i}  \; n\hackscore{i} & \mbox{ on } & \Gamma\hackscore{N} \\
226    p = p^{D} &  \mbox{ on } & \Gamma\hackscore{D} \\
227    \end{array}
228    \end{equation}
229    where $\Gamma\hackscore{N}$ and $\Gamma\hackscore{D}$ are a partition of the boundary of $\Omega$ with $\Gamma\hackscore{D}$ non empty, $n\hackscore{i}$ is the outer normal field of the boundary of $\Omega$, $u^{N}\hackscore{i}$ and $p^{D}$ are given functions on $\Omega$, $g\hackscore{i}$ and $f$ are given source terms and $\kappa\hackscore{ij}$ is the given permability. We assume that $\kappa\hackscore{ij}$ is symmetric (which is not really required) and positive definite, i.e there are positive constants $\alpha\hackscore{0}$ and $\alpha\hackscore{1}$ wich are independent from the location in $\Omega$ such that
230    \begin{equation}
231    \alpha\hackscore{0} \; x\hackscore{i} x\hackscore{i} \le \kappa\hackscore{ij} x\hackscore{i} x\hackscore{j} \le \alpha\hackscore{1} \; x\hackscore{i} x\hackscore{i}
232    \end{equation}
233    for all $x\hackscore{i}$.
234    
 example usage:  
235    
236  solution=StokesProblemCartesian(mesh) \\  \subsection{Solution Method \label{DARCY SOLVE}}
237  solution.setTolerance(TOL) \\  Without loss of generality we can assume that $u^{N}\hackscore{i}  \; n\hackscore{i}=0$ and
238  solution.initialize(fixed\_u\_mask=b\_c,eta=eta,f=Y) \\  $p^{D}$. Otherewise one solves for $u-u^{N}$ and $p-p^{D}$ and sets
239  velocity,pressure=solution.solve(velocity,pressure,max\_iter=max\_iter,solver=solver) \\  \begin{equation}
240    \begin{array}{rcl}
241    g\hackscore{i} & \leftarrow & g\hackscore{i} - u^{N}\hackscore{i} -  \kappa\hackscore{ij} p^{D}\hackscore{,j }\\
242    f & \leftarrow & f - u^{N}\hackscore{k,k}
243    \end{array}
244    \end{equation}
245    We set
246    \begin{equation}
247    V = \{ q \in H^{1}(\Omega) : q=0 \mbox{ on } \Gamma\hackscore{D} \}
248    \end{equation}
249    and
250    \begin{equation}
251    W = \{ v \in (L^2(\Omega))^{d} : v\hackscore{k,k} \in L^2(\Omega) \mbox{ and } u\hackscore{i} \; n\hackscore{i} =0  \mbox{ on } \Gamma\hackscore{N} \}
252    \end{equation}
253    and define the operator $Q: V \rightarrow (L^2(\Omega))^{d}$ defined by
254    \begin{equation}
255    (Qp)\hackscore{i} = \kappa\hackscore{ij} p\hackscore{,j}
256    \end{equation}
257    and the operator $D: W \rightarrow L^2(\Omega)$ defined by
258    \begin{equation}
259    Dv = v\hackscore{k,k}
260    \end{equation}
261    In operator notation the Darcy problem~\ref{DARCY PROBLEM} is written in the form
262    \begin{equation}
263    \begin{array}{rcl}
264    u + Qp & = & g \\
265    Du & = & f
266    \end{array}
267    \end{equation}
268    We solve this equation by minimising the functional
269    \begin{equation}
270    J(u,p):=\|u + Qp - g\|^2\hackscore{0} + \|Du-f\|\hackscore{0}^2
271    \end{equation}
272    over $W \times V$ where $\|.\|\hackscore{0}$ denotes the norm in $L^2(\Omega)$. A simple calculation shows that
273    one has to solve
274    \begin{equation}
275    ( v + Qq , u + Qp - g) + (Dv,Du-f) =0
276    \end{equation}
277    for all $v\in W$ and $q \in V$.which translates back into operator notation
278    \begin{equation}
279    \begin{array}{rcl}
280    (I+D^*D)u + Qp & = & D^*f + g \\
281    Q^*u  + Q^*Q p & = & Q^* g \\
282    \end{array}
283    \end{equation}
284    where $D^*$ and $Q^*$ denote the adjoint operators.
285    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
286    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$)  
287    
288  % \subsection{Benchmark Problem}  The approach we are taking is to eliminate the velocity $u$ from the problem. Assuming that $p$ is known we have
289  %  \begin{equation}
290  % Convection problem  v= (I+D^*D)^{-1} (D^*f + g - Qp)
291    \end{equation}
292    (notice that $(I+D^*D)$ is coercive in $W$) which is inserted into the second equation
293    \begin{equation}
294    Q^* (I+D^*D)^{-1} (D^*f + g - Qp) + Q^* Q p = Q^* g
295    \end{equation}
296    which is
297    \begin{equation}
298    Q^* ( I - (I+D^*D)^{-1} ) Q p = Q^* ( g -(I+D^*D)^{-1} (D^*f + g) )
299    \end{equation}
300    We use the PCG method to solve this. The residual $r$ ($\in V^*$) is given as
301    \begin{equation}
302    \begin{array}{rcl}
303    r & = & Q^* ( g -(I+D^*D)^{-1} (D^*f + g) - Qp + (I+D^*D)^{-1}Q p \\
304    & =&  Q^* \left( (g-Qp) - (I+D^*D)^{-1} (D^*f + g - Qp) \right) \\
305    & =&  Q^* \left( (g-Qp) - v \right)
306    \end{array}
307    \end{equation}
308    So in a partical implementation we use the pair $(g-Qp,v)$ to represent the residual. This will save the
309    reconstruction of the velocity $v$. In this notation the right hand side is given as
310    $(g,(I+D^*D)^{-1} (D^*f + g))$. The evaluation of the iteration operator for a given $p$ is then
311    returning $(Qp,w)$ where $w$ is the solution of
312    \begin{equation}\label{UPDATE W}
313    (I+D^*D)w = Qp
314    \end{equation}
315    We use $Q^*Q$ as a a preconditioner for the iteration operator $Q^* ( I - (I+D^*D)^{-1} ) Q$.
316    
317    \subsection{Functions}
318    \begin{classdesc}{DarcyFlow}{domain}
319    opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.
320    \end{classdesc}
321    
322    \begin{methoddesc}[DarcyFlow]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}}}}}}
323    assigns values to the model parameters. In any call all values must be set.
324    \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,
325    \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.
326    The locations and compontents where the velocity is fixed are set by
327    the values of \var{fixed_u_mask}. The method will try to cast the given values to appropriate
328    \Data class objects.
329    \end{methoddesc}
330    
331  \section{Temperature Cartesian}  
332    \subsection{Example: Gravity Flow}
333    
334    %================================================
335    \section{Temperature Advection Diffusion\label{TEMP ADV DIFF}}
336    
337  \begin{equation}  \begin{equation}
338  \rho c\hackscore{p} \left (\frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T \right ) = k \nabla^{2}T  \rho c\hackscore{p} \left (\frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T \right ) = k \nabla^{2}T
# Line 130  where $\vec{v}$ is the velocity vector, Line 349  where $\vec{v}$ is the velocity vector,
349  \end{classdesc}  \end{classdesc}
350    
351  \subsection{Benchmark Problem}  \subsection{Benchmark Problem}
352    %===============================================================================================================
353    
354    %=========================================================
355    % \section{Level Set Method}
356    
357  \section{Level Set Method}  %\subsection{Description}
358    
359  \subsection{Description}  %\subsection{Method}
   
 \subsection{Method}  
360    
361  Advection and Reinitialisation  %Advection and Reinitialisation
362    
363  \begin{classdesc}{LevelSet}{mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth}  %\begin{classdesc}{LevelSet}{mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth}
364  \end{classdesc}  %\end{classdesc}
365    
366  %example usage:  %example usage:
367    
368  %levelset = LevelSet(mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth)  %levelset = LevelSet(mesh, func\_new, reinit\_max, reinit\_each, tolerance, smooth)
369    
370  \begin{methoddesc}[LevelSet]{update\_parameter}{parameter}  %\begin{methoddesc}[LevelSet]{update\_parameter}{parameter}
371  Update the parameter.  %Update the parameter.
372  \end{methoddesc}  %\end{methoddesc}
373    
374  \begin{methoddesc}[LevelSet]{update\_phi}{paramter}{velocity}{dt}{t\_step}  %\begin{methoddesc}[LevelSet]{update\_phi}{paramter}{velocity}{dt}{t\_step}
375  Update level set function; advection and reinitialization  %Update level set function; advection and reinitialization
376  \end{methoddesc}  %\end{methoddesc}
377    
378  \subsection{Benchmark Problem}  %\subsection{Benchmark Problem}
379    
380  Rayleigh-Taylor instability problem  %Rayleigh-Taylor instability problem
381    
382    
383  % \section{Drucker Prager Model}  % \section{Drucker Prager Model}
# Line 191  for a given power law coefficients $n^{q Line 411  for a given power law coefficients $n^{q
411  Notice that $n^{q}=1$ gives a constant viscosity.  Notice that $n^{q}=1$ gives a constant viscosity.
412  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:
413  \begin{equation}\label{IKM-EQU-6}  \begin{equation}\label{IKM-EQU-6}
414  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}} \;.
415  \end{equation}  \end{equation}
416  With  With
417  \begin{equation}\label{IKM-EQU-8}  \begin{equation}\label{IKM-EQU-8}
# Line 239  D\hackscore{ij}'^{el}=\frac{1}{2 \mu dt Line 459  D\hackscore{ij}'^{el}=\frac{1}{2 \mu dt
459  where $\sigma\hackscore{ij}^{'-}$ is the deviatoric stress at the precious time step.  where $\sigma\hackscore{ij}^{'-}$ is the deviatoric stress at the precious time step.
460  Now we can combine equations~\ref{IKM-EQU-2}, \ref{IKM-EQU-3b} and~\ref{IKM-EQU-6b} to get  Now we can combine equations~\ref{IKM-EQU-2}, \ref{IKM-EQU-3b} and~\ref{IKM-EQU-6b} to get
461  \begin{equation}\label{IKM-EQU-10}  \begin{equation}\label{IKM-EQU-10}
462  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff} D\hackscore{ij}' +  \sigma\hackscore{ij}' =  2 \eta\hackscore{eff}  \left( D\hackscore{ij}' +
463  \frac{\eta\hackscore{eff}}{\mu \; dt}  \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)  \mbox{ with }
 \sigma\hackscore{ij}^{'-} \mbox{ with }  
464  \frac{1}{\eta\hackscore{eff}}=\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}  \frac{1}{\eta\hackscore{eff}}=\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
465  \end{equation}  \end{equation}
466  Notice that $\eta\hackscore{eff}$ is a function of stress.  Notice that $\eta\hackscore{eff}$ is a function of diatoric stress $\sigma\hackscore{ij}'$.
467  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
468  \begin{equation}\label{IKM-EQU-1ib}  \begin{equation}\label{IKM-EQU-1ib}
469  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+  -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})
470    \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+
471  \frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij,j}^{'-}  \frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij,j}^{'-}
472  \end{equation}  \end{equation}
473  Together with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a problem with a form almost identical  Together with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a problem with a form almost identical
474  to the Stokes problem discussed in section~\ref{STOKES SOLVE} but with the difference that $\eta\hackscore{eff}$ is depending on the solution. Analog to the iteration scheme~\ref{SADDLEPOINT iteration 2} we can run  to the Stokes problem discussed in section~\ref{STOKES SOLVE} but with the difference that $\eta\hackscore{eff}$ is depending on the solution. Analog to the iteration scheme~\ref{SADDLEPOINT iteration 2} we can run
475  \begin{equation}  \begin{equation}
476  \begin{array}{rcl}  \begin{array}{rcl}
477  -\left(\eta\hackscore{eff}(dv\hackscore{i,j}+ dv\hackscore{i,j})\right)\hackscore{,j} & = & F\hackscore{i}+ \sigma\hackscore{ij,j}' \\  -\left(\eta\hackscore{eff}(dv\hackscore{i,j}+ dv\hackscore{i,j}
478  \frac{1}{\eta\hackscore{eff}} dp & = & - v\hackscore{i,i}^+ \\  )\right)\hackscore{,j} & = & F\hackscore{i}+ \sigma\hackscore{ij,j}'-p\hackscore{,i} \\
479  d\sigma\hackscore{ij}' & = & 2 \eta\hackscore{eff} D\hackscore{ij}^{+'} + \frac{\eta\hackscore{eff}}{\mu \; dt} \sigma\hackscore{ij}' -  \sigma\hackscore{ij}\\  \frac{1}{\eta\hackscore{eff}} dp & = & - v\hackscore{i,i}^+
480  \end{array}  \end{array}
481  \label{IKM iteration 2}  \label{IKM iteration 2}
482  \end{equation}  \end{equation}
483  where $v^+=v+dv$. As this problem is non-linear the Jacobi-free Newton-GMRES method is used with the norm  where $v^+=v+dv$. As this problem is non-linear the Jacobi-free Newton-GMRES method is used with the norm
484  \begin{equation}  \begin{equation}
485  \|(\sigma, p)\|^2= \int\hackscore{\Omega} \sigma\hackscore{i,j}^2 + p^2 \; dx  \|(v, p)\|^2= \int\hackscore{\Omega} v\hackscore{i,j}^2 + \frac{1}{\bar{\eta}^2\hackscore{eff}} p^2 \; dx
486  \label{IKM iteration 3}  \label{IKM iteration 3}
487  \end{equation}  \end{equation}
488    where  $\bar{\eta}\hackscore{eff}$ is the caracteristic viscosity, for instance:
489    \begin{equation}
490    \frac{1}{\bar{\eta}\hackscore{eff}} = \frac{1}{\tau^{-}}+\sum\hackscore{q}  \frac{1}{\eta^{q}\hackscore{N}}
491    \label{IKM iteration 4}
492    \end{equation}
493    In oder to perform step~\ref{IKM iteration 2} we need to calculate the $\eta\hackscore{eff}$ as well as $\sigma\hackscore{ij}'$ while via $\tau$ the first is a function of the latter. The priority is the
494    calculation of $\eta\hackscore{eff}$ with the Newton-Raphson scheme. This value can then be used to calculate
495    $\sigma\hackscore{ij}'$ via~\ref{IKM-EQU-10}. We need to solve
496    \begin{equation}
497    \tau = \eta\hackscore{eff} \cdot \epsilon \mbox{ with }
498    \epsilon = \sqrt{ 2 \left( D\hackscore{ij}' +
499    \frac{1}{  2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)^2}
500    \label{IKM iteration 5}
501    \end{equation}
502    The Newton scheme takes the form
503    \begin{equation}
504    \tau\hackscore{n+1} = \min(\tau\hackscore{n} - \frac{\tau\hackscore{n} - \eta\hackscore{eff}  \cdot \epsilon}{1 - \eta\hackscore{eff}'  \cdot  \epsilon}, \tau\hackscore{Y} + \beta \; p)
505    = \min(\frac{\eta\hackscore{eff} - \tau\hackscore{n}  \eta\hackscore{eff}'}
506    {1 - \eta\hackscore{eff}'  \cdot  \epsilon}, \frac{\tau\hackscore{Y} + \beta \; p}{\epsilon}) \epsilon
507    \label{IKM iteration 6}
508    \end{equation}
509    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$ or $\epsilon=0$. In fact we have
510    \begin{equation}
511    \eta\hackscore{eff}' = - \eta\hackscore{eff}^2 \left(\frac{1}{\eta\hackscore{eff}}\right)'
512    \mbox{ with }
513    \left(\frac{1}{\eta\hackscore{eff}}\right)' = \sum\hackscore{q} \left(\frac{1}{\eta^{q}} \right)'
514    \label{IKM iteration 7}
515    \end{equation}
516    \begin{equation}\label{IKM-EQU-5XX}
517    \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}}}}
519    = \frac{1-\frac{1}{n^{q}}}{ \tau \eta^{q}}
520    \end{equation}
521    Notice that allways $\eta\hackscore{eff}'\le 0$ which makes the denomionator in~\ref{IKM iteration 6}
522    positive.
523    
524    
525    

Legend:
Removed from v.1966  
changed lines
  Added in v.2108

  ViewVC Help
Powered by ViewVC 1.1.26