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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 9 months ago) by jfenwick
File MIME type: application/x-tex
File size: 22069 byte(s)
Don't panic.
Updating copyright stamps

1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2010 by University of Queensland
5 % Earth Systems Science Computational Center (ESSCC)
6 % http://www.uq.edu.au/esscc
7 %
8 % Primary Business: Queensland, Australia
9 % Licensed under the Open Software License version 3.0
10 % http://www.opensource.org/licenses/osl-3.0.php
11 %
12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14
15 \chapter{Models}
16 \label{MODELS CHAPTER}
17
18 The following sections give a breif overview of the model classes and their corresponding methods.
19
20 \input{stokessolver}
21
22
23 \section{Darcy Flux}
24 We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$ solving
25 the Darcy flux problem \index{Darcy flux}\index{Darcy flow}
26 \begin{equation}\label{DARCY PROBLEM}
27 \begin{array}{rcl}
28 u\hackscore{i} + \kappa\hackscore{ij} p\hackscore{,j} & = & g\hackscore{i} \\
29 u\hackscore{k,k} & = & f
30 \end{array}
31 \end{equation}
32 with the boundary conditions
33 \begin{equation}\label{DARCY BOUNDARY}
34 \begin{array}{rcl}
35 u\hackscore{i} \; n\hackscore{i} = u^{N}\hackscore{i} \; n\hackscore{i} & \mbox{ on } & \Gamma\hackscore{N} \\
36 p = p^{D} & \mbox{ on } & \Gamma\hackscore{D} \\
37 \end{array}
38 \end{equation}
39 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
40 \begin{equation}
41 \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}
42 \end{equation}
43 for all $x\hackscore{i}$.
44
45
46 \subsection{Solution Method \label{DARCY SOLVE}}
47 In practical applications it is an advantage to calculate the pressure $p$ as a correction of a 'static' pressure $p^{ref}$ which is the solution of
48 \begin{equation}
49 -(\kappa\hackscore{ki}\kappa\hackscore{kj} p\hackscore{,j}^{ref})\hackscore{,i} = - (\kappa\hackscore{ki} (g\hackscore{k}- u^{N}\hackscore{k}))\hackscore{,i}
50 \mbox{ with }
51 p^{ref} = p^{D} \mbox{ on } \Gamma\hackscore{D}
52 \end{equation}
53 With setting $u \leftarrow u-u^{N}$ and $p \leftarrow p-p^{ref}$ and
54 \begin{equation}
55 \begin{array}{rcl}
56 g\hackscore{i} & \leftarrow & g\hackscore{i} - u^{N}\hackscore{i} - \kappa\hackscore{ij} p^{ref}\hackscore{,j }\\
57 f & \leftarrow & f - u^{N}\hackscore{k,k}
58 \end{array}
59 \end{equation}
60 we can assume that $u^{N}\hackscore{i} \; n\hackscore{i}=0$ and
61 $p^{D}=0$.
62 We set
63 \begin{equation}
64 V = \{ q \in H^{1}(\Omega) : q=0 \mbox{ on } \Gamma\hackscore{D} \}
65 \end{equation}
66 and
67 \begin{equation}
68 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} \}
69 \end{equation}
70 and define the operator $Q: V \rightarrow (L^2(\Omega))^{d}$ defined by
71 \begin{equation}
72 (Qp)\hackscore{i} = \kappa\hackscore{ij} p\hackscore{,j}
73 \end{equation}
74 and the operator $D: W \rightarrow L^2(\Omega)$ defined by
75 \begin{equation}
76 Dv = v\hackscore{k,k}
77 \end{equation}
78 In operator notation the Darcy problem~\ref{DARCY PROBLEM} is written in the form
79 \begin{equation}
80 \begin{array}{rcl}
81 u + Qp & = & g \\
82 Du & = & f
83 \end{array}
84 \end{equation}
85 We solve this equation by minimising the functional
86 \begin{equation}
87 J(u,p):=\|u + Qp - g\|^2\hackscore{0} + \|Du-f\|\hackscore{0}^2
88 \end{equation}
89 over $W \times V$ where $\|.\|\hackscore{0}$ denotes the norm in $L^2(\Omega)$. A simple calculation shows that
90 one has to solve
91 \begin{equation}
92 ( v + Qq , u + Qp - g) + (Dv,Du-f) =0
93 \end{equation}
94 for all $v\in W$ and $q \in V$.which translates back into operator notation
95 \begin{equation}
96 \begin{array}{rcl}
97 (I+D^*D)u + Qp & = & D^*f + g \\
98 Q^*u + Q^*Q p & = & Q^*g \\
99 \end{array}
100 \end{equation}
101 where $D^*$ and $Q^*$ denote the adjoint operators.
102 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
103 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$)
104
105 The approach we are taking is to eliminate the velocity $u$ from the problem. Assuming that $p$ is known we have
106 \begin{equation}\label{DARCY V FORM}
107 v= (I+D^*D)^{-1} (D^*f + g - Qp)
108 \end{equation}
109 (notice that $(I+D^*D)$ is coercive in $W$) which is inserted into the second equation
110 \begin{equation}
111 Q^* (I+D^*D)^{-1} (D^*f + g - Qp) + Q^* Q p = Q^*g
112 \end{equation}
113 which is
114 \begin{equation}
115 Q^* ( I - (I+D^*D)^{-1} ) Q p = Q^* (g-(I+D^*D)^{-1} (D^*f + g) )
116 \end{equation}
117 We use the PCG \index{linear solver!PCG}\index{PCG} method to solve this. The residual $r$ ($\in V^*$) is given as
118 \begin{equation}
119 \begin{array}{rcl}
120 r & = & Q^* \left( g -(I+D^*D)^{-1} (D^*f + g) - Qp + (I+D^*D)^{-1}Q p \right)\\
121 & =& Q^* \left( g- Qp - (I+D^*D)^{-1} (D^*f + g - Qp) \right) \\
122 & =& Q^* \left( g - Qp - v \right)
123 \end{array}
124 \end{equation}
125 So in a partical implementation we use $\hat{r}=g-Qp-v$ to represent the residual.
126 The evaluation of the iteration operator for a given $p$ is then
127 returning $Qp+v$ where $v$ is the solution of
128 \begin{equation}\label{UPDATE W}
129 (I+D^*D)v = Qp
130 \end{equation}
131 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
132 implemented by solving
133 \begin{equation}\label{UPDATE P}
134 Q^*Q q = Q^*\hat{r}
135 \end{equation}
136 The residual norm used in the PCG is given as
137 \begin{equation}\label{DARCY R NORM}
138 \|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
139 \|\hat{r}\|\hackscore{0}^2
140 \end{equation}
141 The iteration is terminated if
142 \begin{equation}\label{DARCY STOP}
143 \|r\|\hackscore{PCG} \le \mbox{ATOL}
144 \end{equation}
145 where we set
146 \begin{equation}\label{DARCY ATOL DEF}
147 \mbox{ATOL} = \mbox{atol} + \mbox{rtol} \cdot \left(\frac{1}{\|v\|\hackscore{0}} + \frac{1}{\|Qp\|\hackscore{0}} \right)^{-1}
148 \end{equation}
149 where rtol is a given relative tolerance and $\mbox{atol}$ is a given absolute tolerance (typically $=0$).
150 Notice that if $Qp$ and $v$ both are zero, the pair $(0,p)$ is a solution.
151 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$
152 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.
153
154 \subsection{Functions}
155 \begin{classdesc}{DarcyFlow}{domain \optional{, adaptSubTolerance=\True}}
156 opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.
157 If \var{adaptSubTolerance} is set to \True,
158 the relative tolerances for solving~(\ref{DARCY V FORM}),~(\ref{UPDATE W})
159 and~(\ref{UPDATE P}) are set automatically.
160 \end{classdesc}
161
162 \begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None,
163 \\\optional{permeability=None}}}}}}
164 assigns values to the model parameters. Values can be assigned using various calls - in particular
165 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).
166 \var{f} and \var{g} are the corresponding parameters in~\ref{DARCY PROBLEM}.
167 The locations and compontents where the flux is prescribed are set by positive values in
168 \var{location_of_fixed_flux}.
169 The locations where the pressure is prescribed are set by
170 by positive values of \var{location_of_fixed_pressure}.
171 The values of the pressure and flux are defined by the initial guess.
172 Notice that at any point on the boundary of the domain the pressure or the normal component of
173 the flux must be defined. There must be at least one point where the pressure is prescribed.
174 The method will try to cast the given values to appropriate
175 \Data class objects.
176 \end{methoddesc}
177
178 \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e-4}}
179 sets the relative tolerance \mbox{rtol} in \ref{DARCY ATOL DEF}.
180 \end{methoddesc}
181
182 \begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}}
183 sets the absolute tolerance \mbox{atol} in \ref{DARCY ATOL DEF}.
184 \end{methoddesc}
185
186 \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{}
187 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
188 be overwritten before the solver is called.
189 \end{methoddesc}
190
191 \begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{}
192 Returns the solver options used to solve the pressure problems~(\ref{UPDATE P}).
193 Use the returned \SolverOptions object to control the solution algorithms. If the adaption of subtolerance is choosen, the tolerance will
194 be overwritten before the solver is called.
195 \end{methoddesc}
196
197 \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False}}}
198 solves the problem. and returns approximations for the flux $v$ and the pressure $p$.
199 \var{u0} and \var{p0} define initial guess for flux and pressure. Values marked
200 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.
201 \end{methoddesc}
202
203
204 \subsection{Example: Gravity Flow}
205 later
206
207 %\input{levelsetmodel}
208
209
210
211 \section{Isotropic Kelvin Material \label{IKM}}
212 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
213 an elastic part $D\hackscore{ij}^{el}$ and visco-plastic part $D\hackscore{ij}^{vp}$:
214 \begin{equation}\label{IKM-EQU-2}
215 D\hackscore{ij}=D\hackscore{ij}^{el}+D\hackscore{ij}^{vp}
216 \end{equation}
217 with the elastic strain given as
218 \begin{equation}\label{IKM-EQU-3}
219 D\hackscore{ij}^{el'}=\frac{1}{2 \mu} \dot{\sigma}\hackscore{ij}'
220 \end{equation}
221 where $\sigma'\hackscore{ij}$ is the deviatoric stress (Notice that $\sigma'\hackscore{ii}=0$).
222 If the material is composed by materials $q$ the visco-plastic strain can be decomposed as
223 \begin{equation}\label{IKM-EQU-4}
224 D\hackscore{ij}^{vp'}=\sum\hackscore{q} D\hackscore{ij}^{q'}
225 \end{equation}
226 where $D\hackscore{ij}^{q}$ is the strain in material $q$ given as
227 \begin{equation}\label{IKM-EQU-5}
228 D\hackscore{ij}^{q'}=\frac{1}{2 \eta^{q}} \sigma'\hackscore{ij}
229 \end{equation}
230 where $\eta^{q}$ is the viscosity of material $q$. We assume the following
231 betwee the the strain in material $q$
232 \begin{equation}\label{IKM-EQU-5b}
233 \eta^{q}=\eta^{q}\hackscore{N} \left(\frac{\tau}{\tau\hackscore{t}^q}\right)^{1-n^{q}}
234 \mbox{ with } \tau=\sqrt{\frac{1}{2}\sigma'\hackscore{ij} \sigma'\hackscore{ij}}
235 \end{equation}
236 for a given power law coefficients $n^{q}\ge1$ and transition stresses $\tau\hackscore{t}^q$, see~\cite{Muhlhaus2005}.
237 Notice that $n^{q}=1$ gives a constant viscosity.
238 After inserting equation~\ref{IKM-EQU-5} into equation \ref{IKM-EQU-4} one gets:
239 \begin{equation}\label{IKM-EQU-6}
240 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}} \;.
241 \end{equation}
242 and finally with~\ref{IKM-EQU-2}
243 \begin{equation}\label{IKM-EQU-2bb}
244 D\hackscore{ij}'=\frac{1}{2 \eta^{vp}} \sigma'\hackscore{ij}+\frac{1}{2 \mu} \dot{\sigma}\hackscore{ij}'
245 \end{equation}
246 The total stress $\tau$ needs to fullfill the yield condition \index{yield condition}
247 \begin{equation}\label{IKM-EQU-8c}
248 \tau \le \tau\hackscore{Y} + \beta \; p
249 \end{equation}
250 with the Drucker-Prager \index{Druck-Prager} cohesion factor \index{cohesion factor} $\tau\hackscore{Y}$, Drucker-Prager friction $\beta$ and total pressure $p$.
251 The deviatoric stress needs to fullfill the equilibrion equation
252 \begin{equation}\label{IKM-EQU-1}
253 -\sigma'\hackscore{ij,j}+p\hackscore{,i}=F\hackscore{i}
254 \end{equation}
255 where $F\hackscore{j}$ is a given external fource. We assume an incompressible media:
256 \begin{equation}\label{IKM-EQU-2bbb}
257 -v\hackscore{i,i}=0
258 \end{equation}
259 Natural boundary conditions are taken in the form
260 \begin{equation}\label{IKM-EQU-Boundary}
261 \sigma'\hackscore{ij}n\hackscore{j}-n\hackscore{i}p=f
262 \end{equation}
263 which can be overwritten by a constraint
264 \begin{equation}\label{IKM-EQU-Boundary0}
265 v\hackscore{i}(x)=0
266 \end{equation}
267 where the index $i$ may depend on the location $x$ on the bondary.
268
269 \subsection{Solution Method \label{IKM-SOLVE}}
270 By using a first order finite difference approximation wit step size $dt>0$~\ref{IKM-EQU-3} get the form
271 \begin{equation}\label{IKM-EQU-3b}
272 \dot{\sigma}\hackscore{ij}=\frac{1}{dt } \left( \sigma\hackscore{ij} - \sigma\hackscore{ij}^{-} \right)
273 \end{equation}
274 and
275 \begin{equation}\label{IKM-EQU-2b}
276 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}^{-'}
277 \end{equation}
278 where $\sigma\hackscore{ij}^{-}$ is the stress at the precious time step. With
279 \begin{equation}\label{IKM-EQU-2c}
280 \dot{\gamma} = \sqrt{ 2 \left( D\hackscore{ij}' +
281 \frac{1}{ 2 \mu \; dt} \sigma\hackscore{ij}^{-'}\right)^2}
282 \end{equation}
283 we have
284 \begin{equation}
285 \tau = \eta\hackscore{eff} \cdot \dot{\gamma}
286 \end{equation}
287 where
288 \begin{equation}
289 \eta\hackscore{eff}= min( \left(\frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}\right)^{-1}
290 , \eta\hackscore{max}) \mbox{ with }
291 \eta\hackscore{max} = \left\{
292 \begin{array}{rcl}
293 \frac{\tau\hackscore{Y} + \beta \; p}{\dot{\gamma}} & & \dot{\gamma}>0 \\
294 &\mbox{ if } \\
295 \infty & & \mbox{otherwise}
296 \end{array}
297 \right.
298 \end{equation}
299 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
300 \begin{equation}\label{IKM-EQU-10}
301 \sigma\hackscore{ij}' = 2 \eta\hackscore{eff} \left( D\hackscore{ij}' +
302 \frac{1}{ 2 \mu \; dt} \sigma\hackscore{ij}^{'-}\right)
303 \end{equation}
304 After inserting~\ref{IKM-EQU-10} into~\ref{IKM-EQU-1} we get
305 \begin{equation}\label{IKM-EQU-1ib}
306 -\left(\eta\hackscore{eff} (v\hackscore{i,j}+ v\hackscore{i,j})
307 \right)\hackscore{,j}+p\hackscore{,i}=F\hackscore{i}+
308 \left(\frac{\eta\hackscore{eff}}{\mu dt } \sigma\hackscore{ij}^{'-} \right)\hackscore{,j}
309 \end{equation}
310 Combining this with the incomressibilty condition~\ref{IKM-EQU-2} we need to solve a
311 Stokes problem as discussed in section~\ref{STOKES SOLVE} in each time step.
312
313 If we set
314 \begin{equation}\label{IKM-EQU-44}
315 \frac{1}{\eta(\tau)}= \frac{1}{\mu \; dt}+\frac{1}{\eta^{vp}}
316 \end{equation}
317 we need to solve the nonlinear problem
318 \begin{equation}
319 \eta\hackscore{eff} - min(\eta( \dot{\gamma} \cdot \eta\hackscore{eff})
320 , \eta\hackscore{max}) =0
321 \end{equation}
322 We use the Newton-Raphson Scheme \index{Newton-Raphson scheme} to solve this problem
323 \begin{equation}\label{IKM-EQU-45}
324 \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
325 \eta\hackscore{eff}^{(n)} -
326 \frac{\eta\hackscore{eff}^{(n)} - \eta( \tau^{(n)}) }
327 {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
328 =\min(\eta\hackscore{max},
329 \frac{\eta( \tau^{(n)}) -\tau^{(n)} \cdot \eta'( \tau^{(n)} ) }
330 {1-\dot{\gamma} \cdot \eta'( \tau^{(n)} )} )
331 \end{equation}
332 where $\eta'$ denotes the derivative of $\eta$ with respect of $\tau$
333 and $\tau^{(n)} = \dot{\gamma} \cdot \eta\hackscore{eff}^{(n)}$.
334
335 Looking at the evaluation of $\eta$ in~\ref{IKM-EQU-44} it makes sense formulate
336 the iteration~\ref{IKM-EQU-45} using $\Theta=\eta^{-1}$.
337 In fact we have
338 \begin{equation}
339 \eta' = - \frac{\Theta'}{\Theta^2}
340 \mbox{ with }
341 \Theta' = \sum\hackscore{q} \left(\frac{1}{\eta^{q}} \right)'
342 \label{IKM iteration 7}
343 \end{equation}
344 As
345 \begin{equation}\label{IKM-EQU-47}
346 \left(\frac{1}{\eta^{q}} \right)'
347 = \frac{n^{q}-1}{\eta^{q}\hackscore{N}} \cdot \frac{\tau^{n^{q}-2}}{(\tau\hackscore{t}^q)^{n^{q}-1}}
348 = \frac{n^{q}-1}{\eta^{q}}\cdot\frac{1}{\tau}
349 \end{equation}
350 we have
351 \begin{equation}\label{IKM-EQU-48}
352 \Theta' = \frac{1}{\tau} \omega \mbox{ with } \omega = \sum\hackscore{q}\frac{n^{q}-1}{\eta^{q}}
353 \end{equation}
354 which leads to
355 \begin{equation}\label{IKM-EQU-49}
356 \eta\hackscore{eff}^{(n+1)} = \min(\eta\hackscore{max},
357 \eta\hackscore{eff}^{(n)}
358 \frac{\Theta^{(n)} + \omega^{(n)} }
359 {\eta\hackscore{eff}^{(n)} \Theta^{(n)^2}+\omega^{(n)} })
360 \end{equation}
361
362
363 \subsection{Functions}
364
365 \begin{classdesc}{IncompressibleIsotropicFlowCartesian}{
366 domain
367 \optional{, stress=0
368 \optional{, v=0
369 \optional{, p=0
370 \optional{, t=0
371 \optional{, numMaterials=1
372 \optional{, verbose=True
373 \optional{, adaptSubTolerance=True
374 }}}}}}}}
375 opens an incompressible, isotropic flow problem in Cartesian cooridninates
376 on the domain \var{domain}.
377 \var{stress},
378 \var{v},
379 \var{p}, and
380 \var{t} set the initial deviatoric stress, velocity, pressure and time.
381 \var{numMaterials} specifies the number of materials used in the power law
382 model. Some progress information are printed if \var{verbose} is set to
383 \True. If \var{adaptSubTolerance} is equal to True the tolerances for subproblems are set automatically.
384
385 The domain
386 needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for detials~\index{LBB condition}.
387 For instance one can use second order polynomials for velocity and
388 first order polynomials for the pressure on the same element. Alternativly, one can use
389 macro elements\index{macro elements} using linear polynomial for both pressure and velocity bu with a subdivided
390 element for the velocity. Typically, the macro element is more cost effective. The fact that pressure and velocity are represented in different way is expressed by
391 \begin{python}
392 velocity=Vector(0.0, Solution(mesh))
393 pressure=Scalar(0.0, ReducedSolution(mesh))
394 \end{python}
395 \end{classdesc}
396
397 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDomain}{}
398 returns the domain.
399 \end{methoddesc}
400
401 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTime}{}
402 Returns current time.
403 \end{methoddesc}
404
405 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getStress}{}
406 Returns current stress.
407 \end{methoddesc}
408
409 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStress}{}
410 Returns current deviatoric stress.
411 \end{methoddesc}
412
413 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getPressure}{}
414 Returns current pressure.
415 \end{methoddesc}
416
417 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getVelocity}{}
418 Returns current velocity.
419 \end{methoddesc}
420
421 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getDeviatoricStrain}{}
422 Returns deviatoric strain of current velocity
423 \end{methoddesc}
424
425 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getTau}{}
426 Returns current second invariant of deviatoric stress
427 \end{methoddesc}
428
429 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{getGammaDot}{}
430 Returns current second invariant of deviatoric strain
431 \end{methoddesc}
432
433
434 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setTolerance}{tol=1.e-4}
435 Sets the tolerance used to terminate the iteration on a time step.
436 \end{methoddesc}
437
438 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setFlowTolerance}{tol=1.e-4}
439 Sets the relative tolerance for the incompressible solver, see \class{StokesProblemCartesian} for details.
440 \end{methoddesc}
441
442 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
443 Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.
444 \end{methoddesc}
445
446 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setEtaTolerance=}{rtol=1.e-8}
447 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}.
448 \end{methoddesc}
449
450 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setDruckerPragerLaw}
451 {\optional{tau_Y=None, \optional{friction=None}}}
452 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
453 condition is not applied.
454 \end{methoddesc}
455
456 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setElasticShearModulus}{mu=None}
457 Sets the elastic shear modulus $\mu$. If \var{mu} is set to None (default) elasticity is not applied.
458 \end{methoddesc}
459
460
461 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{setPowerLaws}{eta_N, tau_t, power}
462 Sets the parameters of the power-law for all materials as defined in
463 equation~\ref{IKM-EQU-5b}.
464 \var{eta_N} is the list of viscosities $\eta^{q}\hackscore{N}$,
465 \var{tau_t} is the list of reference stresses $\tau\hackscore{t}^q$,
466 and \var{power} is the list of power law coefficients $n^{q}$.
467 \end{methoddesc}
468
469
470 \begin{methoddesc}[IncompressibleIsotropicFlowCartesian]{update}{dt
471 \optional{, iter_max=100
472 \optional{, inner_iter_max=20
473 }}}
474 Updates stress, velocity and pressure for time increment \var{dt}.
475 where \var{iter_max} is the maximum number of iteration steps on a time step to
476 update the effective viscosity and \var{inner_iter_max} is the maximum
477 number of itertion steps in the incompressible solver.
478 \end{methoddesc}
479
480 \subsection{Example}
481 later
482
483 \input{faultsystem}
484
485 % \section{Drucker Prager Model}

  ViewVC Help
Powered by ViewVC 1.1.26