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

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.

 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2010 by University of Queensland 5 % Earth Systems Science Computational Center (ESSCC) 6 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 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}