Helmholtz Solver

PIC

Helmholtz Solver

Department of Aeronautics, Imperial College London, UKScientific Computing and Imaging Institute, University of Utah, USA January 31, 2024

Nektar++ Tutorial

1 Introduction xGoals x1.1 Background x1.2 Problem description 2 Pre-processing x2.1 Configuring the expansion bases and the conditions 3 Running the solver 4 Post-processing 5 Summary x5.1 Additional Exercises

1 Introduction

[next] [prev] [prev-tail] [tail] [up] Chapterx1Introduction Welcome to the tutorial on solving the Helmholtz problem using the Advection-Diffusion-Reaction (ADR) Solver in the Nektar++ framework. This tutorial is aimed to show the main features of the ADR solver in a simple manner. If you have not already downloaded and installed Nektar++, please do so by visiting http://www.nektar.info, where you can also find the User-Guide with the instructions to install the library. This tutorial requires:

Nektar++ ADRSolver and pre- and post-processing tools,

the open-source mesh generator Gmsh,

the visualisation tool Paraview or VisIt

Goals

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 ADR solver; and

the post-processing of the data for a convergence plot and the visualisation of the results in Paraview or VisIt.

Task: 1.1 Prepare for the tutorial. Make sure that you have:

Installed and tested Nektar++ v5.5.0from a binary package, or compiled it from source. By default binary packages will install all executables in /usr/bin. If you compile from source they will be in the sub-directory dist/bin of the build directory you created in the Nektar++ source tree. We will refer to the directory containing the executables as $NEK for the remainder of the tutorial.

Downloaded the tutorial files: basics-helmholtz.tar.gz Unpack it using tar -xzvf basics-helmholtz.tar.gz to produce a directory

basics-helmholtz with subdirectories called tutorial and complete.

We will refer to the tutorial directory as $NEKTUTORIAL.

The tutorial folder contains:

a Gmsh file to generate the mesh, Helm_mesh.geo;

a .msh file containing the mesh in Gmsh format, Helm_mesh.msh;

Task: 1.2 Additionally, you should also install

a visualization package capable of reading VTK files, such as ParaView (which can be downloaded from here) or VisIt (downloaded from here). Alternatively, you can generate Tecplot formatted .dat files for use with Tecplot.

1.1 Background

The ADR solver can solve various problems, including the unsteady advection, unsteady diffusion, unsteady advection diffusion equation, etc. For a more detailed description of this solver, please refer to the User-Guide.

In this tutorial we focus on the Helmholtz equation . ∇^{2}u − λu = f,

(1.1)

where u is the independent variable. The Helmholtz equation can be solved in one, two and three spatial dimensions. We will here consider a two-dimensional problem.

1.2 Problem description

The problem we want to solve consists of known boundary conditions and forcing function which depend on x and y. To model this problem we create a computational domain also referred to as mesh or grid (see section 2) on which we apply the following two-dimensional function with Dirichlet and Neumann boundary conditions.

.  ∇^{2}u − λu = −(2π^{2} + λ) cos(πx) cos(πy), u(x, y) = cos(πx) cos(πy), u(x_{b} = ±1, y_{b}) = cos(πx_{b}) cos(πy_{b}), \frac{∂}{dn}u(x_{b}, y_{b} = ±1) = ±\frac{∂}{dy}u(x_{b}, y_{b} = ±1) = ∓π [cos(πx_{b}) sin(πy_{b})]

(1.2)

where x_{b} and y_{b} represent the boundaries of the computational domain (see section 2.1) and λ is a positive constant.

We will set the boundary conditions and forcing function for this solver (see section 2.1) then, after running the solver (see section 3) we will post-process the data in order to visualise the results (see section 4).

[next] [prev] [prev-tail] [front] [up] 2 Pre-processing

[next] [prev] [prev-tail] [tail] [up]

Chapterx2Pre-processing

The pre-processing step consists of generating the mesh in a Nektar++ compatible format. To do this we can use the open-source mesh generator Gmsh to first define the geometry, which in our case is a square mesh. The resulting mesh is shown in Fig. 2.1. The mesh file format (.msh) generated by Gmsh is not directly compatible with the Nektar++ solvers and, therefore, it needs to be converted. pict Figurex2.1: Mesh generated by Gmsh.

To do so, we need to run the Nektar++ pre-processing routine called NekMesh. This routine requires two command-line arguments: the mesh file generated by Gmsh, Helm_mesh.msh; and the name of the Nektar++-compatible mesh file that NekMesh will generate, for instance Helm_mesh.xml. Task: 2.1 Convert the .msh file into a Nektar++ input from within the $NEKTUTORIAL directory by calling

 $NEK/NekMesh Helm_mesh.msh Helm_mesh.xml

