adr.tex 22.9 KB
Newer Older
1
\chapter{Advection-Diffusion-Reaction Solver}
Daniele de Grazia's avatar
Daniele de Grazia committed
2 3 4 5 6 7

%3.4/UserGuide/Tutorial/ADRSolver
%3.4/UserGuide/Examples/ADRSolver/1DAdvection
%3.4/UserGuide/Examples/ADRSolver/3DAdvectionMassTransport
%3.4/UserGuide/Examples/ADRSolver/Helmholtz2D

8
\section{Synopsis}
Daniele de Grazia's avatar
Daniele de Grazia committed
9 10 11 12 13

The ADRSolver is designed to solve partial differential equations of the form:
\begin{equation}
\alpha \dfrac{\partial u}{\partial t} + \lambda u + \nu \nabla u + \epsilon \nabla \cdot (D \nabla u) = f
\end{equation}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
14 15
in either discontinuous or continuous projections of the solution field.
For a full list of the equations which are supported, and the capabilities of each equation,
Daniele de Grazia's avatar
Daniele de Grazia committed
16 17 18 19 20
see the table below.

\begin{table}[h!]
\begin{center}
\tiny
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
21
\renewcommand\arraystretch{2.2}
22 23 24 25 26
\begin{tabular}{llll}
\toprule
\textbf{Equation to solve} & \textbf{EquationType} & \textbf{Dimensions}   &
\textbf{Projections} \\
\midrule
Dave Moxey's avatar
Dave Moxey committed
27 28
$u = f$ &
    \inltt{Projection} & All & Continuous/Discontinuous \\
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
29
$\nabla^2 u = 0$ &
30
    \inltt{Laplace} & All &  Continuous/Discontinuous	\\
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
31
$\nabla^2 u  =  f$ &
32
    \inltt{Poisson} & All &  Continuous/Discontinuous	\\
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
33
$\nabla^2 u  + \lambda u =  f$ &
34
    \inltt{Helmholtz} & All & Continuous/Discontinuous \\
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
35
$\epsilon \nabla^2 u + \mathbf{V}\nabla u = f$ &
36
    \inltt{SteadyAdvectionDiffusion} & 2D only & Continuous/Discontinuous \\
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
37
$\epsilon \nabla^2 u +  \lambda u = f$ &
38
    \inltt{SteadyDiffusionReaction} & 2D only &  Continuous/Discontinuous \\
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
39 40
$\epsilon \nabla^2 u + \mathbf{V}\nabla u + \lambda u = f$ &
    \inltt{SteadyAdvectionDiffusionReaction} & 2D only &
41 42 43
    Continuous/Discontinuous \\
$ \dfrac{\partial u}{\partial t} + \mathbf{V}\nabla u = f$ &
    \inltt{UnsteadyAdvection} & All & Continuous/Discontinuous \\
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
44
$\dfrac{\partial u}{\partial t}  = \epsilon \nabla^2 u$ &
45
    \inltt{UnsteadyDiffusion} & All & Continuous/Discontinuous \\
Dave Moxey's avatar
Dave Moxey committed
46 47
$\dfrac{\partial u}{\partial t}  = \epsilon \nabla^2 u + R(u)$ &
    \inltt{UnsteadyReactionDiffusion} & All & Continuous \\
48 49 50 51
$\dfrac{\partial u}{\partial t}  + \mathbf{V}\nabla u = \epsilon \nabla^2 u$ &
    \inltt{UnsteadyAdvectionDiffusion} & All & Continuous/Discontinuous \\
$\dfrac{\partial u}{\partial t}  + u\nabla u =  0$ &
    \inltt{UnsteadyInviscidBurger} & 1D only & Continuous/Discontinuous \\
52
\bottomrule
Daniele de Grazia's avatar
Daniele de Grazia committed
53 54 55 56 57 58
\end{tabular}
\end{center}
\caption{Equations supported by the ADRSolver with their capabilities.}
\label{t:ADR1}
\end{table}

59
\section{Usage}
Daniele de Grazia's avatar
Daniele de Grazia committed
60

61
\begin{lstlisting}[style=BashInputStyle]
Daniele de Grazia's avatar
Daniele de Grazia committed
62
ADRSolver session.xml
63
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
64

65
\section{Session file configuration}
Daniele de Grazia's avatar
Daniele de Grazia committed
66

Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
67
The type of equation which is to be solved is specified through the EquationType
Daniele de Grazia's avatar
Daniele de Grazia committed
68 69 70
SOLVERINFO option in the session file. This can be set as in table \ref{t:ADR1}.
At present, the Steady non-symmetric solvers cannot be used in parallel. \\

