
FinTime = 20

TimeStep = 0.01

NumSteps = FinTime/TimeStep

\end{lstlisting}  Daniele de Grazia committed Aug 11, 2014 203   Chris Cantwell committed Sep 04, 2014 204 205 206 207 208 209 We also specify the advection velocity. We first define dummy parameters \begin{lstlisting}[style=XMLStyle]

\end{lstlisting} and then define the actual advection function as  Daniele de Grazia committed Aug 11, 2014 210 \begin{lstlisting}[style=XMLStyle]  Chris Cantwell committed Sep 04, 2014 211 212 213   Daniele de Grazia committed Aug 11, 2014 214 215 \end{lstlisting}  Chris Cantwell committed Sep 04, 2014 216 217 Two boundary regions are defined, one at each end of the domain, and periodicity is enforced  Giacomo Castiglioni committed Nov 07, 2019 218 \begin{lstlisting}[style=XMLStyle]  Chris Cantwell committed Sep 04, 2014 219 220 221 222  C[1] C[2]  Daniele de Grazia committed Aug 11, 2014 223   Chris Cantwell committed Sep 04, 2014 224 225 226 227 228 229 230 231 232 

\end{lstlisting}  Daniele de Grazia committed Aug 11, 2014 233   Chris Cantwell committed Sep 04, 2014 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] \end{lstlisting}  Daniele de Grazia committed Aug 11, 2014 244   245 \subsubsection{Running the code}  Chris Cantwell committed Sep 04, 2014 246 247 248 \begin{lstlisting}[style=BashInputStyle] ADRSolver Advection1D.xml \end{lstlisting}  Daniele de Grazia committed Aug 11, 2014 249   Chris Cantwell committed Sep 04, 2014 250 251 To visualise the output, we can convert it into either TecPlot or VTK formats \begin{lstlisting}[style=BashInputStyle]  Chris Cantwell committed Jul 14, 2015 252 253 FieldConvert Advection1D.xml Advection1D.fld Advection1D.dat FieldConvert Advection1D.xml Advection1D.fld Advection1D.vtu  Chris Cantwell committed Sep 04, 2014 254 255 \end{lstlisting}  Daniele de Grazia committed Aug 11, 2014 256   257 \subsection{2D Helmholtz Problem}  Daniele de Grazia committed Aug 11, 2014 258   Chris Cantwell committed Sep 04, 2014 259 In this example, it will be demonstrated how the Helmholtz equation can be solved on a two-dimensional domain.  Daniele de Grazia committed Aug 11, 2014 260   261 \subsubsection{Helmholtz equation}  Daniele de Grazia committed Aug 11, 2014 262 263 264 265 266 267 268 269 270  We consider the elliptic partial differential equation: \nabla^2 u + \lambda u = f where $\nabla^2$ is the Laplacian and $\lambda$ is a real positive constant.  Giacomo Castiglioni committed Nov 07, 2019 271 \subsubsection{Input file}  Daniele de Grazia committed Aug 11, 2014 272   Chris Cantwell committed Sep 04, 2014 273 274 The input for this example is given in the example file \inlsh{Helmholtz2D\_modal.xml}  Daniele de Grazia committed Aug 11, 2014 275   Chris Cantwell committed Sep 04, 2014 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 committed Aug 11, 2014 280 \begin{lstlisting}[style=XMLStyle]  Chris Cantwell committed Sep 04, 2014 281 282 283 284 285 286 287 288 289 290  Q[22-47] T[0-21] E[0-1] . . E[84,75,69,62,51,40,30,20,6] C[0-1]  Daniele de Grazia committed Aug 11, 2014 291 292 \end{lstlisting}  Chris Cantwell committed Sep 04, 2014 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 committed Aug 11, 2014 295 296 297 298 299 300 301 \begin{lstlisting}[style=XMLStyle] \end{lstlisting}  Chris Cantwell committed Sep 04, 2014 302 Only one parameter is needed for this problem. In this example $\lambda = 1$ and  Giacomo Castiglioni committed Nov 07, 2019 303 the Continuous Galerkin Method is used as projection scheme to solve the  Chris Cantwell committed Sep 04, 2014 304 305 Helmholtz equation, so we need to specify the following parameters and solver information.  Daniele de Grazia committed Aug 11, 2014 306 \begin{lstlisting}[style=XMLStyle]  Chris Cantwell committed Sep 04, 2014 307 308 309 