or

 $NEK/NekMesh Helm_mesh.msh Helm_mesh.xml:xml:uncompress

Note that by default the information is stored in XML blocks of compressed data to reduce the size of large meshes. The second command above tells the converter NekMesh to write the file out in uncompressed ASCII format. The generated .xml mesh file is shown below and can also be found in the completed directory. It contains 5 tags encapsulated within the GEOMETRY tag, which describes the mesh. The first tag, VERTEX, contains the spatial coordinates of the vertices of the various elements of the mesh. The second tag, EDGE contains the lines connecting the vertices. The third tag, ELEMENT, defines the elements (note that in this case we have both triangular - e.g. <T ID="0"> - as well as quadrilateral - e.g. <Q ID="85"> - elements). The fourth tag, COMPOSITE, describes the physical regions of the mesh and each composite may contain only one type of mesh entity. There are three composites composed of triangular or quadrilateral elements which represent the solution sub-domains where we want to solve the linear advection problem. We will use these composites to define expansion bases on each sub-domain in section 2.1. The composites formed by edges are the boundaries of the domain where we will prescribe suitable boundary conditions in section 2.1. Finally, the fifth tag, DOMAIN, formally specifies the overall solution domain as the union of the three element composites. For additional details on the GEOMETRY specification refer to the User-Guide.

   1 <?xmlversion="1.0"encoding="utf-8"?> 2 <NEKTAR> 3 xxx<GEOMETRYDIM="2"SPACE="2"> 4 xxxxxx<VERTEX> 5 xxxxxxxx<VID="0">2.00000000e-01-1.00000000e+000.00000000e+00</V> 6 xxxxxxxx<VID="1">5.09667784e-01-6.15240515e-010.00000000e+00</V> 7 xxxxxxxx... 8 xxxxxxxx<VID="68">-1.00000000e+001.25000000e-010.00000000e+00</V> 9 xxxxxx</VERTEX> 10 xxxxxx<EDGE> 11 xxxxxxxx<EID="0">01</E> 12 xxxxxxxx<EID="1">12</E> 13 xxxxxxxx... 14 xxxxxxxx<EID="153">4068</E> 15 xxxxxx</EDGE> 16 xxxxxx<ELEMENT> 17 xxxxxxxx<TID="0">012</T> 18 xxxxxxxx<TID="1">345</T> 19 xxxxxxxx... 20 xxxxxxxx<QID="85">14693153151</Q> 21 xxxxx</ELEMENT> 22 xxxxxx<COMPOSITE> 23 xxxxxxxx<CID="1">T[0-30]</C> 24 xxxxxxxx<CID="2">Q[62-85]</C> 25 xxxxxxxx<CID="3">T[31-61]</C> 26 xxxxxxxx<CID="100">E[46,12,20,10,45]</C> 27 xxxxxxxx<CID="200">E[50,32,108,111,114,117,87,103]</C> 28 xxxxxxxx<CID="300">E[100,64,74,66,99]</C> 29 xxxxxxxx<CID="400">E[49,33,148,150,152-153,86,104]</C> 30 xxxxxx</COMPOSITE> 31 xxxxxx<DOMAIN>C[1,2,3]</DOMAIN> 32 xxx</GEOMETRY> 33 </NEKTAR>

After generating the mesh file in the Nektar++-compatible format, Helm_mesh.xml, we can visualise the mesh. This step can be done by using the following Nektar++ built-in post-processing routine: Task: 2.2 Convert the .xml file into a .vtu format within the $NEKTUTORIAL directory by calling

 $NEK/FieldConvert Helm_mesh.xml Helm_mesh.vtu

Alternatively a tecplot .dat file can be created by changing the extension of the second file, i.e.

 $NEK/FieldConvert Helm_mesh.xml Helm_mesh.dat

This will produce a Helm_mesh.vtu file which can be directly read by the open-source visualisation tool called Paraview or VisIt. In Fig. 2.2 we show the mesh distribution for the mesh considered in this tutorial, Helm_mesh.xml.

PIC Figurex2.2: Mesh distribution with local polynomial subdivision We can now configure the conditions: initial conditions, boundary conditions, parameters and solver settings.

2.1 Configuring the expansion bases and the conditions

To set the various parameters, the solver settings and the initial and boundary conditions needed, we create a new file called Helm_conditions.xml, which can be found within the completed directory for this tutorial. This new file contains the CONDITIONS tag where we can specify the parameters of the simulations, the solver settings, the initial conditions, the boundary conditions and the exact solution and contains the EXPANSIONS tag where we can specify the polynomial order to be used inside each element of the mesh, the type of expansion bases and the type of points.