71
\subsection{Solver Info}
Daniele de Grazia's avatar
Daniele de Grazia committed
72 73 74 75 76 77 78

The solver info are listed below:
\begin{itemize}
\item \textbf{Eqtype}: This sets the type of equation to solve, according to the table above.
\item \textbf{TimeIntegrationMethod}: The following types of time integration methods have been tested with each solver:
\begin{center}
\footnotesize
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
79
\renewcommand\arraystretch{1.2}
80 81
\begin{tabular}{lcccc}
\toprule
82 83
\textbf{EqType} & \textbf{Explicit} & \textbf{Diagonally Implicit} &
    \textbf{ IMEX} & \textbf{Implicit}   \\
84
\midrule
85 86
\inltt{UnsteadyAdvection} &  \checkmark 	 &	 & 	&\\
\inltt{UnsteadyDiffusion} &  \checkmark 	 & \checkmark 	 & 	&\\
Dave Moxey's avatar
Dave Moxey committed
87
\inltt{UnsteadyReactionDiffusion} &  	 &  	 & \checkmark	&\\
88 89
\inltt{UnsteadyAdvectionDiffusion} &  	 &   	& \checkmark	&\\
\inltt{UnsteadyInviscidBurger} &  \checkmark 	 &	 & 	&\\
90
\bottomrule
Daniele de Grazia's avatar
Daniele de Grazia committed
91 92
\end{tabular}
\end{center}
93

Daniele de Grazia's avatar
Daniele de Grazia committed
94
\item \textbf{Projection}: The Galerkin projection used may be either:
Daniele de Grazia's avatar
Daniele de Grazia committed
95
\begin{itemize}
96 97
	\item \inltt{Continuous} for a C0-continuous Galerkin (CG) projection.
	\item \inltt{Discontinuous} for a discontinous Galerkin (DG) projection.
Daniele de Grazia's avatar
Daniele de Grazia committed
98
\end{itemize}
Daniele de Grazia's avatar
Daniele de Grazia committed
99
\item \textbf{DiffusionAdvancement}: This specifies how to treat the diffusion term. This will be restricted by the choice of time integration scheme:
Daniele de Grazia's avatar
Daniele de Grazia committed
100
\begin{itemize}
101 102
	\item \inltt{Explicit} Requires the use of an explicit time integration
	scheme.
103
	\item \inltt{Implicit} Requires the use of a diagonally implicit, IMEX or
104
	Implicit scheme.
Daniele de Grazia's avatar
Daniele de Grazia committed
105
\end{itemize}
Daniele de Grazia's avatar
Daniele de Grazia committed
106
\item \textbf{AdvectionAdvancement}: This specifies how to treat the advection term. This will be restricted by the choice of time integration scheme:
Daniele de Grazia's avatar
Daniele de Grazia committed
107
\begin{itemize}
108 109 110
	\item \inltt{Explicit} Requires the use of an explicit or IMEX time integration
	scheme.
	\item \inltt{Implicit} Not supported at present.
Daniele de Grazia's avatar
Daniele de Grazia committed
111
\end{itemize}
Daniele de Grazia's avatar
Daniele de Grazia committed
112
\item \textbf{AdvectionType}: Specifies the type of advection:
Daniele de Grazia's avatar
Daniele de Grazia committed
113
\begin{itemize}
114 115
	\item \inltt{NonConservative} (for CG only).
	\item \inltt{WeakDG} (for DG only).
Daniele de Grazia's avatar
Daniele de Grazia committed
116 117 118
\end{itemize}
\item \textbf{DiffusionType}:
\begin{itemize}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
119
	\item \inltt{LDG} (The penalty term is proportional to an optional parameter \inltt{LDGc11} which is by default set to one; proportionality to polynomial order can be manually imposed by setting the parameter \inltt{LDGc11} equal to $p^2$).
Daniele de Grazia's avatar
Daniele de Grazia committed
120 121 122
\end{itemize}
\item \textbf{UpwindType}:
\begin{itemize}
123
	\item \inltt{Upwind}.
Daniele de Grazia's avatar
Daniele de Grazia committed
124 125 126
\end{itemize}
\end{itemize}

127
\subsection{Parameters}
Daniele de Grazia's avatar
Daniele de Grazia committed
128

