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
, andVERTEX
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
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
buttonwill 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$.
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 theNUMMODES
attribute so that the simulation is run at order $p=4$.
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 theSOLVERINFO
tag. This tag allows us to set different options that are passed through to the solver to define the problem. Each option is called aPROPERTY
, and it has aVALUE
associated with it. - To define the equation type, set the
VALUE
for theEQTYPE
property toPoisson
. - For the projection method, create a new tag inside
SOLVERINFO
with aPROPERTY
ofProjection
, and set theVALUE
toContinuous
.
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 beC[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 attributeREF
that refers back to theID
of theBOUNDARYREGIONS
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
of0
and1
respectively. These are Neumann conditions, defined via the tag:<N VAR="u" VALUE= />
. TheN
tag informs Nektar++ that this is a Neumann condition, and theVAR
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
of2
. To do this, duplicate one of theREGION
tags for the Neumann variables, change theREF
to2
, change theN
toD
to represent a Dirichlet condition, and set theVALUE
to0
.
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
toUnsteadyDiffusion
. - Add a new tag with the
PROPERTY
TimeIntegrationMethod
. We’ll use a simple backwards Euler scheme, so set theVALUE
toBackwardEuler
. - Finally we need to tell the
ADRSolver
to use implicit timestepping to solve our diffusion problem. To do this, add a finalPROPERTY
calledDiffusionAdvancement
, for which you should set theVALUE
toImplicit
.
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>
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. UnderneathSOLVERINFO
, add a newPARAMETERS
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 to0.001
.FinTime
: We’ll solve the problem for0.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 to0.05
.
- 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 setIO_InfoSteps
to100
, then every 100 steps you’ll see a status update on the current integration time and walltime required. By default this is set to1
.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 setIO_CheckSteps
to200
, then every 200 steps the solution will be written in a fileAdrMesh_#.chk
where#
is the number of that step. For example at step 400, the output isAdrMesh_2.chk
. This also causes a fileAdrMesh_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
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:
!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.
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.
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 theTimeDependent
parameter for the variable at the boundary of interest.
- Note that the left and right boundaries are now time-dependent. We need to use the
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.
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>
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 useDisconinuous Galerkin
as well. Change theProjection
in theSolverInfo
toDisContinuous
, 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 to0.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.