| 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} \} |
| 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 |
| 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 |
|
|
| 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 |
| 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 |
%================================================ |
%================================================ |
| 384 |
\section{Temperature Advection Diffusion\label{TEMP ADV DIFF}} |
% \section{Temperature Advection Diffusion\label{TEMP ADV DIFF}} |
| 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 |
%========================================================= |
%========================================================= |
Parent Directory
| 
Revision Log
| 
Patch