129
The following parameters can be specified in the \inltt{PARAMETERS} section of
130
the session file:
Daniele de Grazia's avatar
Daniele de Grazia committed
131
\begin{itemize}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
132
\item \inltt{epsilon}: sets the diffusion coefficient $\epsilon$.\\
Daniele de Grazia's avatar
Daniele de Grazia committed
133 134
\textit{Can be used} in: SteadyDiffusionReaction, SteadyAdvectionDiffusionReaction, UnsteadyDiffusion, UnsteadyAdvectionDiffusion. \\
\textit{Default value}: 0.
135 136
\item  \inltt{d00}, \inltt{d11}, \inltt{d22}: sets the diagonal entries of the
diffusion tensor $D$. \\
Daniele de Grazia's avatar
Daniele de Grazia committed
137
\textit{Can be used in}: UnsteadyDiffusion \\
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
138
\textit{Default value}: All set to 1 (i.e. identity matrix).
139
\item  \inltt{lambda}: sets the reaction coefficient  $\lambda$. \\
Daniele de Grazia's avatar
Daniele de Grazia committed
140 141 142 143
\textit{Can be used in}: SteadyDiffusionReaction, Helmholtz, SteadyAdvectionDiffusionReaction\\
\textit{Default value}: 0.
\end{itemize}

144
\subsection{Functions}
Daniele de Grazia's avatar
Daniele de Grazia committed
145

146
The following functions can be specified inside the \inltt{CONDITIONS} section
147
of the session file:
Daniele de Grazia's avatar
Daniele de Grazia committed
148 149

\begin{itemize}
150 151 152
\item \inltt{AdvectionVelocity}: specifies the advection velocity $\mathbf{V}$.
\item \inltt{InitialConditions}: specifies the initial condition for unsteady problems.
\item \inltt{Forcing}: specifies the forcing function f.
Daniele de Grazia's avatar
Daniele de Grazia committed
153 154
\end{itemize}

155
\section{Examples}
156 157
Example files for the ADRSolver are provided in
\inltt{solvers/ADRSolver/Examples}
158

159
\subsection{1D Advection equation}
Daniele de Grazia's avatar
Daniele de Grazia committed
160

161 162
In this example, it will be demonstrated how the Advection equation can be
solved on a one-dimensional domain.
Daniele de Grazia's avatar
Daniele de Grazia committed
163

164
\subsubsection{Advection equation}
Daniele de Grazia's avatar
Daniele de Grazia committed
165 166 167 168
We consider the hyperbolic partial differential equation:
\begin{equation}
\dfrac{\partial u}{\partial t} + \dfrac{\partial f}{\partial x} = 0,
\end{equation}
169
where $f =  a u$ is the advection flux.
Daniele de Grazia's avatar
Daniele de Grazia committed
170

171
\subsubsection{Input file}
172
The input for this example is given in the example file \inlsh{Advection1D.xml}
Daniele de Grazia's avatar
Daniele de Grazia committed
173

174 175 176
The geometry section defines a 1D domain consisting of $10$ segments. On each
segment an expansion consisting of $4$ Lagrange polynomials on the
Gauss-Lobotto-Legendre points is used as specified by
Daniele de Grazia's avatar
Daniele de Grazia committed
177
\begin{lstlisting}[style=XMLStyle]
178 179 180
<EXPANSIONS>
    <E COMPOSITE="C[0]" FIELDS="u" TYPE="GLL_LAGRANGE_SEM" NUMMODES="4"/>
</EXPANSIONS>
Daniele de Grazia's avatar
Daniele de Grazia committed
181 182
\end{lstlisting}

183 184 185
Since we are solving the unsteady advection problem, we must specify this in the
solver information. We also choose to use a discontinuous flux-reconstruction
projection and use a Runge-Kutta order 4 time-integration scheme.
Daniele de Grazia's avatar
Daniele de Grazia committed
186
\begin{lstlisting}[style=XMLStyle]
187 188 189 190 191 192 193
<SOLVERINFO>
    <I PROPERTY="EQTYPE"                VALUE="UnsteadyAdvection"    />
    <I PROPERTY="Projection"            VALUE="DisContinuous"        />
    <I PROPERTY="AdvectionType"         VALUE="FRDG"                 />
    <I PROPERTY="UpwindType"            VALUE="Upwind"               />
    <I PROPERTY="TimeIntegrationMethod" VALUE="ClassicalRungeKutta4" />
</SOLVERINFO>
Daniele de Grazia's avatar
Daniele de Grazia committed
194 195
\end{lstlisting}

