Learning Objectives
1. Introduction
Nektar++ session files
1. Defining a mesh
2. Visualising the mesh
Steady diffusion problem
Unsteady diffusion problem
Unsteady advection-diffusion problem

This tutorial can be run as a Jupyter notebook. Try it on our cloud resources!

Learning Objectives

Welcome to the tutorial on the Advection-Diffusion-Reaction (ADR) solver in the Nektar++ framework. This tutorial is aimed to show the main features of the ADRSolver in a simple manner. After the completion of this tutorial, you will be familiar with:

the generation of a simple mesh in Gmsh and its conversion into a Nektar++ compatible format;
the visualisation of the mesh in Paraview or VisIt
the setup of the initial and boundary conditions, the parameters and the solver settings;
running a simulation with the ADRSolver; and
the post-processing of the data for a convergence plot and the visualisation of the results in Paraview or VisIt.

Introduction

The ADRSolver is designed to solve partial differential equations of the form

\begin{align} \alpha \frac{\partial u}{\partial t} + \lambda u + \nu \nabla u + \epsilon \nabla \cdot (D\nabla u) = f \label{eq1} \end{align}

using either continuous or discontinuous Galerkin methods. Using different combinations of the coefficients, this equation can represents various quasi-linear problems involving diffusion, unsteady advection, unsteady diffusion, unsteady advection-diffusion equation, and soforth.

In this tutorial, we focus on three specific types of this equation, through which we will demonstrate the setup and usage of Nektar++ for both steady and time-dependent problems. In particular, we will solve:

Type	Equation to solve
Steady diffusion	$ \nabla^2 u = f$
Unsteady diffusion	$ \displaystyle\frac{\partial u}{\partial t} = \epsilon \nabla^2 u $
Unsteady advection-diffusion	$ \displaystyle\frac{\partial u}{\partial t} + V \nabla u = \epsilon \nabla^2 u $

We will start with a simple steady-state diffusion problem, extend it to a unsteady-diffusion and finally will change our setup to unsteady advection-diffusion problem. For a more detailed description of this solver, please refer to the Nektar++ user guide.

Although Nektar++ is capable of 3D simulations, for simplicity in this tutorial we will consider two-dimensional case on a simple domain $ \Omega = [0,1]^2 $. Additionally to make each problem well-posed, we require known boundary conditions and a forcing function which depend on the spatial coordinates $x$ and $y$. To model these problem we create a computational domain also referred to as mesh or grid on which we apply the two-dimensional function with Dirichlet and Neumann boundary conditions.

Nektar++ session files

Any problem in Nektar++ requires a session file that defines:

a mesh that represents the domain of interest;
information on the basis to be used on each element within the domain: i.e. what basis functions to use on each element, what order they should be;
solver-specific information such as the boundary and initial conditions to be used.

This information is encoded in an XML file, in order to make it portable and structured. All session files contain a root NEKTAR block, inside which we define the information above:

<NEKTAR>
    <GEOMETRY>
        <!-- Information regarding the mesh goes in this block. -->
    </GEOMETRY>
    <EXPANSIONS>
        <!-- The basis to be used in each element. -->
    </EXPANSIONS>
    <CONDITIONS>
        <!-- Solver-specific information, including boundary and initial conditions. -->
    </CONDITIONS>
</NEKTAR>

In the rest of this tutorial, we’ll walk through how to populate some of this information using templates that you will complete to see how adjust simple parameters. However, one thing to note is that solvers can use more than a single session file. For example, if you have a mesh that never changes, this can be defined in one file; the solver conditions can then be defined in another file. Multiple files are then automatically merged at runtime.

Defining a mesh

To keep everything as simple as possible for this tutorial, to start with we will supply a mesh of $\Omega$. In Nektar++, meshes comprise of either:

quadrilateral or triangular elements in two dimensions; or
tetrahedra, triangular prisms, square-based pyramids or hexahedra in three dimensions.

Optionally, where the underlying geometry is curved, elements can be curved to align with it. If you are interested, later in this tutorial, we will guide you on how to adjust and change the mesh. The mesh file is AdrMesh.xml, which you can open up to take a look at the structure. It defines:

The geometry in a hierarchical format: elements are constructed from edges, which are constructed from vertices. These are represented in the ELEMENT, EDGE, and VERTEX blocks. Note that in this mesh, this information is compressed for space storage and is not designed to be human-readable.
Collections of elements that are called COMPOSITES. These composites can be used to group elements together and thus define boundaries or the domain. In this example, you can see that we define ```xml

T[0-34]

which tells us that triangles with IDs 0-34 make up a composite with ID 1. 

- This can then be used to define the domain:
```xml
<DOMAIN> C[1-3] </DOMAIN>

All of this will seem relatively arbitrary, so let’s have a look at the mesh before we continue further. To do this, we need to convert the mesh to an appropriate format for visualization, for which we can use the FieldConvert executable. Execute the following cell to convert the mesh:

!FieldConvert -f AdrMesh.xml AdrMesh.vtu

Visualising the mesh

In the above command, the AdrMesh.xml file is converted to a VTK file, AdrMesh.vtu. For visualisation, there are a number of options:

interactively, in the notebook: you can test this by executing the cell directly below this text;
download the files generated by the commands in this notebook, and then use software such as https://www.paraview.org or https://wci.llnl.gov/simulation/computer-codes/visit to display them locally on your computer.

In this instance we suggest executing the cell below, which runs some basic Python code to read and display the mesh. You can use the mouse scroller to zoom in and out, and click/drag to rotate.

If you click on the Z plane button , this will snap the view into two dimensions.
The Surface with edges button will show cell boundaries and allow you to see the underlying mesh.

import pyvista as pv

# First read the VTK file
mesh = pv.read('AdrMesh.vtu')

# Now create an itkwidget to visualise this in-browser.
pl = pv.PlotterITK()

# Add the mesh to the plotter.
pl.add_mesh(mesh, smooth_shading=True)

# Show the viewer.
pl.show(False)

Some important points to note in the visualisation of high-order meshes and fields:

This visualisation shows a linear representation of the mesh. That is, inside each element, we are sampling the coordinates at discrete points and then performing linear interpolation between these. Therefore lower-order simulations may look jagged when the underlying error is actually rather small.
Owing to this, in the visualisation above, you are viewing all integration points of the element. In triangular elements in particular, this can be somewhat visually odd: in Nektar++, we utilise a collapsed coordinate system to represent triangles, where two corners of a quadrilateral are merged into one. Thus the distribution of points within a triangular element is ‘squashed’ in one direction.

Steady diffusion problem

In this section we now move to looking at the solver itself, starting with the steady diffusion or Poisson equation. Physical example of the problems that can be described by this equation are steady-state heat transfer, steady Stokes flow or steady deformation. The equation takes the form

\[\nabla^2 u = f\]

where $u$ is the independent variable (for example, temperature) and $f$ is the forcing function (sometimes called source term).

To construct exact solutions to this equation, we select a non-polynomial forcing function $f(x,y) = -8\pi^2 \sin(2 \pi x) \sin(2\pi y)$, so that solutions are given by $u(x,y) = \sin(2 \pi x) \sin(2\pi y)$. These are consistent with the boundary conditions

\[\begin{array}{l} u(x_{b}, y_{b} = \pm 1) = 0,\\[1em] \dfrac{\partial }{dn} u(x_{b} = \pm 1, y_{b}) = \pm \dfrac{\partial }{dx} u(x_{b} = \pm 1, y_{b}) = \pm 2 \pi \sin(2\pi y_{b}) \end{array}\]

where $(x_{b}, y_{b})\in\partial\Omega$ represent coordinates within the boundaries of the computational domain. In the rest of this section, we will set the boundary conditions and forcing function for this solver; then, after running the solver we will post-process the data in order to visualise the results.

Setting up the session file

In this first section we will interactively setup the session file, starting from a template that has some sections missing. We will then run a simulation at polynomial order $p=4$.

**Note:** The complete condition file is also provided. If you run into problem you can have a look at the complete session file [sDiffSession-Complete.xml](sDiffSession-Complete.xml), and compare it to yours.

Setting the polynomial order

The first thing we need to set is the EXPANSIONS tag, which will allow us to specify the order of polynomial $p$ that is used within each of our elements. This is done by defining NUMMODES for each variable: the number of modes in the expansion. For a simulation of order $p$ the number of modes is one greater than $p$ and is given by $p+1$.

Open up the session file and navigate to the EXPANSIONS tag.
For each of the tags starting E, set the NUMMODES attribute so that the simulation is run at order $p=4$.

**Note:** Each of the tags starting with `E` defines an expansion on a given composite and variable name. For example ```xml ``` says that composite `C[1]` should use a $p=9$ represenation of the _modified_ polynomial basis (i.e. that defined in Karniadakis and Sherwin) for the variable $u$.

Setting the equations to be solved

Since the ADRSolver specifies many different equation forms to be solved, we need to tell it which one we’d like it to use. We also need to tell it the kind of discretisation we would like to impose: i.e. whether the underlying function space should be continuous or discontinuous. This is done by defining an equation type and a projection type.

Open up the session file and navigate to the CONDITIONS tag.
Inside the CONDITIONS tag, find the SOLVERINFO tag. This tag allows us to set different options that are passed through to the solver to define the problem. Each option is called a PROPERTY, and it has a VALUE associated with it.
To define the equation type, set the VALUE for the EQTYPE property to Poisson.
For the projection method, create a new tag inside SOLVERINFO with a PROPERTY of Projection, and set the VALUE to Continuous.

When you are done, the SOLVERINFO section should look like:

<SOLVERINFO>
    <I PROPERTY="EQTYPE"     VALUE="Poisson"     /> 
    <I PROPERTY="Projection" VALUE="Continuous"  />
</SOLVERINFO>

Setting the boundary conditions

We now need to set each of our boundaries. For the top and bottom boundaries where $y=0,1$, we require a Dirichlet condition; at the left and right sides where $x=0,1$ we will impose a Neumann condition. In both cases these will coincide with the exact solution.

Boundaries in Nektar++ are defined using two tags:

BOUNDARYREGIONS define the composites that define up a boundary region;
BOUNDARYCONDITIONS define the type and value of the boundary conditions on each boundary region.

The reason for these two definitions is to separate the concept of a geometric composite (which is just a collection of elements), and a boundary on which we need to impose some condition. This also allows multiple composites to share the same boundary condition.

In our mesh, composites 100, 200, 300 and 400 denote the bottom, right, top and left boundaries respectively. Let’s define the boundary regions, and then we’ll define the conditions themselves.

Open up the session file and navigate to the BOUNDARYREGIONS tag.
Note that the boundary regions for the left and right boundaries have been defined already. Each region is defined with a <B> tag, and each region has to have a unique ID. The text inside the tag defines the composite(s) that make up the region, so e.g. C[100] is composite 100.
Add a new boundary region for the bottom and top surfaces which has ID set to 2. The text inside the tag should be C[100,300], which denotes the combination of composites 100 and 300.

Now let’s set the boundary conditions themselves.

Open up the session file and navigate to the BOUNDARYCONDITIONS tag.
Note that each boundary condition is defined inside a REGION tag, which has an attribute REF that refers back to the ID of the BOUNDARYREGIONS tag. The condition itself is defined inside this tag.
Let’s first set the boundary conditions on the right and left, which have and ID of 0 and 1 respectively. These are Neumann conditions, defined via the tag: <N VAR="u" VALUE= />. The N tag informs Nektar++ that this is a Neumann condition, and the VAR tag that it’s for the variable $u$.
For the outwards facing normal $n$, we have that $\partial u(1,y)/\partial n = -2\pi\sin(2\pi y)$. So, for the right boundary, set the VALUE string to "-2*PI*sin(2*PI*y)". Similarly, for the left, set it to "2*PI*sin(2*PI*y)".
Finally, add the Dirichlet condition for the bottom and top region with ID of 2. To do this, duplicate one of the REGION tags for the Neumann variables, change the REF to 2, change the N to D to represent a Dirichlet condition, and set the VALUE to 0.

Finally, since the boundary conditions are set, we need to set the forcing term for the equation.

