 # Diff of /trunk/doc/user/Models.tex

revision 2207 by gross, Mon Dec 15 05:09:02 2008 UTC revision 2208 by gross, Mon Jan 12 06:37:07 2009 UTC
# Line 249  f & \leftarrow & f - u^{N}\hackscore{k,k Line 249  f & \leftarrow & f - u^{N}\hackscore{k,k
249  \end{array}  \end{array}
250  \end{equation}  \end{equation}
251  we can assume that $u^{N}\hackscore{i} \; n\hackscore{i}=0$ and  we can assume that $u^{N}\hackscore{i} \; n\hackscore{i}=0$ and
252  $p^{D}=0$. Notice that with this setting  $p^{D}=0$.
\begin{equation}\label{DIV FREE DARCY FLUX}
(\kappa\hackscore{ki} g\hackscore{k})\hackscore{,i} = 0
\end{equation}
253  We set  We set
254  \begin{equation}  \begin{equation}
255  V = \{ q \in H^{1}(\Omega) : q=0 \mbox{ on } \Gamma\hackscore{D} \}  V = \{ q \in H^{1}(\Omega) : q=0 \mbox{ on } \Gamma\hackscore{D} \}
# Line 276  u + Qp & = & g \\ Line 273  u + Qp & = & g \\
273  Du & = & f  Du & = & f
274  \end{array}  \end{array}
275  \end{equation}  \end{equation}
Notice that because of~\ref{DIV FREE DARCY FLUX} we have
\begin{equation} \label{ABSTRACT DIV FREE DARCY FLUX}
Q^*g =0
\end{equation}
where $Q^*$ denote the adjoint operators of $Q$.
276  We solve this equation by minimising the functional  We solve this equation by minimising the functional
277  \begin{equation}  \begin{equation}
278  J(u,p):=\|u + Qp - g\|^2\hackscore{0} + \|Du-f\|\hackscore{0}^2  J(u,p):=\|u + Qp - g\|^2\hackscore{0} + \|Du-f\|\hackscore{0}^2
# Line 294  for all $v\in W$ and $q \in V$.which tra Line 286  for all $v\in W$ and $q \in V$.which tra
286  \begin{equation}  \begin{equation}
287  \begin{array}{rcl}  \begin{array}{rcl}
288  (I+D^*D)u + Qp & = & D^*f + g \\  (I+D^*D)u + Qp & = & D^*f + g \\
289  Q^*u  + Q^*Q p & = & 0 \\  Q^*u  + Q^*Q p & = & Q^*g \\
290  \end{array}  \end{array}
291  \end{equation}  \end{equation}
292  where $D^*$ and $Q^*$ denote the adjoint operators. We use the fact that $Q^*g=$.  where $D^*$ and $Q^*$ denote the adjoint operators.
293  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{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
294  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$)
295
# Line 307  v= (I+D^*D)^{-1} (D^*f + g - Qp) Line 299  v= (I+D^*D)^{-1} (D^*f + g - Qp)
299  \end{equation}  \end{equation}
300  (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
301  \begin{equation}  \begin{equation}
302  Q^* (I+D^*D)^{-1} (D^*f + g - Qp) + Q^* Q p = 0  Q^* (I+D^*D)^{-1} (D^*f + g - Qp) + Q^* Q p = Q^*g
303  \end{equation}  \end{equation}
304  which is  which is
305  \begin{equation}  \begin{equation}
306  Q^* ( I - (I+D^*D)^{-1} ) Q p = - Q^* (I+D^*D)^{-1} (D^*f + g) )  Q^* ( I - (I+D^*D)^{-1} ) Q p = Q^* (g-(I+D^*D)^{-1} (D^*f + g) )
307  \end{equation}  \end{equation}
308  We use the PCG method to solve this. The residual $r$ ($\in V^*$) is given as  We use the PCG \index{linear solver!PCG}\index{PCG} method to solve this. The residual $r$ ($\in V^*$) is given as
309  \begin{equation}  \begin{equation}
310  \begin{array}{rcl}  \begin{array}{rcl}
311  r & = & Q^*  \left( -(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)\\
312  & =&  Q^* \left( - Qp - (I+D^*D)^{-1} (D^*f + g - Qp) \right) \\  & =&  Q^* \left( - Qp - (I+D^*D)^{-1} (D^*f + g - Qp) \right) \\
313  & =&  - Q^* \left( Qp + v \right)  & =&  Q^* \left( g - Qp - v \right)
314  \end{array}  \end{array}
315  \end{equation}  \end{equation}
316  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 the pair $(Qp,v)$ to represent the residual. This will save the
# Line 330  returning $(Qp,w)$ where $w$ is the solu Line 322  returning $(Qp,w)$ where $w$ is the solu
322  \end{equation}  \end{equation}
323  We use $Q^*Q$ as a a preconditioner for the iteration operator $Q^* ( I - (I+D^*D)^{-1} ) Q$.  We use $Q^*Q$ as a a preconditioner for the iteration operator $Q^* ( I - (I+D^*D)^{-1} ) Q$.
324
325    The iteration PCG \index{linear solver!PCG}\index{PCG} is terminated if
326    \begin{equation}\label{DARCY STOP}
327    \int r \cdot (Q^*Q)^{-1} r \; dx \le \mbox{ATOL}^2
328    \end{equation}
329    where ATOL is a given absolute tolerance.
330    The initial residual $r\hackscore{0}$ is
331    \begin{equation}\label{DARCY STOP 2}
332    r\hackscore{0}=Q^* \left( g - v\hackscore{ref} \right) \mbox{ with } v\hackscore{ref} = (I+D^*D)^{-1} (D^*f + g)
333    \end{equation}
334    so the
335    \begin{equation}\label{DARCY NORM 0}
336    \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)
337    \end{equation}
338    So we set
339    \begin{equation}\label{DARCY NORM 1}
340    ATOL = atol + rtol \cdot \max(|g - v\hackscore{ref}|\hackscore{0}, |Q p\hackscore{ref} |\hackscore{0} )
341    \end{equation}
342    where atol and rtol a given absolute and relative tolerances, respectively. The reference flux $v\hackscore{ref}$
343    and reference pressure $p\hackscore{ref}$ may be calcualated from their definition which would require to solve to
344    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).
345
346  \subsection{Functions}  \subsection{Functions}
347  \begin{classdesc}{DarcyFlow}{domain}  \begin{classdesc}{DarcyFlow}{domain}
348  opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.  opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.
349  \end{classdesc}  \end{classdesc}
350    \begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None, \optional{permeability=None}}}}}}
351  \begin{methoddesc}[DarcyFlow]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}}}}}}  assigns values to the model parameters. Values can be assigned using various calls - in particular
352  assigns values to the model parameters. In any call all values must be set.  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).
353  \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,  \var{f} and \var{g} are the corresponding parameters in~\ref{DARCY PROBLEM}.
354  \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.  The locations and compontents where the flux is prescribed are set by positive values in
355  The locations and compontents where the velocity is fixed are set by  \var{location_of_fixed_flux}.
356  the values of \var{fixed_u_mask}. The method will try to cast the given values to appropriate  The locations where the pressure is prescribed are set by
357    by positive values of \var{location_of_fixed_pressure}.
358    The values of the pressure and flux are defined by the initial guess.
359    Notice that at any point on the boundary of the domain the pressure or the normal component of
360    the flux must be defined. There must be at least one point where the pressure is prescribed.
361    The method will try to cast the given values to appropriate
362  \Data class objects.  \Data class objects.
363  \end{methoddesc}  \end{methoddesc}
364
365    \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{atol=0,\optional{rtol=1e-8,\optional{p_ref=None,\optional{v_ref=None}}}}}
366    sets the absolute tolerance ATOL according to~\ref{DARCY NORM 1}. If \var{p_ref} is not present $0$ is used.
367    If \var{v_ref} is not present $0$ is used. If the final result ATOL is not positive an exception is thrown.
368    \end{methoddesc}
369
370
371
372    \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False \optional{sub_rtol=1.e-8}}}}
373    solves the problem. and returns approximations for the flux $v$ and the pressure $p$.
374    \var{u0} and \var{p0} define initial guess for flux and pressure. Values marked
375    by positive values \var{location_of_fixed_flux} and \var{location_of_fixed_pressure}, respectively, are kept unchanged.
376    \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.
377    \end{methoddesc}
378
379
380  \subsection{Example: Gravity Flow}  \subsection{Example: Gravity Flow}
381    later
382
383  %================================================  %================================================
385
386  \begin{equation}  %\begin{equation}
387  \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
388  \label{HEAT EQUATION}  % \label{HEAT EQUATION}
389  \end{equation}  % \end{equation}

where $\vec{v}$ is the velocity vector, $T$ is the temperature, $\rho$ is the density, $\eta$ is the viscosity, $c\hackscore{p}$ is the specific heat at constant pressure and $k$ is the thermal conductivity.
390
391  \subsection{Description}  % where $\vec{v}$ is the velocity vector, $T$ is the temperature, $\rho$ is the density, $\eta$ is the viscosity, % % $c\hackscore{p}$ is the specific heat at constant pressure and $k$ is the thermal conductivity.
392
393  \subsection{Method}  % \subsection{Description}
394
395  \begin{classdesc}{TemperatureCartesian}{dom,theta=THETA,useSUPG=SUPG}  % \subsection{Method}
396  \end{classdesc}  %
397    % \begin{classdesc}{TemperatureCartesian}{dom,theta=THETA,useSUPG=SUPG}
398    % \end{classdesc}
399
400  \subsection{Benchmark Problem}  % \subsection{Benchmark Problem}
401  %===============================================================================================================  %===============================================================================================================
402
403  %=========================================================  %=========================================================

Legend:
 Removed from v.2207 changed lines Added in v.2208