196 197 198 199 200 201 202
We choose to advect our solution for $20$ time units with a time-step of $0.01$
and so provide the following parameters
\begin{lstlisting}[style=XMLStyle]
<P> FinTime         = 20                   </P>
<P> TimeStep        = 0.01                 </P>
<P> NumSteps        = FinTime/TimeStep     </P>
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
203

204 205 206 207 208 209
We also specify the advection velocity. We first define dummy parameters
\begin{lstlisting}[style=XMLStyle]
<P> advx            = 1                     </P>
<P> advy            = 0                     </P>
\end{lstlisting}
and then define the actual advection function as
Daniele de Grazia's avatar
Daniele de Grazia committed
210
\begin{lstlisting}[style=XMLStyle]
211 212 213
<FUNCTION NAME="AdvectionVelocity">
    <E VAR="Vx" VALUE="advx" />
</FUNCTION>
Daniele de Grazia's avatar
Daniele de Grazia committed
214 215
\end{lstlisting}

216 217
Two boundary regions are defined, one at each end of the domain, and periodicity
is enforced
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
218
\begin{lstlisting}[style=XMLStyle]
219 220 221 222
<BOUNDARYREGIONS>
    <B ID="0"> C[1] </B>
    <B ID="1"> C[2] </B>
</BOUNDARYREGIONS>
Daniele de Grazia's avatar
Daniele de Grazia committed
223

224 225 226 227 228 229 230 231 232
<BOUNDARYCONDITIONS>
    <REGION REF="0">
        <P VAR="u" VALUE="[1]" />
    </REGION>
    <REGION REF="1">
        <P VAR="u" VALUE="[0]" />
    </REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
233

234 235 236 237 238 239 240 241 242 243
Finally, we specify the initial value of the solution on the domain
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="InitialConditions">
    <E VAR="u" VALUE="exp(-20.0*x*x)" />
</FUNCTION>

<FUNCTION NAME="ExactSolution">
    <E VAR="u" VALUE="exp(-20.0*x*x)" />
</FUNCTION>
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
244

245
\subsubsection{Running the code}
246 247 248
\begin{lstlisting}[style=BashInputStyle]
ADRSolver Advection1D.xml
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
249

250 251
To visualise the output, we can convert it into either TecPlot or VTK formats
\begin{lstlisting}[style=BashInputStyle]
252 253
FieldConvert Advection1D.xml Advection1D.fld Advection1D.dat
FieldConvert Advection1D.xml Advection1D.fld Advection1D.vtu
254 255
\end{lstlisting}

Daniele de Grazia's avatar
Daniele de Grazia committed
256

257
\subsection{2D Helmholtz Problem}
Daniele de Grazia's avatar
Daniele de Grazia committed
258

259
In this example, it will be demonstrated how the Helmholtz equation can be solved on a two-dimensional domain.
Daniele de Grazia's avatar
Daniele de Grazia committed
260

261
\subsubsection{Helmholtz equation}
Daniele de Grazia's avatar
Daniele de Grazia committed
262 263 264 265 266 267 268 269 270

We consider the elliptic partial differential equation:

\begin{equation}
\nabla^2 u  + \lambda u =  f
\end{equation}

where $\nabla^2$ is the Laplacian and $\lambda$ is a real positive constant.

Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
271
\subsubsection{Input file}
Daniele de Grazia's avatar
Daniele de Grazia committed
272

273 274
The input for this example is given in the example file
\inlsh{Helmholtz2D\_modal.xml}
Daniele de Grazia's avatar
Daniele de Grazia committed
275

276 277 278 279
The geometry for this problem is a two-dimensional octagonal plane containing
both triangles and quadrilaterals. Note that a mesh composite may only contain
one type of element. Therefore, we define two composites for the domain, while
the rest are used for enforcing boundary conditions.
Daniele de Grazia's avatar
Daniele de Grazia committed
280
\begin{lstlisting}[style=XMLStyle]
281 282 283 284 285 286 287 288 289 290
<COMPOSITE>
     <C ID="0"> Q[22-47] </C>
     <C ID="1"> T[0-21] </C>
     <C ID="2"> E[0-1] </C>
     .
     .
     <C ID="10"> E[84,75,69,62,51,40,30,20,6] </C>
</COMPOSITE>

<DOMAIN> C[0-1] </DOMAIN>
Daniele de Grazia's avatar
Daniele de Grazia committed
291 292
\end{lstlisting}