Setting the forcing function

To set the forcing function $f(x,y) =-8\pi^2 \sin(2 \pi x) \sin(2\pi y)$, we need to define it within our session file. Different solvers require different functions to be set in order to work. Where a function is required for a solver, we define it within the CONDITIONS tag with a FUNCTION. The ADRSolver expects the Forcing function to be defined.

The template is given in the XML file and you need to set the forcing value:

<FUNCTION NAME="Forcing">
    <E VAR="u" VALUE="-8*PI*PI*sin(2*PI*x)*sin(2*PI*y)" />
</FUNCTION>

Running the solver

We’re now set up to run the solver. If everything went to plan, you can execute the cell below to run the solver with your session file. Note that we pass in the mesh file first, which defines the GEOMETRY tag, and the the sDiffSession.xml file which defines the rest of the session file.

Note: If things don’t run or you get errors, have a look at the completed session file and check to see if you’ve got any differences or mistakes!

!ADRSolver AdrMesh.xml sDiffSession.xml

Post-processing

For the post-processing in the notebook environment we will use the python and vtk.js as will be shown shortly. However, do note that in the practical use of Nektar++, output files will be considerably larger; in this case we typically use ParaView or VisIt.

First, we need to convert the Nektar++ output to something readable by ParaView. As with the mesh above, we can use the FieldConvert utility. It is possible to convert the Nektar++ outputs to a number of different formats, including:

.vtu (unstructured VTK format);
.dat (ASCII Tecplot output) or .plt (binary Tecplot output).

The format output will automatically be determined via the filename supplied to FieldConvert. All three formats can be visualised uisng ParaView or VisIt, while Tecplot only supports .plt and .dat.

Execute the cell below to run FieldConvert. Note that:

-f is supplied to force output, otherwise you would be prompted whether to overwrite an existing file;
the first three arguments are the mesh, session file and field file;
the final argument is the name of the output file.

!FieldConvert -f AdrMesh.xml sDiffSession.xml AdrMesh.fld steadyDiffusion.vtu

Now we can use PyVista, combined with itkwidgets to visualise this in-browser as follows:

import pyvista as pv

# First read the VTK file
mesh = pv.read('steadyDiffusion.vtu')

# Now create an itkwidget to visualise this in-browser.
pl = pv.PlotterITK()
pl.add_mesh(mesh, scalars=mesh['u'], smooth_shading=True)
pl.show(False)

Running in parallel

Nektar++ can easily be run in parallel using MPI, as long as the solver has been installed with MPI enabled. To do this, run with mpirun in front of the execution command and provide the desired number of cores. For example:

!mpirun -np 2 ADRSolver AdrMesh.xml sDiffSession.xml

Convergence analysis

In this example we solved the laplacian equation $\nabla^2 u = f$ with forcing function $f(x,y) = -8\pi^2 \sin(2 \pi x) \sin(2\pi y)$ and the following boundary conditions:

\[\begin{array}{l} u(x_{b}, y_{b} = \pm 1) = 0,\\[1em] \dfrac{\partial }{dn} u(x_{b} = \pm 1, y_{b}) = \pm \dfrac{\partial }{dx} u(x_{b} = \pm 1, y_{b}) = \pm 2 \pi \sin(2\pi y_{b}) \end{array}\]

The exact solution to this problem that satisfies the given boundary conditions is

$u_e(x,y) = \sin(2 \pi x) \sin(2\pi y)$

We can provide this exact solution in the session file and Nekar++ reports $L^2$ and $L^\infty$ errors with respect to this exact solution at the end of the computation. The exact solution is defined in the session file inside the <CONDITIONS>...</CONDITIONS> as follows

<FUNCTION NAME="ExactSolution">
    <E VAR="u"  VALUE="sin(2*PI*x)*sin(2*PI*y)" />
</FUNCTION>

To perform convergence and error analysis, here we will conduct two sets of simulations. The first one is the h-type convergence analysis where we refine the mesh while keeping polynomila order constant. This is actually study of h-refinement. On the other hand, in the second experiment, we use a constant element size, i.e. using one mesh with a specific element size and changing the polynomial order $P$.

For simplicity and have a better control on the number of mesh and elemet sizes, a new set of structured mesh is provided with 1, 4, 9, 16, 25, 36, 49, 81 and 100 elements. These meshes are named mesh-#elmnt.xml where # is the number of element. For example, mesh-100elmnt.xml is the mesh with 100 elements.

Additionally, a new session file named errorAnalysisSession.xml is provided that is consistent with these meshes.

To perform $P$ refinement study, we need to change the polynomial order (or equivalently NUMMODES) in the session file. To make the experiment more automatic, series of session files are also provided wehre they are named exp#.xml where # is the number of NUMMODES (note NUMMODES = P + 1 ). For example exp5.xml is a session file that sets the NUMMODES=5 and contains only the followings:

<?xml version="1.0" encoding="utf-8" ?>
<NEKTAR>  
    <EXPANSIONS>     
        <E COMPOSITE="C[1]" NUMMODES="5" TYPE="MODIFIED" FIELDS="u" />
    </EXPANSIONS>
</NEKTAR>

This is the only tag defined in these session files and hence they cannot be used independently to run the simulations as they lack other necessary settings. To run the simulaitons we need to provide the desired exp#.xml session at the end of command line (after the session file) so the nummodes defined in these files overwritten the nummodes defined in the session file.

import subprocess
import re
import numpy as np
import matplotlib.pyplot as plt

############ h-type convergence tests ################################
meshList = [ "mesh-1elmnt.xml", "mesh-4elmnt.xml", "mesh-25elmnt.xml", "mesh-36elmnt.xml","mesh-49elmnt.xml","mesh-81elmnt.xml", "mesh-100elmnt.xml"]

# we will use Polynomial order 4 (NUMMODES=5) for h-type study. You can change this to any polynomial order you wish by adding the 
# associate exp sessions at the end of command line. For example, if you want to use NUMMODES=8 modify the commands as follwos
#
# FieldConvert -m dof mesh.xml session.xml exp8.xml a.stdout
#
# ADRSolver mesh.xml session.xml exp8.xml
#
# note: the name of mesh and session in the above commands should be replaced with appropriate ones you are using.