Lambda = 1

 Daniele de Grazia committed Aug 11, 2014 310   Chris Cantwell committed Sep 04, 2014 311   Andrea Cassinelli committed Sep 20, 2019 312 313   Chris Cantwell committed Sep 04, 2014 314   Daniele de Grazia committed Aug 11, 2014 315 316 \end{lstlisting}  Chris Cantwell committed Sep 04, 2014 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 committed Aug 11, 2014 321 322 \begin{lstlisting}[style=XMLStyle]  Chris Cantwell committed Sep 04, 2014 323 324 325 326 327 328 329 330 331 332 333  C[2] . . C[10]  Giacomo Castiglioni committed Nov 07, 2019 334   Chris Cantwell committed Sep 04, 2014 336 337   Chris Cantwell committed Sep 11, 2014 338 339   Chris Cantwell committed Sep 04, 2014 340   Daniele de Grazia committed Aug 11, 2014 341 342  . .  Chris Cantwell committed Sep 04, 2014 343   Daniele de Grazia committed Aug 11, 2014 344 345 \end{lstlisting}  Giacomo Castiglioni committed Nov 07, 2019 346 We know that for $f = -(\lambda + 2 \pi^2)sin(\pi x)cos(\pi y)$, the exact  Chris Cantwell committed Sep 04, 2014 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 committed Aug 11, 2014 350 351 \begin{lstlisting}[style=XMLStyle]  Chris Cantwell committed Sep 04, 2014 352 353   Daniele de Grazia committed Aug 11, 2014 354   Chris Cantwell committed Sep 04, 2014 355 356 357   Daniele de Grazia committed Aug 11, 2014 358 359 360 \end{lstlisting}  361 \subsubsection{Running the code}  Chris Cantwell committed Sep 04, 2014 362 \begin{lstlisting}[style=BashInputStyle]  Dave Moxey committed Apr 01, 2015 363 ADRSolver Test_Helmholtz2D_modal.xml  Chris Cantwell committed Sep 04, 2014 364 \end{lstlisting}  Daniele de Grazia committed Aug 11, 2014 365   Giacomo Castiglioni committed Nov 07, 2019 366 This execution should print out a summary of input file, the $L_2$ and  Chris Cantwell committed Sep 04, 2014 367 $L_{inf}$ errors and the time spent on the calculation.  Daniele de Grazia committed Aug 11, 2014 368   369 \subsubsection{Post-processing}  Chris Cantwell committed Sep 04, 2014 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 committed Mar 04, 2016 373 FieldConvert Helmholtz2D_modal.xml Helmholtz2D_modal.fld Helmholtz2D_modal.vtu  Chris Cantwell committed Sep 04, 2014 374 \end{lstlisting}  Chris Cantwell committed Mar 04, 2016 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 committed Aug 11, 2014 377   Chris Cantwell committed Sep 04, 2014 378 \begin{figure}  Daniele de Grazia committed Aug 11, 2014 379 \begin{center}  Chris Cantwell committed Jul 25, 2015 380 \includegraphics[width=6cm]{img/Helmholtz2D}  Daniele de Grazia committed Aug 11, 2014 381 \caption{Solution of the 2D Helmholtz Problem.}  Chris Cantwell committed Sep 04, 2014 382 \label{f:adrsolver:helmholtz2D}  Daniele de Grazia committed Aug 11, 2014 383 384 385 386 \end{center} \end{figure}  387 \subsection{Advection dominated mass transport in a pipe}  Daniele de Grazia committed Aug 11, 2014 388   Chris Cantwell committed Sep 04, 2014 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}  Chris Cantwell committed Sep 04, 2014 397 398 399 The governing equation for modelling mass transport is the unsteady advection diffusion equation: \begin{align*}  Daniele de Grazia committed Aug 11, 2014 400 \dfrac{\partial u}{\partial t} + v\nabla u + \epsilon \nabla^2 u = 0  Chris Cantwell committed Sep 04, 2014 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 committed Nov 07, 2019 420 S h(z) = \dfrac{2^{4/3}(Pe R/z)^{1/3}}{g^{1/3}\Gamma(4/3)} ,  Chris Cantwell committed Sep 04, 2014 421 \end{align*}  Giacomo Castiglioni committed Nov 07, 2019 422 where $z$ is the streamwise coordinate, $R$ the pipe radius, $\Gamma(4/3)$ an incomplete  Daniele de Grazia committed Aug 11, 2014 423 Gamma function and $Pe$ the Peclet number given by:  Chris Cantwell committed Sep 04, 2014 424 \begin{align*}  Daniele de Grazia committed Aug 11, 2014 425 Pe = \dfrac{2 U R}{\epsilon}  Chris Cantwell committed Sep 04, 2014 426 \end{align*}  Daniele de Grazia committed Aug 11, 2014 427   Chris Cantwell committed Sep 04, 2014 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 committed Aug 11, 2014 432   433 \subsubsection{Input file}  Chris Cantwell committed Sep 04, 2014 434 435 The geometry under consideration is a pipe of radius, $R = 0.5$ and length $l = 0.5$  Daniele de Grazia committed Aug 11, 2014 436 437 438  \begin{figure}[h!] \begin{center}  Chris Cantwell committed Jul 25, 2015 439 \includegraphics[width=6cm]{img/pipe}  Daniele de Grazia committed Aug 11, 2014 440 441 442 443 \caption{Pipe.} \end{center} \end{figure}  Chris Cantwell committed Sep 04, 2014 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 committed Aug 11, 2014 448   Chris Cantwell committed Sep 04, 2014 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 committed Aug 11, 2014 454 455 456 457 458 459 460 461 462 463 464 \begin{lstlisting}[style=XMLStyle] \end{lstlisting}  Chris Cantwell committed Sep 04, 2014 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 committed Aug 11, 2014 468   Chris Cantwell committed Sep 04, 2014 469 We choose to use a continuous projection and an first-order implicit-explicit  Chris Cantwell committed Sep 04, 2014 470 471 time-integration scheme. The \inltt{DiffusionAdvancement} and \inltt{AdvectionAdvancement} parameters specify how these terms are treated.  Daniele de Grazia committed Aug 11, 2014 472 \begin{lstlisting}[style=XMLStyle]  Andrea Cassinelli committed Sep 20, 2019 473 474 475 476 477 478 479 480   Daniele de Grazia committed Aug 11, 2014 481 482 \end{lstlisting}  Chris Cantwell committed Sep 04, 2014 483 We integrate for a total of $30$ time units with a time-step of $0.0005$,  Chris Cantwell committed Sep 04, 2014 484 necessary to keep the simulation numerically stable.  Daniele de Grazia committed Aug 11, 2014 485 \begin{lstlisting}[style=XMLStyle]  Chris Cantwell committed Sep 04, 2014 486 487 488 