293 294
For both the triangular and quadrilateral elements, we use the modified Legendre
basis with $7$ modes (maximum polynomial order is $6$).
Daniele de Grazia's avatar
Daniele de Grazia committed
295 296 297 298 299 300 301
\begin{lstlisting}[style=XMLStyle]
<EXPANSIONS>
    <E COMPOSITE="C[0]" NUMMODES="7" FIELDS="u" TYPE="MODIFIED" />
    <E COMPOSITE="C[1]" NUMMODES="7" FIELDS="u" TYPE="MODIFIED" />
</EXPANSIONS>
\end{lstlisting}

302
Only one parameter is needed for this problem. In this example $\lambda = 1$ and
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
303
the Continuous Galerkin Method is used as projection scheme to solve the
304 305
Helmholtz equation, so we need to specify the following parameters and solver
information.
Daniele de Grazia's avatar
Daniele de Grazia committed
306
\begin{lstlisting}[style=XMLStyle]
307 308 309
<PARAMETERS>
    <P> Lambda = 1 </P>
</PARAMETERS>
Daniele de Grazia's avatar
Daniele de Grazia committed
310

311
<SOLVERINFO>
312 313
    <I PROPERTY="EQTYPE"        VALUE="Helmholtz"  />
    <I PROPERTY="Projection"    VALUE="Continuous" />
314
</SOLVERINFO>
Daniele de Grazia's avatar
Daniele de Grazia committed
315 316
\end{lstlisting}

317 318 319 320
All three basic boundary condition types have been used in this example:
Dirichlet, Neumann and Robin boundary. The boundary regions are defined, each of
which corresponds to one of the edge composites defined earlier. Each boundary
region is then assigned an appropriate boundary condition.
Daniele de Grazia's avatar
Daniele de Grazia committed
321 322
\begin{lstlisting}[style=XMLStyle]
<BOUNDARYREGIONS>
323 324 325 326 327 328 329 330 331 332 333
    <B ID="0"> C[2] </B>
    .
    .
    <B ID="8"> C[10] </B>
</BOUNDARYREGIONS>

<BOUNDARYCONDITIONS>
    <REGION REF="0">
        <D VAR="u" VALUE="sin(PI*x)*sin(PI*y)" />
    </REGION>
    <REGION REF="1">
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
334
        <R VAR="u" VALUE="sin(PI*x)*sin(PI*y)-PI*sin(PI*x)*cos(PI*y)"
335
           PRIMCOEFF="1" />
336 337
    </REGION>
    <REGION REF="2">
338 339
        <N VAR="u" VALUE="(5/sqrt(61))*PI*cos(PI*x)*sin(PI*y)-
                          (6/sqrt(61))*PI*sin(PI*x)*cos(PI*y)" />
340
    </REGION>
Daniele de Grazia's avatar
Daniele de Grazia committed
341 342
        .
        .
343
</BOUNDARYCONDITIONS>
Daniele de Grazia's avatar
Daniele de Grazia committed
344 345
\end{lstlisting}

Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
346
We know that for $f = -(\lambda + 2 \pi^2)sin(\pi x)cos(\pi y)$, the exact
347 348 349
solution of the two-dimensional Helmholtz equation is $u = sin(\pi x)cos(\pi
y)$. These functions are defined specified to initialise the problem and verify
the correct solution is obtained by evaluating the $L_2$ and $L_{inf}$ errors.
Daniele de Grazia's avatar
Daniele de Grazia committed
350 351
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="Forcing">
352 353
    <E VAR="u" VALUE="-(Lambda + 2*PI*PI)*sin(PI*x)*sin(PI*y)" />
</FUNCTION>
Daniele de Grazia's avatar
Daniele de Grazia committed
354

355 356 357
<FUNCTION NAME="ExactSolution">
    <E VAR="u" VALUE="sin(PI*x)*sin(PI*y)" />
</FUNCTION>
Daniele de Grazia's avatar
Daniele de Grazia committed
358 359 360
\end{lstlisting}


361
\subsubsection{Running the code}
362
\begin{lstlisting}[style=BashInputStyle]
363
ADRSolver Test_Helmholtz2D_modal.xml
364
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
365

Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
366
This execution should print out a summary of input file, the $L_2$ and
367
$L_{inf}$ errors and the time spent on the calculation.
Daniele de Grazia's avatar
Daniele de Grazia committed
368