# first we need to set up the array of Ndof (number of degree of freedom)
Ndof_htype = []
for mesh in meshList:
    cmd = "FieldConvert -m dof " + mesh + " errorAnalysisSession.xml a.stdout"
    res = subprocess.check_output(cmd, shell=True, text=True)
    ndof = res.split('Total number of DOF:', maxsplit=1)[-1]\
               .split(maxsplit=1)[0]
    Ndof_htype.append(ndof)
    
Ndof_htype=np.array(Ndof_htype, dtype=int)
print("Ndof array for h type analysis:")
print(Ndof_htype)
print("Dof count are finished")
print(" ")
print("Starting simulations for h type convergence")
# now run the simulations and store the errors

L2Err_htype=[]
LinfErr_htype=[]
# now solving the equations using several polynomial orders and save L2 and Linf errors
for mesh in meshList:
    cmd="ADRSolver " + mesh + " errorAnalysisSession.xml "
    # grab the std out
    stdout = subprocess.check_output(cmd, shell=True, text=True)
    
    # get L2 error
    l2 = stdout.split('L 2 error (variable u) :', maxsplit=1)[-1]\
               .split(maxsplit=1)[0]
    L2Err_htype.append(l2)
    
    # get the Linf error
    lifn = stdout.split('L inf error (variable u) :', maxsplit=1)[-1]\
               .split(maxsplit=1)[0]
    LinfErr_htype.append(l2)

# convert the results to numeric arrays
L2Err_htype = np.array(L2Err_htype, dtype=np.float64)
LinfErr_htype = np.array(LinfErr_htype, dtype=np.float64)


print("h type convergence simulatins are finished")
print("h type L2 errors: ")
print(L2Err_htype)
print("")
print("h type Linf erros:")
print(LinfErr_htype)
print("\n")
print("Starting The P type convergence analysis")
print("\n")
############ P-type convergence tests ################################
expOrders = ["exp2.xml","exp3.xml","exp4.xml","exp5.xml","exp6.xml","exp7.xml","exp8.xml","exp9.xml","exp10.xml","exp11.xml"]

# For this experiment, we will use the mesh with 25 elements. you can change the mesh to any mesh you desire

# first set up the array of Ndof (number of degree of freedom)
Ndof_ptype=[]
for exp in expOrders:
    cmd = "FieldConvert -m dof mesh-25elmnt.xml errorAnalysisSession.xml " + exp + " a.stdout"
    res = subprocess.check_output(cmd, shell=True, text=True)
    ndof = res.split('Total number of DOF:', maxsplit=1)[-1]\
               .split(maxsplit=1)[0]
    Ndof_ptype.append(ndof)
    
Ndof_ptype=np.array(Ndof_ptype, dtype=int)
print("Ndof array for P type analysis:")
print(Ndof_ptype)
print("Dof count finished")
print(" ")
print("Starting simulations for P type convergence")

L2Err_ptype=[]
LinfErr_ptype=[]
# now solving the equations using several polynomial orders and save L2 and Linf errors
for exp in expOrders:
    cmd="ADRSolver mesh-25elmnt.xml errorAnalysisSession.xml " + exp
    # grab the std out
    stdout = subprocess.check_output(cmd, shell=True, text=True)
    
    # get L2 error
    l2 = stdout.split('L 2 error (variable u) :', maxsplit=1)[-1]\
               .split(maxsplit=1)[0]
    L2Err_ptype.append(l2)
    
    # get the Linf error
    lifn = stdout.split('L inf error (variable u) :', maxsplit=1)[-1]\
               .split(maxsplit=1)[0]
    LinfErr_ptype.append(l2)

# convert the results to numeric arrays
L2Err_ptype = np.array(L2Err_ptype, dtype=np.float64)
LinfErr_ptype = np.array(LinfErr_ptype, dtype=np.float64)

print("P type simulations are finished.")
print("P type L2 errors: ")
print(L2Err_ptype)
print("")
print("P type Linf erros:")
print(LinfErr_ptype)
print("\n\n")

# Creating figure and axes
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=[7, 11])

# plot h_type L2 error
ax1.loglog(Ndof_htype, L2Err_htype, '-or', linewidth=2, markersize=10, label='h Convergence')
# Plot p_type L2 error
ax1.loglog(Ndof_ptype, L2Err_ptype, '-sb', linewidth=2, markersize=10, label='P Convergence')

ax1.set_title('L2 Error', fontsize=15)
ax1.set_xlabel('Number of Dof', fontsize=13)
ax1.set_ylabel('||E||2', fontsize=13)
ax1.legend()
ax1.grid(True, which="both", ls="-")


# plot h_type Linf error
ax2.loglog(Ndof_htype, LinfErr_htype, '-or', linewidth=2, markersize=10, label='h Convergence')
# Plot p_type Linf error
ax2.loglog(Ndof_ptype, LinfErr_ptype, '-sb', linewidth=2, markersize=10, label='P Convergence')

ax2.set_title('Linf Error', fontsize=15)
ax2.set_xlabel('Number of Dof', fontsize=13)
ax2.set_ylabel('||E||inf', fontsize=13)
ax2.legend()
ax2.grid(True, which="both", ls="-")

plt.tight_layout()
plt.show()

Extra tasks

Note: For these tasks, you might find it helpful to create new cells to run your code in. You can insert a new Jupyter cell using the + in the toolbar above, and run the command from inside that cell. You can have as many cells as you wish. Alternatively, create a new terminal and run the commands inside there by going to File > New > Terminal.

Process the result of the parallel run and confirm the results are identical to the serial run. Create a new cell below, and use the code in the post-processing section above to visualize the results.
Note that the $L^2$ and $L^\infty$ errors are reported in the solver output, e.g.
```
 L 2 error (variable u) : 0.00153213
 L inf error (variable u) : 0.00319081
```
Change the EXPANSIONS tag to run at polynomial orders $p=3$ through $p=9$ and observe the trend in error.

Unsteady diffusion problem