TimeStep = 0.0005

FinalTime = 30

NumSteps = FinalTime/TimeStep

 Daniele de Grazia committed Aug 11, 2014 489 490 \end{lstlisting}  Chris Cantwell committed Sep 04, 2014 491 The value of the $\epsilon$ parameter is $\epsilon = 1/Pe$  Chris Cantwell committed Sep 04, 2014 492 493 494 \begin{lstlisting}[style=XMLStyle]

epsilon = 1.33333e-6

\end{lstlisting}  Daniele de Grazia committed Aug 11, 2014 495   Chris Cantwell committed Sep 04, 2014 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 committed Aug 11, 2014 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] C[3] C[4] C[2] C[5] \end{lstlisting}  Chris Cantwell committed Sep 04, 2014 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 committed Aug 11, 2014 532 533 534 535 536 537 538 539 \begin{lstlisting}[style=XMLStyle] \end{lstlisting}  Chris Cantwell committed Sep 04, 2014 540 We assume that the initial domain concentration is uniform everywhere and the  Giacomo Castiglioni committed Nov 07, 2019 541 same as the inlet. This is defined by,  Daniele de Grazia committed Aug 11, 2014 542 543 544 545 546 547 \begin{lstlisting}[style=XMLStyle] \end{lstlisting}  548 \subsubsection{Results}  Chris Cantwell committed Sep 04, 2014 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 committed Aug 11, 2014 553 554 555  \begin{figure}[h!] \begin{center}  Chris Cantwell committed Jul 25, 2015 556 \includegraphics[width=7cm]{img/graetz-nusselt}  Daniele de Grazia committed Aug 11, 2014 557 \caption{Concentration gradient at the surface of the pipe.}  Chris Cantwell committed Sep 04, 2014 558 \label{f:adrsolver:masstransport}  Daniele de Grazia committed Aug 11, 2014 559 560 561 \end{center} \end{figure}  Dave Moxey committed Aug 16, 2019 562 \subsection{Unsteady reaction-diffusion systems}  Daniele de Grazia committed Aug 11, 2014 563   Dave Moxey committed Aug 16, 2019 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 committed Aug 11, 2014 573   Dave Moxey committed Aug 16, 2019 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] \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] \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] BodyForce \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}.