We begin to describe the Helm_conditions.xml file from the CONDITIONS tag, and in particular from the boundary conditions, initial conditions and exact solution sections:

   1 <CONDITIONS> 2 xxx... 3 xxx... 4 xxx... 5 xxx<VARIABLES> 6 xxxxxx<VID="0">u</V> 7 xxx</VARIABLES> 8 9 xxx<BOUNDARYREGIONS> 10 xxxxxx<BID="0">C[100]</B> 11 xxxxxx<BID="1">C[200]</B> 12 xxxxxx<BID="2">C[300]</B> 13 xxxxxx<BID="3">C[400]</B> 14 xxx</BOUNDARYREGIONS> 15 16 xxx<BOUNDARYCONDITIONS> 17 xxxxxx<REGIONREF="0"> 18 xxxxxxxx<NVAR="u"VALUE="-PI*cos(PI*x)*sin(PI*y)"/> 19 xxxxxx</REGION> 20 xxxxxx<REGIONREF="1"> 21 xxxxxxxx<DVAR="u"VALUE="cos(PI*x)*cos(PI*y)"/> 22 xxxxxx</REGION> 23 xxxxxx<REGIONREF="2"> 24 xxxxxxxx<NVAR="u"VALUE="PI*cos(PI*x)*sin(PI*y)"/> 25 xxxxxx</REGION> 26 xxxxxx<REGIONREF="3"> 27 xxxxxxxx<DVAR="u"VALUE="cos(PI*x)*cos(PI*y)"/> 28 xxxxxx</REGION> 29 xxx</BOUNDARYCONDITIONS> 30 31 xxx<FUNCTIONNAME="Forcing"> 32 xxxxxx<EVAR="u"VALUE="-(Lambda␣+␣2*PI*PI)*cos(PI*x)*cos(PI*y)"/> 33 xxx</FUNCTION> 34 35 xxx<FUNCTIONNAME="ExactSolution"> 36 xxxxxx<EVAR="u"VALUE="cos(PI*x)*cos(PI*y)"/> 37 xxx</FUNCTION> 38 </CONDITIONS>

In the above piece of .xml, we first need to specify the non-optional tag called VARIABLES that sets the solution variable (in this case u).

The second tag that needs to be specified is BOUNDARYREGIONS through which the user can specify the regions where to apply the boundary conditions. For instance, <B ID="0"> C[100] </B> indicates that composite 100 (which has been introduced in section 2) has a boundary ID equal to 0. This boundary ID is successively used to prescribe the boundary conditions.

The third tag is BOUNDARYCONDITIONS by which the boundary conditions are actually specified for each boundary ID specified in the BOUNDARYREGIONS tag. The syntax <D VAR="u" corresponds to a Dirichlet boundary condition for the variable u, while <N VAR="u" corresponds to Neumann boundary conditions. For additional details on the various options possible in terms of boundary conditions refer to the User-Guide.

Finally, <FUNCTION NAME="Forcing"> allows the specification of the Forcing function and <FUNCTION NAME="ExactSolution"> permits us to provide the exact solution, against which the L_{2} and L_{∞} errors are computed.

After having configured the VARIABLES tag, the boundary conditions, the forcing function and the exact solution we can complete the tag CONDITIONS prescribing the parameters necessary (PARAMETERS)and the solver settings (SOLVERINFO):

   1 <CONDITIONS> 2 xxx<PARAMETERS> 3 xxxxxx<P>Lambda=2.5</P> 4 xxx</PARAMETERS> 5 6 xxx<SOLVERINFO> 7 xxxxxx<IPROPERTY="EQTYPE"VALUE="Helmholtz"/> 8 xxxxxx<IPROPERTY="Projection"VALUE="Continuous"/> 9 xxx</SOLVERINFO> 10 xxx... 11 xxx... 12 xxx...

In the PARAMETERS tag, Lambda is the Helmholtz constant λ.

In the SOLVERINFO tag, EQTYPE is the type of equation to be solved, Projection is the spatial projection operator to be used (which in this case is specified to be ‘Continuous’) For additional solver-setting options refer to the User-Guide.

Finally, we need to specify the expansion bases we want to use in each of the three composites or sub-domains (COMPOSITE="..") introduced in section 2:

   1 <EXPANSIONS> 2 xxx<ECOMPOSITE="C[1]"NUMMODES="5"TYPE="MODIFIED"FIELDS="u"/> 3 xxx<ECOMPOSITE="C[2]"NUMMODES="5"TYPE="MODIFIED"FIELDS="u"/> 4 xxx<ECOMPOSITE="C[3]"NUMMODES="5"TYPE="MODIFIED"FIELDS="u"/> 5 </EXPANSIONS>