In this section we will extend the previous example to an unsteady diffusion, in order to highlight the timestepping framework within Nektar++. This takes the form:

\[\frac{\partial u}{\partial t} = \epsilon \nabla^2 u\]

where $\epsilon$ is a constant diffusion coefficient. This is equipped with homogeneous Dirichlet boundary conditions

\[u(x_{b} = \pm 1, y_{b}) = u(x_{b}, y_{b} = \pm 1) = 0\]

The exact solution to this problem is:

\[u(x,y) = \sin(2\pi x) \sin(2\pi y) \exp (- 8 \pi^2 \epsilon t)\]

Changing the session file

For our starting point we consider the file uDiffSession.xml, which is the same as the completed steady diffusion session file sDiffSession-Complete.xml.

Setting the solver type and timestepping scheme

Similar to the steady diffusion, for the solver, we need to define the equation type and projection method inside SOLVERINFO tag as we did before. This section is also used to select a timestepping scheme.

Open up the session file and navigate to the SOLVERINFO tag.
Change the EQTYPE to UnsteadyDiffusion.
Add a new tag with the PROPERTY TimeIntegrationMethod. We’ll use a simple backwards Euler scheme, so set the VALUE to BackwardEuler.
Finally we need to tell the ADRSolver to use implicit timestepping to solve our diffusion problem. To do this, add a final PROPERTY called DiffusionAdvancement, for which you should set the VALUE to Implicit.

When you’re done, the SOLVERINFO section should look like the following:

<SOLVERINFO>
    <I PROPERTY="EQTYPE"                VALUE="UnsteadyDiffusion" /> 
    <I PROPERTY="Projection"            VALUE="Continuous"        />
    <I PROPERTY="TimeIntegrationScheme" VALUE="BackwardEuler"     />
    <I PROPERTY="DiffusionAdvancement"  VALUE="Implicit"          />
</SOLVERINFO>

**Note:** Although in general many different combinations of implicit/explicit timestepping are support inside `ADRSolver`, not _every_ combination is implemented. Therefore when considering problems outside the scope of this tutorial, consult the Nektar++ user guide to explore the possible options.

Adding solution parameters

For this unsteady problem we need to define several parameters for the problem, such as the timestep size, final integration time and a diffusion coefficient $\epsilon$. These parameters will be defined inside a PARAMETERS tag inside the CONDITIONS tag, where for each of them we use the following template:

<P> ParamName = ParamValue  </P>

Open up the session file and navigate to the SOLVERINFO section. Underneath SOLVERINFO, add a new PARAMETERS section, i.e.:

</SOLVERINFO> <!-- end of the SOLVERINFO section -->
<PARAMETERS>
    <!-- we'll add parameters here -->
</PARAMETERS>

Define the following parameters:
- TimeStep: Defines the timestep size, which we’ll set moderately to 0.001.
- FinTime: We’ll solve the problem for 0.3 time units: you can choose a higher value, but this will obviously lead to increased execution time!
- epsilon: This defines the physical diffusion coefficient $\epsilon$, which you should set to 0.05.

**Note:** You can define any two of `TimeStep`, `FinTime` and a third parameter `NumSteps`. Nektar++ will automatically determine the remaining parameter for you.

Let’s also define two additional parameter that control output:
- IO_InfoSteps: Controls the frequency that the solver prints out the information to the console. For example if we set IO_InfoSteps to 100, then every 100 steps you’ll see a status update on the current integration time and walltime required. By default this is set to 1.
- IO_CheckSteps: Sets the frequency that the solution will be written out, i.e. the checkpoint frequency. These files can be used to restart a simulation if it is long-running, or to visualise an animation of the unsteady solution field. For example, if we set IO_CheckSteps to 200, then every 200 steps the solution will be written in a file AdrMesh_#.chk where # is the number of that step. For example at step 400, the output is AdrMesh_2.chk. This also causes a file AdrMesh_0.chk to be output with the initial condition.

After you set all of these parameters in the <PARAMETERS> tag, your session file should look something like this:

<PARAMETERS>   
    <P> FinTime         = 0.3    </P> <!-- Simulation final time -->
    <P> TimeStep        = 0.001  </P> <!-- Time step -->
    <P> epsilon         = 0.05   </P> <!-- diffusion coefficent   -->
    <P> IO_InfoSteps    = 5     </P> <!-- Frequency of printing information for the user -->
    <P> IO_CheckSteps   = 20    </P> <!-- Frequency of writing the checkpoints, the results -->
<PARAMETERS>

Setting the boundary conditions

Similar to the steady diffusion, now we are going to set the boundary conditions. As explained in the problem description we want to set the Dirchlet boundary for all the boundaries with the value $u=0$. Armed with your knowledge from the steady diffusion section, open up your session file and set these boundary conditions. You can do this by either:

changing the two Neumann conditions to homogeneous Dirichlet conditions;
removing the Neumann conditions, and their boundary regions, and putting all of the boundary composites in the same boundary region.

Removing the forcing term

In contrast to the steady diffusion problem of the previous example, here, since the equation we are considering is

\[\frac{\partial u}{\partial t} = \epsilon \nabla^2 u\]

we do not need to set the forcing! You can open your session file and remove the Forcing function.

Setting initial conditions

However, we do need to set the initial condition. Open your session file, add a new FUNCTION with the NAME set to InitialConditions, and adjust the expression to correspond to our initial condition

\[u(x,y,t=0) = \sin(2\pi x) \sin(2\pi y)\]

The completed function should look like:

<FUNCTION NAME="InitialConditions">
    <E VAR="u"  VALUE="sin(2*PI*x)*sin(2*PI*y)" />
</FUNCTION>

Running the solver

As everything is now set you can run the solver by executing the following cell:

**Note:** If things don't run or you get errors, [have a look at the completed session file](uDiffSession-Complete.xml) and check to see if you've got any differences or mistakes!

!ADRSolver AdrMesh.xml uDiffSession.xml

Post-processing