369
\subsubsection{Post-processing}
370 371 372
Simulation results are written in the file Helmholtz2D\_modal.fld. We can choose
to visualise the output in Gmsh
\begin{lstlisting}[style=BashInputStyle]
Chris Cantwell's avatar
Chris Cantwell committed
373
FieldConvert Helmholtz2D_modal.xml Helmholtz2D_modal.fld Helmholtz2D_modal.vtu
374
\end{lstlisting}
Chris Cantwell's avatar
Chris Cantwell committed
375 376
which generates the file Helmholtz2D\_modal.vtu which can be visualised and is
shown in Fig.~\ref{f:adrsolver:helmholtz2D}
Daniele de Grazia's avatar
Daniele de Grazia committed
377

378
\begin{figure}
Daniele de Grazia's avatar
Daniele de Grazia committed
379
\begin{center}
380
\includegraphics[width=6cm]{img/Helmholtz2D}
Daniele de Grazia's avatar
Daniele de Grazia committed
381
\caption{Solution of the 2D Helmholtz Problem.}
382
\label{f:adrsolver:helmholtz2D}
Daniele de Grazia's avatar
Daniele de Grazia committed
383 384 385 386
\end{center}
\end{figure}


387
\subsection{Advection dominated mass transport in a pipe}
Daniele de Grazia's avatar
Daniele de Grazia committed
388

389 390 391 392 393 394 395
The following example demonstrates the application of the ADRsolver for
modelling advection dominated mass transport in a straight pipe.
Such a transport regime is encountered frequently when modelling mass transport
in arteries. This is because the diffusion coefficient of small blood borne
molecules, for example oxygen or adenosine triphosphate, is very small
$O(10^{-10})$.

396
\subsubsection{Background}
397 398 399
The governing equation for modelling mass transport is the unsteady advection
diffusion equation:
\begin{align*}
Daniele de Grazia's avatar
Daniele de Grazia committed
400
\dfrac{\partial u}{\partial t}  + v\nabla u +  \epsilon \nabla^2 u = 0
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419
\end{align*}

For small diffusion coefficient, $\epsilon$, the transport is dominated by
advection and this leads to a very fine boundary layer adjacent to the surface
which must be captured in order to get a realistic representation of the wall
mass transfer processes. This creates problems not only from a meshing
perspective, but also numerically where classical oscillations are observed in
the solution due to under-resolution of the boundary layer.

The Graetz-Nusselt solution is an analytical solution of a developing mass (or
heat) transfer boundary layer in a pipe. Previously this solution has been used
as a benchmark for the accuracy of numerical methods to capture the fine
boundary layer which develops for high Peclet number transport (the ratio of
advection to diffusion). The solution is derived based on the assumption that
the velocity field within the mass transfer boundary layer is linear i.e. the
Schmidt number (the relative thickness of the momentum to mass transfer boundary
layer) is sufficiently large. The analytical solution for the non-dimensional
mass transfer at the wall is given by:
\begin{align*}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
420
S h(z) = \dfrac{2^{4/3}(Pe R/z)^{1/3}}{g^{1/3}\Gamma(4/3)} ,
421
\end{align*}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
422
where $z$ is the streamwise coordinate, $R$ the pipe radius, $\Gamma(4/3)$ an incomplete
Daniele de Grazia's avatar
Daniele de Grazia committed
423
Gamma function and $Pe$ the Peclet number given by:
424
\begin{align*}
Daniele de Grazia's avatar
Daniele de Grazia committed
425
Pe = \dfrac{2 U R}{\epsilon}
426
\end{align*}
Daniele de Grazia's avatar
Daniele de Grazia committed
427

428 429 430 431
In the following we will numerically solver mass transport in a pipe and compare
the calculated mass transfer at the wall with the Graetz-Nusselt solution. The
Peclet number of the transport regime under consideration is $750000$, which is
physiologically relevant.
Daniele de Grazia's avatar
Daniele de Grazia committed
432

433
\subsubsection{Input file}
434 435
The geometry under consideration is a pipe of radius, $R = 0.5$ and length $l =
0.5$
Daniele de Grazia's avatar
Daniele de Grazia committed
436 437 438

\begin{figure}[h!]
\begin{center}
439
\includegraphics[width=6cm]{img/pipe}
Daniele de Grazia's avatar
Daniele de Grazia committed
440 441 442 443
\caption{Pipe.}
\end{center}
\end{figure}