In particular, for all the composites, COMPOSITE="C[i]" with i=1,2,3 we select identical basis functions and polynomial order, where NUMMODES is the number of coefficients we want to use for the basis functions (that is commonly equal to P+1 where P is the polynomial order of the basis functions), TYPE allows selecting the basis functions FIELDS is the solution variable of our problem and COMPOSITE are the mesh regions created by Gmsh. For additional details on the EXPANSIONS tag refer to the User-Guide.

Task: 2.3 Generate the file Helm_conditions.xml in the directory $NEKTUTORIAL with λ = 2.5 or copy it from the completed directory.

[next] [prev] [prev-tail] [front] [up] 3 Running the solver

[next] [prev] [prev-tail] [tail] [up]

Chapterx3Running the solver

Now that we have the mesh file compatible with Nektar++ which will support periodic boundary conditions, Helm_mesh.xml, and we have completed the condition file, Helm_conditions.xml, we can run the solver by using the following command: Task: 3.1 Run the ADRSolver in the directory $NEKTUTORIAL using the command:

 $NEK/ADRSolver Helm_mesh.xml Helm_conditions.xml

Note that we have written the mesh in a separate file from the conditions. This is generally more efficient because it allows reopening just the condition file which is much smaller in size than the mesh file (especially for large problems). However, we could also have written both the mesh and the conditions in unique file and used the same command as above for running the solver (in this case with just one file instead of two as line argument).

As soon as the file finishes running, we should see the following screen output:

    =======================================================================EquationType:HelmholtzSessionName:Helm_meshSpatialDim.:2MaxSEMExp.Order:5ExpansionDim.:2ProjectionType:ContinuousGalerkinLambda:2.5Forcingfunc[0]:-(Lambda+2*PI*PI)*cos(PI*x)*cos(PI*y)=======================================================================Writing:"Helm_mesh.fld"Writing:"Helm_mesh.fld"-------------------------------------------TotalComputationTime=0s-------------------------------------------L2error(variableu):0.000159378Linferror(variableu):0.000454467

where the L2 and L inf errors are evaluated against the <FUNCTION NAME="ExactSolution"> provided in the Helm_conditions.xml file. To have a more detailed view on the solver settings and parameters used, it is possible to use the -v option (which stands for verbose) as follows: Task: 3.2 Rerun the ADRSolver with the verbose option:

 $NEK/ADRSolver -v Helm_mesh.xml Helm_conditions.xml

The simulation has now produced a final .fld binary file.

[next] [prev] [prev-tail] [front] [up] 4 Post-processing

[next] [prev] [prev-tail] [tail] [up]

Chapterx4Post-processing

Now that the simulation has been completed, we need to post-process the file in order to visualise the results. In order to do so, we can use the built-in post-processing routines within Nektar++. In particular, we can use the following command Task: 4.1 In the $NEKTUTORIAL directory convert the .xml and .chk files into a .vtu format by calling

 $NEK/FieldConvert Helm_mesh.xml Helm_conditions.xml Helm_mesh.fld
 Helm_mesh.vtu

which generates a .vtu file that is a readable format for the open-source package Paraview or VisIt. Note that we typically have to specify both the mesh .xml file and the condition .xml file. We can now open the .vtu file just generated. This produces the image in Fig.x(4.1).

PIC Figurex4.1: Helmholtz solution [next] [prev] [prev-tail] [front] [up] 5 Summary

[prev] [prev-tail] [tail] [up] Chapterx5Summary You should be now familiar with the following topics:

Generate a simple mesh in Gmsh and convert it in a Nektar++-compatible format; Visualise the mesh in Paraview; Setup the boundary conditions, the parameters and the solver settings; Run the ADR solver; and Post-process the data in order to visualise results in Paraview/VisIt.

5.1 Additional Exercises

Increase the polynomial order and plot the L_{2} error vs. the polynomial order in a semilogarithmic scale. If the solver is compiled with the MPI option, then try running the case in parallel

with mpirun -np 2.

Change the Projection Operator to DisContinuous to see the same problem running with an HDG solver.

Change one of the boundary conditiosn to a Robin (Mixed) boundary condition of the form ∂u/∂n = αu + β whith α = 1 and β determined from the exact solution.

Tip: To check the additional settings and parameters that can be used for this solver, check the folder: $NEK/solvers/ADRSolver/Tests/ where you can find several tests associated to the ADR solver.

[prev] [prev-tail] [front] [up] [next]