The basics of the post-processing remain the same as we dissucssed earlier for the steady diffusion problem. However, now we have 50 output files (since the number of steps to be taken is FinTime/TiemStep = 0.3/0.0001=3000 and IO_ChekcSteps=200 so we have 3000/200=15` checkpoints). Lets try to convert all the results first.

to convert serires of the results we are going to usse a combination of shell command and FieldConvert as follows. Lets name our output files out_ud_0.vut , out_ud_1.vtu …

for n in range(16):
    !FieldConvert -f AdrMesh.xml uDiffSession.xml AdrMesh.xml AdrMesh_$n\.chk out_ud_$n\.vtu

Now that you converted all the results, you can visualise the evolution of u over the time by executing the cell below.

**Note:** You will need to stop the cell's execution below by pressing the stop button in the Jupyter interface.

import pyvista as pv
import time
import ipywidgets as widgets

# Here we read the output data. Note that in the FieldConvert command above,
# we set the output file names to be out_ud_0.vtu, out_ud_1.vtu ...
# the " 'out_ud_%d.vtu' % i " will produce the file names.
#
# In the subsequent sections, you should replace the "out_ud_" part by the 
# name you give for the outputs.
meshes = [ pv.read('out_ud_%d.vtu' % i) for i in range(1,15) ]

plotter = pv.PlotterITK(notebook=True)
plotter.add_mesh(meshes[0], scalars='u')
viewer = plotter.show(False)

# we have 15 output files
nfiles = 15

i = 1
while True:
    viewer.geometries = meshes[i-1]
    i = (i+1) % nfiles
    time.sleep(0.25)

The best to visualize the results is using Paraview, Visit or Tecplot. Since we converted our results to vtu it can be visualized with either Paraview or VisIt. We recommand you to install Paraview, then download the vtu files you converted here, Open the files in Paraview together and see the time-evolution of the solution.

**Note:** At this point, you have learned how to setup a simulation and visualize the results. you can continue to the next seb-section and practice more with the extra tasks, or you can jump to the next section `Unsteady advection-diffusion problem`. The other choice is to skip the rest of this tutorial to start the tutorial on the `Incompressible flow over cylinder` where you will learn how to do simulation using `Incompressible Navier-Stokes`.

Extra tasks

Task 1: Boundary conditions

Let’s try to make the problem a bit more complex. Previously we set all the boundary conditions to u=0. Now, lets change the bondary conditions as follows

$u=5\,(1.-x^2)$ for the top boundary
$u = 2\pi \sin(2\pi y)\exp(-\epsilon \pi t)$ for the left and right boundaries
- Note that the left and right boundaries are now time-dependent. We need to use the USERDEFINEDTYPE with the TimeDependent parameter for the variable at the boundary of interest.

open the session file and modify the boundary conditions for these three boundaries. The new conditions will look like the following:

     <REGION REF="0"> <!-- Left Boundary-->
        <D VAR="u"  USERDEFINEDTYPE="TimeDependent" 
          VALUE="2*PI*sin(2*PI*y)*exp(-epsilon*PI*t)" />
      </REGION>

     <REGION REF="1"> <!-- Right Boundary-->
        <D VAR="u"  USERDEFINEDTYPE="TimeDependent" 
           VALUE="2*PI*sin(2*PI*y)*exp(-epsilon*PI*t)" />
     </REGION>

     <REGION REF="2"> <!-- Top Boundary-->
        <D VAR="u"  VALUE="5*(1-x*x)" />
     </REGION>

     <REGION REF="3"> <!-- bottom Boundary-->
        <D VAR="u"  VALUE="0" />
     </REGION>

Note that, in the boundary region you need to separate the top and bottom boundaries. boundary ids 0 and 1 are already asing to the left and right boundaries, and till now, ID=2 is assigne collectively to top and bottom. Now ID=2 will be used for top and ID=3 for bottom:

        <BOUNDARYREGIONS>
            <B ID="2"> C[100] </B>  <!-- top  -->
            <B ID="3"> C[300] </B>  <!-- bottom   -->
        </BOUNDARYREGIONS>

Moreover, lets set the initial condition to u=0 as follows:

    <FUNCTION NAME="InitialConditions">
      <E VAR="u"  VALUE="0" />
    </FUNCTION>

Now, update these changes in the session file, run the problem and visualize the results.

Task 2: change the timestepping type

Backwards Euler is a very simple timestepping scheme that is first order. Nektar++ supports a wide range of timestepping schemes: try switching to a higher-order scheme such as the third-order diagonally-implicit Runge-Kutta scheme, by changing the TimeIntegrationMethod to DIRKOrder3. Reducing the temporal error in this manner should allow you to also validate the convergence order of the solver by adjusting NUMMODES, as in the steady-state solver.

Task 3: change to explicit timestepping

Explore the use of explicit time-advancement for this problem, by changing the DiffusionAdvancement value to Explicit. This should be combined with an appropriate TimeIntegrationMethod, such as ClassicalRungeKutta4 for a classical 4th order Runge-Kutta method. Note that you will need to considerably reduce the timestep size (by a factor of 10) – and probably you will require to also increase the IO_CheckSteps parameter to avoid unnecessary output!

Unsteady advection-diffusion problem

Finally, lets extend our practice to an unsteady advection-diffusion equation, which takes the form

\[\frac{\partial u}{\partial t} + V \nabla u = \epsilon \nabla^2 u\]

given a constant velocity $ V = (V_x, V_y) $ with periodic boundary conditions

\[\begin{array}{l} u(x_{b} = 1, y_{b}) = u(x_{b} = -1, y_{b}),\\[1em] u(x_{b}, y_{b} = 1) = u(x_{b}, y_{b} = -1) \end{array}\]

and initial conditions

\[u(x,y,t=0) = \sin(2\pi x)\sin(2\pi y)\]

whose exact solution is

\[u(x,y) = \sin \left( 2\pi (x - V_x t) \right) \sin \left( 2\pi (y - V_y t) \right) \exp (- 8 \pi^2 \epsilon t)\]

Generating a new mesh

So far, we have used a mesh that has been provided with the tutorial. For this problem, we need to have periodic boundary conditions where left and right, and top and bottom boundaries are paired respectively. Hence, it is a good practice to setup the mesh together.

In this instance, we will use Gmsh to generate a mesh and convert it to the Nektar++ format. For this, we require a .geo file to define the geometry, which is provided in this tutorial as AdvDiff.geo. In this file, the domain, element sizes and physical boundaries are defined.

First we need to make the mesh using the following command:

!gmsh -2 AdvDiff.geo

Gmsh will produce the mesh and save it with the name AdvDiff.msh. The resulting mesh AdvDiff.msh cannot be directly used with the Nektar++ solvers. We need to convert the msh format to the appropriate format for Nektar++ using the NekMesh utility. Let’s convert the mesh and save it to a temporary file AdvDiff-tmp.xml as follows:

!NekMesh AdvDiff.msh AdvDiff-tmp.xml

If you open the AdvDiff-tmp.xml file to have a look inside, you will notice that the information regarding definition of verticies and elements are in binary. However, at the end of the file you will see these lines:

        <COMPOSITE>
            <C ID="1"> T[0-34] </C>
            <C ID="2"> Q[35-58] </C>
            <C ID="3"> T[59-93] </C>
            <C ID="100"> E[46,42,0,43-44] </C>
            <C ID="200"> E[50,24,60,63,66,69,139,155] </C>
            <C ID="300"> E[153,148,114,149,151] </C>
            <C ID="400"> E[157,136,111,109,107,105,27,48] </C>
        </COMPOSITE>
        <DOMAIN> C[1-3] </DOMAIN>

Where the references for various composites and domain is given. It is easy to figure out which are the define the domain and which the boundaries. domain is defined by <DOMAIN> C[1-3] </DOMAIN>, hence, the composites with ID=1, ID=2 and ID=3 are defining the domain. The rest, i.e. ID=100,…,ID=400 defines the boundaries. Further, looking at the AdvDiff.geo file, we see will apreciate that composite ID=100 is the bottom boundary, ID=200 is the right, ID=300 is the top and finally ID=400 is the left boundary.

To be able to define the periodic boundary conditions, we need to pair the left and right boundaries as well as pair the bottom and top. This can be achieved using the peralign module for the NekMesh as follows (the -v option is for verbose so the NekMesh prints out information about the process):

!NekMesh -v -m peralign:dir=x:surf1=200:surf2=400 -m peralign:dir=y:surf1=100:surf2=300 AdvDiff.msh AdvDiffPeriodic.xml

Now we have the mesh file which is saved in AdvDiffPeriodic.xml and its left and right, and bottom and top boundaries are paired respectively. We will use this mesh for our simulation.

Changing the session file

As before, the starting point for our advection-diffusion problem is the completed unsteady diffusion problem.

Adjusting the solver settings

Since in the previous section we already set up our diffusion advancement and time-integration methods, we are pretty much good to go. We only need to define a AdvectionAdvancement to tell the ADRSolver how to integrate the advection term in time.

For the AdvectionAdvancement, only Explicit time-integration is available. However, we can combine this with the Implicit integration for the DiffusionAdvancement to have a better stability for the simulation. Having this combination requires the use of an implicit-explicit or semi-implicit timestepping scheme. In Nektar++, we can use IMEX schemes up to third order: for this example, let’s select a 2nd order scheme.

Open up the session file and change the TimeIntegrationMethod to IMEXOrder2. Your SOLVERINFO section should now look like:

<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="IMEXOrder2"                 />
</SOLVERINFO>

Solver parameters

Similar to the unsteady diffusion problem, we need to define timestep, final time and the diffusion coefficient $\epsilon$. We can use all of the same simulation settings from the previous solver.

In addition to these parameters, we need to define two more parameters that will denote a constant velocity in x and y direction, which we will call advx and advy and set to 2.0 and 1.0 respectively. Define and set the solution parameter as shown below in the session file:

<PARAMETERS>   
    <P> FinTime         = 1.0               </P>
    <P> TimeStep        = 0.005             </P>
    <P> NumSteps        = FinTime/TimeStep  </P>
    <P> IO_CheckSteps   = 20                </P>
    <P> IO_InfoSteps    = 10                </P>
    <P> epsilon         = 0.01              </P>
    <P> advx            = 2.0               </P>
    <P> advy            = 1.0               </P>
</PARAMETERS>

Setting the boundary conditions

The boundary conditions here are rather different from what we had before. We want to have periodic boundary conditions, hence we need to use a P tag for the boundaries. Further, as shown below, we gave each composite on the boundaries a unique ID, where Nektar++ solvers identify the boundaries using these IDs. First, open your session file and set the BOUNDARYREGIONS so that every edge is unique:

<BOUNDARYREGIONS>
    <B ID="0"> C[100] </B>  <!-- Bottom -->
    <B ID="1"> C[200] </B>  <!-- Right  -->
    <B ID="2"> C[300] </B>  <!-- Top    -->
    <B ID="3"> C[400] </B>  <!-- Left   -->
</BOUNDARYREGIONS>

To define a periodic boundary condition, we need to give each boundary the ID of its pair boundary. For example, for the left and right boundaries, when defining the boundary condition on the left boundary, the ID of the right boundary is given:

<REGION REF="1">
    <P VAR="u" VALUE="[3]" />
</REGION>

and for the right boundary, the ID of left boundary should be given:

<!--Left-->
<REGION REF="3">
    <P VAR="u" VALUE="[1]" />
</REGION>

We should do the same for top and bottom. After you set the boundary conditions, your session file should looks lik:

<BOUNDARYCONDITIONS>
    <!-- Bottom -->
    <REGION REF="0">
        <P VAR="u" VALUE="[2]" />
    </REGION>
    <!-- Right -->
    <REGION REF="1">
        <P VAR="u" VALUE="[3]" />
    </REGION>
    <!-- Top -->
    <REGION REF="2">
        <P VAR="u" VALUE="[0]" />
    </REGION>
    <!-- Left -->
    <REGION REF="3">
        <P VAR="u" VALUE="[1]" />
    </REGION>
</BOUNDARYCONDITIONS>

Setting the advection velocity

The ADRSolver is now expecting an advection velocity to be set that defines advection at each point in the domain. To do this, we must set an AdvectionVelocity function. Since we already parameterized their values by defining advx and advy in <PARAMETERS>, we can use these in the definition AdvectionVelocity function as shown below:

<FUNCTION NAME="AdvectionVelocity">
    <E VAR="Vx" VALUE="advx" />
    <E VAR="Vy" VALUE="advy" />
</FUNCTION>

Now you can set these velocities in the session file.

Initial conditions

We can leave the initial conditions as they were in the previous example:

<FUNCTION NAME="InitialConditions">
    <E VAR="u"  VALUE="sin(2*PI*x)*sin(2*PI*y)" />
</FUNCTION>

Exact solution

Finally, update the exact solution to reflect the advection term, corresponding to the solution:

\[u(x,y) = \sin \left( 2\pi (x - V_x t) \right) \sin \left( 2\pi (y - V_y t) \right) \exp (- 8 \pi^2 \epsilon t)\]

Running the solver

Every thing is now set, so you can run the solver by executing the following command:

!ADRSolver AdvDiffPeriodic.xml advDiffSession.xml

Post-processing

The visualizaton is exactly as we did in the previous section. Open a new cell using + in the toolbar above, convert the results and visualise them.

**Note:** You can either continue to the next seb-section and practice more with the extra tasks, or you can skip the rest of this tutorial to start the tutorial on the `Incompressible flow over cylinder` where you will learn how to do simulation using `Incompressible Navier-Stokes`.

Extra tasks

Task A: Lets try to make the problem a bit more complex. Previously we used the periodic boundary conditions for top and bottom. Lets change the top boundary condition to vary over the x and be time-dependent as well, i.e set it to $u=5\,(1.-x^2)sin(2\pi t/4)$ as follows:

      <REGION REF="2"> <!-- Top Boundary-->
        <D VAR="u" USERDEFINEDTYPE="TimeDependent" VALUE="5*(1-x*x)*sin(2*PI*t/4)" />
      </REGION>

      <REGION REF="0"> <!-- Bottop Boundary-->
        <D VAR="u" VALUE="0" />
      </REGION>

**Note:** Please observe the `USERDEFINEDTYPE=TimeDependent` parameter used to make the boundary condition unsteady.

Moreover, lets set the initial condition to u=0:

    <FUNCTION NAME="InitialConditions">
      <E VAR="u"  VALUE="0" />
    </FUNCTION>

Finally, increase the final time to FinalTime=2.0 and change the direction of y-velocity,i.e. avy=-1.0. Now, update these changes in the session file and run the problem.

A1 Before running the problem, lets clear out all the previous resutls. open a new cell (using the + sing in top menue), type and execute the following:

!rm *.chk *.vtu

A2 to run the problem, open another cell, type the following and run:

!ADRSolver AdvDiff.xml advDiffSession.xml

A3 convert all the solution files. See how many checkpoint files you have, the AdvDiff_0.chk, AdvDiff_1.chk, … . Lets assume you have n_chk file (if you have 10 output n_chk is 10 , n_chk=10 ). Further, lets assume you will name the output files as out_0.vtu, out_1.vtu, … . Open another cell and run the following command:

for n in range(16): !FieldConvert -f AdvDiffPeriodic.xml advDiffSession.xml AdvDiffPeriodic_$n\.chk out_$n\.vtu

 Now, open another cell and you can use the following to visualize your data over time:

    
    import pyvista as pv
    import time
    import ipywidgets as widgets
    
    meshes = [ pv.read('out_%d.vtu' % i) for i in range(1,15) ]
    
    plotter = pv.PlotterITK(notebook=True)
    plotter.add_mesh(meshes[0], scalars='u')
    viewer = plotter.show(False)
    
    # replace n_chk with the number of outputs you have, for example if you have
    # 50 outputs change the following line to nfiles=50
    nfiles = n_chk
    
    i = 1
    while True:
        viewer.geometries = meshes[i-1]
        i = (i+1) % nfiles
        time.sleep(0.25)
        

Task B We run our problem using the continuous Glerkin projection method. We can use Disconinuous Galerkin as well. Change the Projection in the SolverInfo to DisContinuous, run and visualize the problem.
Task C: Similar to what we did for steady-state problem and unsteady diffusion, lets increase the order of polynomial expansions by setting NUMMODES=9 as follows:
```
<EXPANSIONS>     
  <E COMPOSITE="C[1,2,3]" NUMMODES="9" TYPE="MODIFIED" FIELDS="u" />
</EXPANSIONS> 
```
Note, since the NUMMODES=9 is used, we need to decrease the timestep to 0.001.

Next Step

You have finished the tutorials on the ADR solver where you have learned basic terminology and solver settings for Nektar++. Most of these settings such as parameters, boundary conditions, initial condition, etc are common among most of the solvers. The SolverInfo settings are also very similar (maight not be exactly the same) for other Nektar++ solvers. Comprehensive information on these setting for various solvers can be found in Nektar++ user guide.

we recommend you to continue to the next step by doing the tutoiral on the **Incompressible Flow Simulation** where flow over a cylinder in 2D and Quasi-3D are solved using **Nektar++**'s **Incompressible Navier-Stokes Solver**.

Tutorials

Getting started

Theory and fundamentals

Pre and post processing

Incompressible Flow

Compressible Flow

Learning Objectives

Introduction

Nektar++ session files

Defining a mesh

Visualising the mesh

Steady diffusion problem

Setting up the session file

Setting the polynomial order

Setting the equations to be solved

Setting the boundary conditions

Setting the forcing function

Running the solver

Post-processing

Running in parallel

Convergence analysis

Extra tasks

Unsteady diffusion problem

Changing the session file

Setting the solver type and timestepping scheme

Adding solution parameters

Setting the boundary conditions

Removing the forcing term

Setting initial conditions

Running the solver

Post-processing

Extra tasks

Task 1: Boundary conditions

Task 2: change the timestepping type

Task 3: change to explicit timestepping

Unsteady advection-diffusion problem

Generating a new mesh

Changing the session file

Adjusting the solver settings

Solver parameters

Setting the boundary conditions

Setting the advection velocity

Initial conditions

Exact solution

Running the solver

Post-processing

Extra tasks

Next Step

End of the tutorial.