444 445 446 447
Since the mass transport boundary layer will be confined to a very small layer
adjacent to the wall we do not need to mesh the interior region, hence the mesh
consists of a layer of ten prismatic elements over a thickness of 0.036R. The
elements progressively grow over the thickness of domain.
Daniele de Grazia's avatar
Daniele de Grazia committed
448

449 450 451 452 453
In this example we utilise heterogeneous polynomial order, in which the
polynomial order normal to the wall is higher so that we avoid unphysical
oscillations, and hence the incorrect solution, in the mass transport boundary
layer. To do this we specify explicitly the expansion type, points type and
distribution in each direction as follows:
Daniele de Grazia's avatar
Daniele de Grazia committed
454 455 456 457 458 459 460 461 462 463 464
\begin{lstlisting}[style=XMLStyle]
<EXPANSIONS>
  <E COMPOSITE="C[0]"
     NUMMODES="3,5,3"
     BASISTYPE="Modified_A,Modified_A,Modified_B"
     NUMPOINTS="4,6,3"
     POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre,GaussRadauMAlpha1Beta0"
     FIELDS="u" />
</EXPANSIONS>
\end{lstlisting}

465 466 467
The above represents a quadratic polynomial order in the azimuthal and
streamwise direction and 4th order polynomial normal to the wall for a prismatic
element.
Daniele de Grazia's avatar
Daniele de Grazia committed
468

469
We choose to use a continuous projection and an first-order implicit-explicit
470 471
time-integration scheme. The \inltt{DiffusionAdvancement} and
\inltt{AdvectionAdvancement} parameters specify how these terms are treated.
Daniele de Grazia's avatar
Daniele de Grazia committed
472
\begin{lstlisting}[style=XMLStyle]
473 474 475 476 477 478 479 480
<SOLVERINFO>
    <I PROPERTY="EQTYPE"                VALUE="UnsteadyAdvectionDiffusion" />
    <I PROPERTY="Projection"            VALUE="Continuous"                 />
    <I PROPERTY="DiffusionAdvancement"  VALUE="Implicit"                   />
    <I PROPERTY="AdvectionAdvancement"  VALUE="Explicit"                   />
    <I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder1"                 />
    <I PROPERTY="GlobalSysSoln"         VALUE="IterativeStaticCond"        />
</SOLVERINFO>
Daniele de Grazia's avatar
Daniele de Grazia committed
481 482
\end{lstlisting}

483
We integrate for a total of $30$ time units with a time-step of $0.0005$,
484
necessary to keep the simulation numerically stable.
Daniele de Grazia's avatar
Daniele de Grazia committed
485
\begin{lstlisting}[style=XMLStyle]
486 487 488
<P> TimeStep = 0.0005              </P>
<P> FinalTime = 30                 </P>
<P> NumSteps = FinalTime/TimeStep  </P>
Daniele de Grazia's avatar
Daniele de Grazia committed
489 490
\end{lstlisting}

491
The value of the $\epsilon$ parameter is $\epsilon = 1/Pe$
492 493 494
\begin{lstlisting}[style=XMLStyle]
<P> epsilon = 1.33333e-6           </P>
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
495

496 497 498 499 500 501 502 503 504 505
The analytical solution represents a developing mass transfer boundary layer in
a pipe. In order to reproduce this numerically we assume that the inlet
concentration is a uniform value and the outer wall concentration is zero; this
will lead to the development of the mass transport boundary layer along the
length of the pipe. Since we do not model explicitly the mass transfer in the
interior region of the pipe we assume that the inner wall surface concentration
is the same as the inlet concentration; this assumption is valid based on the
large Peclet number meaning the concentration boundary layer is confined to the
region in the immediate vicinity of the wall. The boundary conditions are
specified as follows in the input file:
Daniele de Grazia's avatar
Daniele de Grazia committed
506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529
\begin{lstlisting}[style=XMLStyle]
<BOUNDARYREGIONS>
   <B ID="0"> C[3] </B>  <!-- inlet -->
   <B ID="1"> C[4] </B>  <!-- outlet -->
   <B ID="2"> C[2] </B>  <!-- outer surface -->
   <B ID="3"> C[5] </B>  <!-- inner surface -->
</BOUNDARYREGIONS>

<BOUNDARYCONDITIONS>
  <REGION REF="0">
    <D VAR="u" VALUE="1" />
  </REGION>
  <REGION REF="1">
    <N VAR="u" VALUE="0" />
  </REGION>
  <REGION REF="2">
    <D VAR="u" VALUE="0" />
  </REGION>
  <REGION REF="3">
    <D VAR="u" VALUE="1" />
  </REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}

530 531
The velocity field within the domain is fully devqeloped pipe flow (Poiseuille
flow), hence we can define this through an analytical function as follows:
Daniele de Grazia's avatar
Daniele de Grazia committed
532 533 534 535 536 537 538 539
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="AdvectionVelocity">
  <E VAR="Vx" VALUE="0" />
  <E VAR="Vy" VALUE="0" />
  <E VAR="Vz" VALUE="2.0*(1-(x*x+y*y)/0.25)" />
</FUNCTION>
\end{lstlisting}

540
We assume that the initial domain concentration is uniform everywhere and the
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
541
same as the inlet. This is defined by,
Daniele de Grazia's avatar
Daniele de Grazia committed
542 543 544 545 546 547
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="InitialConditions">
  <E VAR="u" VALUE="1" />
</FUNCTION>
\end{lstlisting}

548
\subsubsection{Results}
549 550 551 552
To compare with the analytical expression we numerically calculate the
concentration gradient at the surface of the pipe. This is then plotted against
the analytical solution by extracting the solution along a line in the
streamwise direction, as shown in Fig.~\ref{f:adrsolver:masstransport}.
Daniele de Grazia's avatar
Daniele de Grazia committed
553 554 555

\begin{figure}[h!]
\begin{center}
556
\includegraphics[width=7cm]{img/graetz-nusselt}
Daniele de Grazia's avatar
Daniele de Grazia committed
557
\caption{Concentration gradient at the surface of the pipe.}
558
\label{f:adrsolver:masstransport}
Daniele de Grazia's avatar
Daniele de Grazia committed
559 560 561
\end{center}
\end{figure}

Dave Moxey's avatar
Dave Moxey committed
562
\subsection{Unsteady reaction-diffusion systems}
Daniele de Grazia's avatar
Daniele de Grazia committed
563

Dave Moxey's avatar
Dave Moxey committed
564 565 566 567 568 569 570 571 572
Reaction-diffusion systems are prevalent in a number of areas relating to the
modelling of various physical phenomena, and are particularly prevalent in the
study of chemical interactions and pattern formation. The ADRSolver supports the
solution of a single-variable system
\[
  \frac{\partial u}{\partial t} = \epsilon\nabla^2{u}{x} + R(u)
\]
where the diffusion coefficient $\epsilon$ and reaction term $R(u)$ are defined
using the session file.
Daniele de Grazia's avatar
Daniele de Grazia committed
573

Dave Moxey's avatar
Dave Moxey committed
574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615
\subsubsection{Numerical restrictions}

The reaction-diffusion system is only supported in a selected configuration,
which is mostly defined inside the \inltt{SOLVERINFO} block:

\begin{itemize}
  \item use of a continuous Galerkin discretisation;
  \item use an implicit-explicit (IMEX) timestepping scheme, such as
  \inltt{IMEXOrder3};
\end{itemize}

This naturally leads to the following \inltt{SOLVERINFO} configuration:
\begin{lstlisting}[style=XMLStyle]
<SOLVERINFO>
    <I PROPERTY="EQTYPE"                VALUE="UnsteadyReactionDiffusion" />
    <I PROPERTY="Projection"            VALUE="Continuous"                />
    <I PROPERTY="DiffusionAdvancement"  VALUE="Implicit"                  />
    <I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder3"                />
</SOLVERINFO>
\end{lstlisting}
Further to this, the reaction term $R(u)$ is imposed by the definition of a body
forcing function. For example, the reaction term $R(u) = 0.1u$ may be defined
using the function:
\begin{lstlisting}[style=XMLStyle]
<!-- Body force to enforce reaction term -->
<FUNCTION NAME="BodyForce">
    <E VAR="u" EVARS="u" VALUE="0.1*u" />
</FUNCTION>
\end{lstlisting}
Note in particular the use of the \inltt{EVARS} (equation variables) attribute,
which permits the definition of this function in terms of the scalar variable
$u$. This function should be used together with an appropriate \inltt{FORCING}
block (as described in section~\ref{sec:xml:forcing}):
\begin{lstlisting}[style=XMLStyle]
<FORCING>
    <FORCE TYPE="Body">
        <BODYFORCE> BodyForce </BODYFORCE>
    </FORCE>
</FORCING>
\end{lstlisting}
An example of a simple unsteady reaction-diffusion problem can be found in the
\inltt{Tests} directory in the session file \inltt{ReactionDiffusion2D.xml}.