A Programmer’s Guide to Nektar++

PIC

A Programmer’s Guide to Nektar++

Department of Aeronautics, Imperial College London, UKScientific Computing and Imaging Institute, University of Utah, USA December 9, 2019

Nektar++ Developer’s Guide

xx

Contents 1 Preface 2 Introduction x2.1 The Ethos of Nektar++ x2.2 The Structure of Nektar++ x2.3 Assumed Proficiencies x2.4 Other Software Implementations and Frameworks x2.5 How to Use This Document 3 Preliminaries x3.1 Summary of Development Tools and Best Practices x3.2 Documentation and Tutorials x3.3 Compiling Tutorials x3.4 Core Nektar++ Programming Concepts x3.5 Design Patterns x3.6 Software Testing Approaches IxxBuilding-Blocks of Our Framework (Inside the Library) 4 Inside the Library: LibUtilities x4.1 BasicConst x4.2 BasicUtils x4.3 Communication x4.4 FFT x4.5 Foundations x4.6 Interpreter x4.7 Kernel x4.8 Linear Algebra x4.9 Memory x4.10 Polylib x4.11 Time Integration 5 Inside the Library: StdRegions x5.1 The Fundamentals Behind StdRegions x5.2 The Fundamental Data Structures within StdRegions x5.3 The Fundamental Algorithms within StdRegions 6 Inside the Library: SpatialDomains x6.1 The Fundamentals Behind SpatialDomains x6.2 The Fundamental Data Structures within SpatialDomains x6.3 The Fundamental Algorithms within SpatialDomains 7 Inside the Library: LocalRegions x7.1 The Fundamentals Behind LocalRegions x7.2 The Fundamental Data Structures within LocalRegions x7.3 The Fundamental Algorithms within LocalRegions 8 Inside the Library: Collections x8.1 The Fundamentals Behind Collections x8.2 The Fundamental Data Structures within Collections x8.3 The Fundamental Algorithms within Collections 9 Inside the Library: MultiRegions x9.1 The Fundamentals Behind MultiRegions x9.2 The Fundamental Data Structures within MultiRegions x9.3 The Fundamental Algorithms within MultiRegions x9.4 Preconditioners 10 Inside the Library: GlobalMapping x10.1 The Fundamentals Behind GlobalMapping x10.2 The Fundamental Data Structures within GlobalMapping x10.3 The Fundamental Algorithms within GlobalMapping 11 Inside the Library: FieldUtils x11.1 The Fundamentals Behind FieldUtils x11.2 The Fundamental Data Structures within FieldUtils x11.3 The Fundamental Algorithms within FieldUtils 12 Inside the Library: SolverUtils x12.1 The Fundamentals Behind SolverUtils x12.2 The Fundamental Data Structures within SolverUtils x12.3 The Fundamental Algorithms within SolverUtils IIxxSolvers 13 ADRSolver: Solving the Advection-Reaction-Diffusion Equation 14 IncNavierStokesSolver: Solving the Incompressible Navier-Stokes Equations 15 CompressibleFlowSolver: Solving the Compressible Navier-Stokes Equations x15.1 Fundamental Theories of CompressibleFlowSolver x15.2 Data Structure of CompressibleFlowSolver x15.3 Flow Chart of CompressibleFlowSolver IIIxxUtilities 16 FieldConvert 17 NekMesh IVxxNekPy: Python interface to Nektar++ 18 Introduction x18.1 Features and functionality 19 Installing NekPy x19.1 Compiling and installing Nektar++ x19.2 Using the bindings 20 Package structure 21 NekPy wrapping guide x21.1 Defining a library x21.2 Basic class wrapping 22 Documentation 23 Memory management in NekPy x23.1 Memory management in C++ and Python x23.2 Passing C++ memory to Python 24 FieldConvert in NekPy x24.1 Idea and motivation x24.2 Design and implementation x24.3 Further development and improvement Bibliography Bibliography Contents

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

1 Preface 2 Introduction x2.1 The Ethos of Nektar++ x2.2 The Structure of Nektar++ x2.3 Assumed Proficiencies x2.4 Other Software Implementations and Frameworks x2.5 How to Use This Document 3 Preliminaries x3.1 Summary of Development Tools and Best Practices xx3.1.1 General Resources xx3.1.2 Version Control (git) xx3.1.3 Building (CMake) xx3.1.4 Testing (CTest) xx3.1.5 Merge Requests (Gitlab) x3.2 Documentation and Tutorials xx3.2.1 Dependencies xx3.2.2 Compiling the User Guide xx3.2.3 Developers Guide xx3.2.4 Compiling the code documentation x3.3 Compiling Tutorials x3.4 Core Nektar++ Programming Concepts xx3.4.1 Namespaces xx3.4.2 C++ Standard Template Library (STL) xx3.4.3 typedefs xx3.4.4 Forward Declarations xx3.4.5 Templated Classes and Specialization xx3.4.6 Multiple Inheritance and the virtual Keyword xx3.4.7 Virtual Functions and Inheritance xx3.4.8 Const keyword xx3.4.9 Function pointers and bind xx3.4.10 Memory Pools and NekArray x3.5 Design Patterns xx3.5.1 Template pattern xx3.5.2 Abstract Factory Pattern x3.6 Software Testing Approaches xx3.6.1 Unit Tests xx3.6.2 Integration, System and Regression Tests xx3.6.3 Continuous Integration IxxBuilding-Blocks of Our Framework (Inside the Library) 4 Inside the Library: LibUtilities x4.1 BasicConst x4.2 BasicUtils x4.3 Communication x4.4 FFT x4.5 Foundations xx4.5.1 Points xx4.5.2 Basis x4.6 Interpreter x4.7 Kernel x4.8 Linear Algebra x4.9 Memory x4.10 Polylib x4.11 Time Integration xx4.11.1 General Linear Methods xx4.11.2 Introduction xx4.11.3 Formulation xx4.11.4 Matrix notation xx4.11.5 General Linear Methods in Nektar++ xx4.11.6 Types of time Integration Schemes xx4.11.7 Usage xx4.11.8 Implementation of a time-dependent problem xx4.11.9 Strongly imposed essential boundary conditions xx4.11.10 How to add a new time-stepping method xx4.11.11 Examples of already implemented time stepping schemes 5 Inside the Library: StdRegions x5.1 The Fundamentals Behind StdRegions xx5.1.1 Reference Element Transformations That Facilitate Separability xx5.1.2 Reference Elements On Primitive Geometric Types x5.2 The Fundamental Data Structures within StdRegions xx5.2.1 Variables at the Level of StdExpansion xx5.2.2 Variables at the Level of StdExpansion$D for various Dimensions xx5.2.3 Variables at the Level of Shape-Specific StdExpansions xx5.2.4 General Layout of the Basis Functions in Memory xx5.2.5 General Layout xx5.2.6 2D Geometries xx5.2.7 3D Geometries x5.3 The Fundamental Algorithms within StdRegions 6 Inside the Library: SpatialDomains x6.1 The Fundamentals Behind SpatialDomains xx6.1.1 Vertex and Curvature Mapping Information xx6.1.2 Regular Mappings: Geometric Factors xx6.1.3 Deformed Mappings: Geometric Factors x6.2 The Fundamental Data Structures within SpatialDomains xx6.2.1 Variables at the Level of Geometry xx6.2.2 Variables at the Level of Geometry$D for various Dimensions xx6.2.3 Variables at the Level of Shape-Specific Geometry Information xx6.2.4 Reference to World-Space Mapping xx6.2.5 MeshGraph x6.3 The Fundamental Algorithms within SpatialDomains 7 Inside the Library: LocalRegions x7.1 The Fundamentals Behind LocalRegions x7.2 The Fundamental Data Structures within LocalRegions xx7.2.1 Variables at the Level of Expansion xx7.2.2 Variables at the Level of Expansion$D for various Dimensions xx7.2.3 Variables at the Level of Shape-Specific Expansions x7.3 The Fundamental Algorithms within LocalRegions 8 Inside the Library: Collections x8.1 The Fundamentals Behind Collections x8.2 The Fundamental Data Structures within Collections x8.3 The Fundamental Algorithms within Collections 9 Inside the Library: MultiRegions x9.1 The Fundamentals Behind MultiRegions x9.2 The Fundamental Data Structures within MultiRegions xx9.2.1 Variables at the Level of ExpList xx9.2.2 Variables at the Level of ExpList$D for various Dimensions xx9.2.3 Variables at the Level of Discontinuous Field Expansions xx9.2.4 Variables at the Level of Continuous Field Expansions x9.3 The Fundamental Algorithms within MultiRegions x9.4 Preconditioners xx9.4.1 Mathematical formulation xx9.4.2 Preconditioners xx9.4.3 Low energy 10 Inside the Library: GlobalMapping x10.1 The Fundamentals Behind GlobalMapping xx10.1.1 Divergence operator xx10.1.2 Laplacian operator xx10.1.3 Anisotropic diffusion xx10.1.4 Anisotropic Laplacian operator x10.2 The Fundamental Data Structures within GlobalMapping x10.3 The Fundamental Algorithms within GlobalMapping 11 Inside the Library: FieldUtils x11.1 The Fundamentals Behind FieldUtils x11.2 The Fundamental Data Structures within FieldUtils x11.3 The Fundamental Algorithms within FieldUtils 12 Inside the Library: SolverUtils x12.1 The Fundamentals Behind SolverUtils x12.2 The Fundamental Data Structures within SolverUtils x12.3 The Fundamental Algorithms within SolverUtils IIxxSolvers 13 ADRSolver: Solving the Advection-Reaction-Diffusion Equation 14 IncNavierStokesSolver: Solving the Incompressible Navier-Stokes Equations 15 CompressibleFlowSolver: Solving the Compressible Navier-Stokes Equations x15.1 Fundamental Theories of CompressibleFlowSolver x15.2 Data Structure of CompressibleFlowSolver x15.3 Flow Chart of CompressibleFlowSolver IIIxxUtilities 16 FieldConvert 17 NekMesh IVxxNekPy: Python interface to Nektar++ 18 Introduction x18.1 Features and functionality 19 Installing NekPy x19.1 Compiling and installing Nektar++ xx19.1.1 macOS xx19.1.2 Linux: Ubuntu/Debian xx19.1.3 Compiling the wrappers x19.2 Using the bindings xx19.2.1 Examples 20 Package structure 21 NekPy wrapping guide x21.1 Defining a library x21.2 Basic class wrapping

xx21.2.1 py::class_<>

xx21.2.2 Wrapping member functions xx21.2.3 Inheritance xx21.2.4 Wrapping enums 22 Documentation 23 Memory management in NekPy x23.1 Memory management in C++ and Python xx23.1.1 C++ x23.2 Passing C++ memory to Python xx23.2.1 Passing Python data to C++ 24 FieldConvert in NekPy x24.1 Idea and motivation x24.2 Design and implementation

xx24.2.1 FieldConverter class

xx24.2.2 User workflow xx24.2.3 Conversion process x24.3 Further development and improvement Bibliography [next] [prev] [prev-tail] [front] [up] 1 Preface

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

Chapterx1Preface

Like with any software project, people want to know it origins: what motivated it, what and who drove it, and what constrained it. Nektar++ officially started as a project idea in 2004 in Salt Lake City, UT, USA, and we registered the first commit to SVN in 2006. The basic backstory is as follows: Mike Kirby (University of Utah) and Spencer Sherwin (Imperial College London) had both studied under George Karniadakis. Though Mike and Spencer did not overlap in terms of their studies (Spencer was a Princeton graduate while Mike was a Brown graduate), they both worked on Nektar, a spectral/hp element code supervised by George. George’s research group has had a long history of involvement in various software projects, e.g. PRISM (which has continued its existence under Professor Hugh Blackburn at Monash University) and Nektar. Spencer and George initiated the Nektar code with triangular (2D) and tetrahedral (3D) spectral/hp elements in the early 1990s, and the code grew and evolved with additions by Dr. Igor Lomtev, Dr. Tim Warburton, Dr. Mike Kirby, etc., all under the PhD advisor direction of George. Spencer graduated from Princeton under George’s supervision in 1995 and went on to Imperial College London as a faculty member; Mike graduated from Brown under George’s supervision in 2002 and went on to the University of Utah in 2002. In 2004, Mike and Spencer teamed up to re-write Nektar in light of modern C++ programming practices and in light of what had been learned in the decade or so since its foundation.

What were the observations that motivated Nektar++? The first observation was that Nektar’s origin was triangular and tetrahedral spectral/hp elements as applied to incompressible fluid mechanics problems (i.e., the incompressible Navier-Stokes equations). These ideas were extended to the compressible Navier-Stokes by Lomtev and the two parallel paths joined and extended into what is often called Hybrid Nektar by Tim Warburton. Hybrid Nektar (referred to as Nektar from here on out) used the C++ programming paradigm to facilitate hybrid elements: triangles and quadrilaterals in 2D and tetrahedra, hexahedra, prisms and pyramids in 3D. In addition, Warburton structured the code in a way that allowed extension to the Arbitrary Lagrangian Eulerian (ALE) formulation within Nektar, as well as various other features such as the Nektar Magneto-Hydrodynamics (MHD) solver. Mike was mentored (as a student) by Warburton and continued the expansion of Nektar (e.g. compressible ALE solver, fluid-structure interaction capabilities, etc.). The expansion of Nektar’s capabilities, under the direction of George Karniadakis, continues to this day (e.g., Stress-Nektar). The upside of this expansion was that Nektar could be expanded and used for solving more and more engineering problems; the downside, in our opinion, was that continually expanding and extending Nektar without re-evaluating its fundamental design meant that some components became quite cumbersome from the programming perspective.

The second observation that occurred was that spectral/hp elements, as used within the incompressible fluid mechanics world, could be viewed as a special case of the broader set of high-order finite element methods as applied to various fluid and solid mechanics systems. In fact, in a much wider context, these methods represent ways of discretizing various partial differential equations (PDEs) using their weak (variational) form. From this vantage point, one can see the commonality between strong-form methods such as spectral collocation and flux reconstruction, and weak-form methods such as traditional finite elements, spectral elements, and even discontinuous Galerkin methods. During the 1990s, it became more and more apparent that there was a broader context and a broader community for discussing and disseminating the ideas surrounding high-order finite elements, and correspondingly a fruitful environment for cross-pollination of ideas and programming practices.

With all this in mind, Mike and Spencer set out to re-architect Nektar from the ground up, and in homage to its C++ core, this new software suite was called Nektar++. Nektar++ is an open-source software framework designed to support the development of high-performance scalable solvers for partial differential equations using the spectral/hp element method. High-order methods are gaining prominence in several engineering and biomedical applications due to their improved accuracy over low-order techniques at reduced computational cost for a given number of degrees of freedom. However, their proliferation is often limited by their complexity, which makes these methods challenging to implement and use. Nektar++ is an initiative to overcome this limitation by encapsulating the mathematical complexities of the underlying method within an efficient C++ framework, making the techniques more accessible to the broader scientific and industrial communities. Given the commonalities and connections between various strong-form and weak-form methods and their implementations, we use the term spectral/hpmethods to refer to the entire family of methods, from traditional FEM to discontinuous Galerkin (dG) to Flux Reconstruction (FR) and beyond.

The software supports a variety of discretization techniques and implementation strategies, supporting methods research as well as application-focused computation, and the multi-layered structure of the framework allows the user to embrace as much or as little of the complexity as they need. The libraries capture the mathematical constructs of spectral/hp element methods, while the associated collection of pre-written PDE solvers provides out-of-the-box application-level functionality and a template for users who wish to develop solutions for addressing questions in their own scientific domains.

After five years of laying the groundwork for Nektar++, Dr. Chris Cantwell and Dr. David Moxey joined Spencer’s group at ICL. Their involvement in the project has greatly impacted the structure and capabilities of Nektar++, and has played a significant role in its success. In 2017, we formalized our Development Team structure. Kirby and Sherwin as Founders, along with Cantwell and Moxey, form the key Project Leaders. We established the following roles within the Nektar++ community:

User: Individuals or teams who use Nektar++ as part of their research and who may interact with the community through the mailing list but do not directly contribute code.

Contributor: Individuals or teams who use Nektar++ as part of their work but also contribute modifications back into the code which arise as a direct consequence of their research.

Developer: Individuals who use Nektar++ for their research, and make code contributions which not only benefit their own research goals but also benefits the wider needs of the Nektar++ community. Such contributions typically benefit multiple application domains, and developers will make the extra effort to generalize new functionality beyond their own needs. They also fix bugs, identified by others, in areas of the code with which they are familiar.

Senior Developer: Senior Developers are involved in the development of Nektar++ beyond their individual research area and interact in more of a transcendent way, making contributions widely across the codebase. Senior developers are entrusted with the tasks of reviewing and merging contributions made by others and maintaining the integrity of the code.

Project Leader: These individuals meet all the requirements of Senior Developers but in addition direct how Nektar++ evolves in terms of applications, solvers, library and educational outreach.

This developer’s guide, like most of Nektar++, is not due to one sole person but an army of people each with different talents, skills and motivations. As we conclude this preface, we will now provide a short biography for the four editors of this volume, and also provide a listing of all the people who have contributed to this document.

Robert M. (Mike) Kirby : Prof. Mike Kirby is a Professor and the Associate Director of the School of Computing at the University of Utah.

Spencer J. Sherwin

: Prof. Spencer Sherwin is a Professor of Computational Fluid Mechanics in the Department of Aeronautics at Imperial College London (UK).

Chris D. Cantwell

: Dr. Chris Cantwell is a Senior Lecturer in Aeronautics in the Department of Aeronautics at Imperial College London (UK).

David Moxey

: Dr. David Moxey is a Senior Lecturer in Engineering in the College of Engineering, Mathematics and Physical Sciences at the University of Exeter (UK).

Contributors

: The thank all the various participants in the Nektar++ Project who have contributed to this

document1 1Affiliation information is in reference to where the person was employed at the time of their contribution. : Dr. Peter Vos (ICL) and Mr. Dav de St. Germain (Utah).

To see the full team and other current details about Nektar++, visit our website.

We hope that you find this developer’s guide of use to you, that you find Nektar++ helpful to your educational or industrial pursuits, and that like us, you will actively want to engage and contribute back to the Nektar++ community.

Mike Kirby Spencer Sherwin Chris Cantwell David Moxey [next] [prev] [prev-tail] [front] [up] 2 Introduction

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

Chapterx2Introduction

Nektar++ [15], is an open-source software framework designed to support the development of high-performance scalable solvers for partial differential equations using the spectral/hp element method. What in part makes Nektar++ unique is that it is an initiative to overcome the mathematical complexities of the methods by encapsulating them within an efficient cross-platform C++ software environment, consequently making the techniques more accessible to the broader scientific and industrial communities. The software supports a variety of discretization techniques and implementation strategies [61,x17,x16,x12], supporting methods research as well as application-focused computation, and the multi-layered structure of the framework allows the user to embrace as much or as little of the complexity as they need. The libraries capture the mathematical constructs of spectral/hp element methods in the form of a hierarchy of C++ components, while the associated collection of pre-written PDE solvers provides out-of-the-box functionality across a range of application areas, as well as a template for users who wish to develop solutions for addressing questions in their own scientific domains.

The cross-platform nature of the software libraries enables the rapid development of solvers for use in a wide variety of computing environments. The code accommodates both small research problems, suitable for desktop computers, and large-scale industrial simulations, requiring modern HPC infrastructure, where there is a need to maintain efficiency and scalability up to many thousands of processor cores.

Nektar++ provides a single codebase with the following key features:

Arbitrary-order spectral/hp element discretisations in one, two and three dimensions;

Support for variable polynomial order in space and heterogeneous polynomial order within two- and three-dimensional elements;

High-order representation of the geometry;

Continuous Galerkin, discontinuous Galerkin and hybridised discontinuous Galerkin projections; Support for a Fourier extension of the spectral element mesh;

Support for a range of linear solvers and preconditioners;

Multiple implementation strategies for achieving linear algebra performance on a range of platforms;

Efficient parallel communication using MPI showing strong scaling up to 8192-cores on Archer, the UK national HPC system;

A range of time integration schemes implemented using generalised linear methods; and

Cross-platform support for Linux, OS X and Windows operating systems.

In addition to the core functionality listed above, Nektar++ includes a number of solvers covering a range of application areas. A range of pre-processing and post-processing utilities are also included with support for popular mesh and visualization formats, and an extensive test suite ensures the robustness of the core functionality.

x2.1 The Ethos of Nektar++

x2.2 The Structure of Nektar++

x2.3 Assumed Proficiencies

x2.4 Other Software Implementations and Frameworks

x2.5 How to Use This Document

The Ethos of Nektar++

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

2.1 The Ethos of Nektar++ As with any research effort, one is required to decide on a set of guiding principles that will drive the investigation. Similarly with a software development effort of this form, we spent considerable time early on considering which aspects we wanted to be distinctive about Nektar++ and also what guiding principles would be most suitable to support both the goals and the boundaries of what we wanted to do. We did this for at least two reasons: (1) We acknowledged then and now that there are various software packages and open-source efforts that deal with finite element frameworks, and so we wanted to be able to understand and express to people those things we thought were distinctive to us – that is, our “selling points”. (2) We also acknowledged, from our own experience on software projects, that if we did not set up some collection of guiding principles for our work, that we would gravitate towards trying to be “all things to all men", and in doing so be at odds with the first item. Below are a list of the guiding principles, the “ethos”, of the Nektar++ software development effort.

Principles:

The following are our three guiding principles for Nektar++ which respect (1) above:

Efficiently: Nektar++ was to be a “true” high-order code. “True” is put in quotations because we acknowledge that high-order means different things to different communities. Based upon a review of the literature, we came to the conclusion that part of our h-to-p philosophy should be that we accommodate polynomial degrees ranging from zero (finite volumes) or one (traditional linear finite elements) up to what is considered “spectral” (pseudospectral) orders of 16^{th} degree. As part of our early work [61], we established that in order to span this range of polynomial degrees and attempt to maintain some level of computational efficiency, we would need to develop order-aware algorithms: that is, we would need to utilize different (equivalent) algorithms appropriate for a particular order. This principle was the starting point of our h-to-p efficiently branding and continues to be a driving principle of our work.

Transparently: Nektar++ was to be agnostic as to what the “right” way to discretize a partial differential equation (PDE). “Right” is put in quotations to acknowledge that like the issue of polynomial order, there appear to be different “camps” who hold very strong views as to which discretization method should be used. Some might concede that continuous Galerkin methods are very natural for elliptic and parabolic PDEs and then work very laboriously to shoehorn all mixed-type PDEs (such as the incompressible Navier-Stokes equations) into this form. This holds true (similarly) for finite volume and dG proponents with regards to hyperbolic problems (and those systems that are dominated by hyperbolicity). Our goal was to be a high-order finite element framework that allowed users to experiment with continuous Galerkin (cG) and discontinuous Galerkin (dG) methods. After our original developments, we incorporated Flux-Reconstruction Methods (FR Methods) into our framework; this confirmed that generality of our approach as our underlying library components merely needed new features added to accommodate the FR perspective on dG methods. As stated earlier, this is in part why we use the term spectral/hp as it allows us to capture all these various methods (e.g., cG, dG and FR) under one common element-and-order terminology.

Seamlessly: Nektar++ was to be a code that could run from the desktop (or laptop) to exascale, seamlessly. It was our observation that many grand challenge efforts target petascale and exascale computing, with the idea that in the future what is done now on a supercomputer will be done on the desktop. The supercomputer of today is in a rack within five years and on our desktop in ten years. However, we also acknowledged that many of our end users were interested in computing now. That is, following from their engineering tradition they had an engineering or science problem to solve, and they wanted the ability to run on machines ranging from their laptops to exascale as the problem demanded. Not all problems fit on your laptop, and yet not all problems are exascale problems. Like with our h-to-p efficiently approach, we wanted our algorithms and code development to allow a range of choices in the hands of the engineer.

Based upon these principles, our long-term vision is a software framework that allows the engineer the flexibility to make all these critical choices (elemental discretization, polynomial order, discretization methodology, algorithms to employ, and architecture-specific details) while at the same time having smart and adaptive defaults. That is, we want to accommodate both the engineer that wants control over all the various knobs that must be set to actualize a simulation run and the engineer that wants to remain at the high-level and let the software system choose what is best in terms of discretization, polynomial order, algorithm choice, etc. We want the flexibility to self tune, while targeting auto-tuning.

The three items above denote the principles of our development efforts. We now present the “guardrails”. By guardrails, we mean the guidelines we use to help steer us along towards our goal as stated above. These are not meant to be absolutes necessarily, but are considerations we often use to try to keep us on track in terms of respecting our three guiding principles above.

Guardrails:

We are principally concerned with advection-reaction-diffusion and conservation law problems. There are many simulation codes that are deal with solid mechanics, and we did not see ourselves as being a competitor with them. As such, we did not initially construct our framework to inherently deal with H − div and H − curl spaces (and their respective element types). Our perspective was to hold a system of scalar fields and, as needed, constrain them to respect certain mathematical properties.

We rely heavily on the ability to do tensor-product-based quadrature. Although our more recent additions to the Point and Basis information within our library allows for unstructured quadrature, many of our algorithms are designed to try to capitalize on tensor-product properties.

As we said earlier, we have designed our code to run at a range of orders, from one to sixteen. However, we acknowledge that for many simulation regimes, we often do not run at the lowest and/or the highest orders. For many applications, the sweet spot seems to be polynomials of degree four through eight. At the present time, a reasonable amount of time and effort has been spent optimizing for this range.

The robustness of our testing is often dictated by the application areas that have funded our code development. At the present time, our incompressible Navier-Stokes solver is probably the most tested of our solvers, followed by our ADR solver. We make every attempt to test our codes thoroughly; however we benefit greatly from various users “stress-testing” our codes and providing feedback on things we can improve (or better yet, suggesting coding solutions).

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

The Structure of Nektar++

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

2.2 The Structure of Nektar++

When we speak of Nektar++, it often means different things to different people, all under an umbrella of code. For the founders of Nektar++, the view was that the core of Nektar++ was the library, and that everything else was to be built around or on top of our basic spectral/hp framework. For others, the heart of Nektar++ are its solvers – the collection of simulation codes, built upon the library, that enable users to solve science and engineering problems. For a smaller group of people, it is all the add-ons that we provide in our utilities, from our mesh generation techniques to our visualization software ideas.

For the purposes of this document, we will structure our discussions into three main parts: library functionality, solvers and utilities. In this section, we will provide an overview of the structure of the libraries. We will start by giving a quick overview of the basic subdirectories contained within library and their purpose. We will then provide a bottom-up description of how the library can be viewed, as well as a top-down perspective. Each perspective (bottom-up or top-down) is fully consistent with each other; the advantage of these approaches is that they help the future developer understand the library as someone trying to build up towards our solvers, or conversely someone trying to understand our solver functionality having already been a solver user and now trying to understand the library components on which it was built.

The basic subdirectories with the library are as follows:

LibUtilities: This library contains all the basic mathematical and computer science building blocks of the Nektar++ code.

StdRegions: This library contains the objects that express “standard region” data and operations. In one dimension, this is the StdSegment. In two dimensions, this is the StdTri (Triangle) and StdQuad (Quadrilateral). In three dimensions, this is the StdTet (Tetrahedra), StdHex (Hexahedra), StdPrism (Prism) and StdPyr (Pyramid). These represent the seven different standardized reference regions over which we support differentiation and integration.

SpatialDomains: This library contains the mesh and elemental geometric information. In particular, this part of the library deals with the basic mesh data structures, and the mapping information (such as Jacobians) from StdRegions to LocalRegions.

LocalRegions: This library contains objects that express data and operations on individual physical elements of the computational mesh. Local regions are spectral/hp elements in world-space (either straight-/planar-sided or curved-sided). Using C++ terminology, a local region object is-a standard region object and has-a spatial domain object.

MultiRegions: This library holds the data structures that represent sets of elements (local regions). At the most fundamental level, these represent the union of local regions into a (geometrically-contiguous) space. By specifying the interaction of these elements, the function space they represent and/or the approximation method is defined. It is at this point in the hierarchy that a set of (local region) elements can be thought of as representing a dG or cG field.

Collections: In this library we amalgamate, in a linear algebra sense, the action of key operators on multiple (standard region or local region) elements into a single, memory-efficient block. These strategies depend on external factors such as BLAS implementation and the geometry of interest.

GlobalMapping: This library supports the analytical mapping of complex physical domains to simpler computational domains.

NekMeshUtils: This library contains processing modules for the generation (potentially from CAD geometries), conversion and manipulation of high-order meshes.

SolverUtils: This library contains data structures and algorithms which form the basis of solvers, or provide auxiliary functionality.

UnitTests: This part of the library contains unit tests that allow us to verify the correctness of the core functionality within Nektar++. These are useful for verifying that new additions or modifications to the lower-levels of the code do not compromise existing functionality or correctness.

Bottom-Up Perspective

The bottom-up perspective on the library is best understood from Figures 2.1 – 2.2. In Figure 2.1, we take the view of understanding the geometric regions over which we build approximations. Our starting point is within the StdRegions library, in which we define our canonical standard regions. There are seven fundamental regions supported by Nektar++: segments (1D), triangles and quadrilaterals (2D) and tetrahedra, hexahedra, prisms and pyramids (3D). Since we principally employ Gaussian quadrature, these regions are defined by various tensor-product and collapsing of the compact interval [−1, 1]. For the purposes of illustration, let us use a quadrilateral as our example. The StdQuadExp is a region defined on [−1, 1] × [−1, 1] over which we can build approximations ϕ^{e}(ξ_{1}, ξ_{2}). Typically ϕ^{e} is based upon polynomials in each coordinate direction; using linear functions in both directions yields the traditional Q(1) space in traditional finite elements.

pict

Figurex2.1: Diagram showing the LocalRegion (yellow) is-a StdRegion (blue) which has-a SpatialDomain (orange) in terms of coordinate systems and mapping functions. This diagram highlights how points (positions in space) are mapped between different (geometric) regions.

Since ϕ^{e}(ξ_{1}, ξ_{2}) lives on Q = [−1, 1] × [−1, 1] (i.e. ϕ^{e} : Q → ℝ) and is for this example polynomial, we can integrate it exactly (to machine precision) using Gaussian integration, and we can differentiate it by writing it in a Lagrange basis and forming a differential operator matrix to act on values of the function evaluated at points. If the function were not polynomial but instead only a smooth function, we could approximate it with quadrature and decide an appropriate basis by which to approximate its derivatives. All the routines needed for differentiating and integrating polynomials over various standard regions are contained within the StdRegions directory (and will be discussed in Chapter 5). A local region expansion, such as a basis defined on a quadrilateral element, QuadExp, is a linear combination of basis functions over its corresponding standard region as mapped by information contained within its spatial domain mapping. Local region class definitions are in the LocalRegions directory (and will be discussed in Chapter 7). Using the inheritance language of C++, we would say that a local region is-a standard region and has-a spatial domain object. The SpatialDomains directory contains information that expresses the mapping function χ_{e}(·) from the standard region to a local region. SpatialDomains is explored in Chapter 6. In the case of our quad example, the SpatialDomain object held by a QuadExp would connect the local region to its StdRegion parent, and correspondingly would allow integration and differentiation in world space (i.e., the natural coordinates in which the local expansion lives). If E denotes our geometric region in world space and if F : E → ℝ is built upon polynomials over its standard region, then we obtain F (x_{1}, x_{2}) = ϕ^{e}(χ_{e}^{−1} (x_{1}, x_{2})). Note that even though ϕ^{e} is polynomial and χ_{e} is polynomial, the composition using the inverse of χ_{e} is not guaranteed to be polynomial: it is only guaranteed to be a smooth function.

Putting this in the context of MultiRegions, we arrive at Figure 2.2. The MultiRegions directory (which will be discussed in Chapter 9) contains various data structures that combine local regions. One can think of a multi-region as being a set of local region objects in which some collection of geometric and/or function properties are enforced. Conceptually, one can have a set of local region objects that have no relationship to each other in space. This is a set in the mathematical sense, but not really meaningful to us for solving approximation properties. Most often we want to think of sets of local regions as being collections of elements that are geometrically contiguous – that is, given any two elements in the set, we expect that there exists a path that allows us to trace from one element to the next. Assuming a geometrically continuous collection of elements, we can now ask if the functions built over those elements form a piece-wise discontinuous approximation of a function over our Multiregion or a piece-wise continuous C^{0} approximation of a function over our Multiregion.

PIC

Figurex2.2: Diagram showing how a MultiRegion (green) contains a collection of LocalRegions (yellow), where a LocalRegion is-a StdRegion (blue) which has-a SpatialDomain (orange) in terms of coordinate systems and mapping functions. This diagram highlights the expansions of the solution formed over each geometric domain.

Top-Down Perspective

The top-down perspective on the library is best understood from Figure 2.3. From this perspective, we are interested in understanding Nektar++ from the solvers various people have contributed.

PIC

Figurex2.3: Diagram showing the top-down perspective on library components. Various solvers (along the top) can capitalize on generic SolverUtils and are built upon the core Nektar++ libraries consisting of MultiRegions (green), Collections (red), LocalRegions (yellow), SpatialDomains (orange) and StdRegions (blue). The most generic mathematical and computer science features are contained within LibUtilities, which in part draws from various community resources such as Boost, Metis, etc.

Assumed Proficiencies

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

2.3 Assumed Proficiencies

This developer’s guide is designed for the experienced spectral/hp user who wishes to go beyond using various Nektar++ solvers and possibly to add new features or capabilities at the library, solver or utilities levels. Since the focus of this document is Nektar++, we cannot recapitulate all relevant mathematical or computer science concepts upon which our framework is built. In this section we provide a listing of areas and/or topics of assumed knowledge, and we provide a non-exhaustive list of references to help the reader see the general areas of additional reading they may need to benefit fully from this manual.

We assume the reader has a familiarity and comfort-level with the following areas:

Finite Element Methods (FEM) x[40,x54,x9] and more generally the mathematical ideas surrounding continuous Galerkin (cG) methods. This includes basic calculus of variation concepts, basic partial differential equation knowledge, and general forms of discretization and approximation.

Polynomial Methodsx[19,x32,x36], and in particular concepts surrounding polynomial spaces, basis functions and numerical differentiation and quadrature.

Spectral Element Methods (SEM)x[25,x44].

Discontinuous Galerkin Methods [20,x39].

Scientific Computingx[35,x43]. We assume a graduate-level proficiency in basic computational techniques such as dealing with numerical precision (accuracy and conditioning), root-finding algorithms, differentiation and integration, and basic optimization.

Linear Algebra [59,x24]. We assume a strong knowledge of linear algebra and a reasonable comfort level with the various numerical algorithms used for the solution of symmetric and non-symmetric linear systems.

Other Software Implementations and Frameworks

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

2.4 Other Software Implementations and Frameworks

In the last ten years a collection of software frameworks has been put forward to try to bridge the gap between the mathematics of high-order methods and their implementation. A number of software packages already exist for fluid dynamics which implement high-order finite element methods, although these packages are typically targeted at a specific domain or provide limited high-order capabilities as an extension. A major challenge many practitioners have with spectral/hp elements and high-order methods, in general, is the complexity (in terms of algorithmic design) they encounter. In this section, we give an incomplete but representative summary of several of these attempts to overcome this challenge.

The Nektar flow solver is the predecessor to Nektar++ and implements the spectral/hp element method for solving the incompressible and compressible Navier-Stokes equations in both 2D and 3D. While it is widely used and the implementation is computationally efficient on small parallel problems, achieving scaling on large HPC clusters is challenging. Semtex [11] implements the 2D spectral element method coupled with a Fourier expansion in the third direction. The implementation is highly efficient, but can only be parallelised through Fourier-mode decomposition. Nek5000 [29] is a 3D spectral element code, based on hexahedral elements, which has been used for massively parallel simulations up to 300,000 cores. The Non-hydrostatic Unified Model of the Atmosphere (NUMA) [33] is a spectral element framework that employs continuous and discontinuous Galerkin strategies for solving a particular problem of interest, but in a way on which others could adopt and build. Hermes [60] implements hp-FEM for two-dimensional problems and has been used in a number of application areas. Limited high-order finite element capabilities are also included in a number of general purpose PDE frameworks including the DUNE project [23] and deal.II [10]. FEniCS [47] is a collaborative project for the development of scientific computing tools, with a particular focus on the automated solution of differential equations by finite element methods (FEM). Through the use of concepts such as meta-programming, FEniCS tries to keep the solving of PDEs with FEM, from the application programmers’ perspective, as close to the mathematical expressions as possible without sacrificing computational efficiency. A number of codes also implement high-order finite element methods on GPGPUs including nudg++, which implments a nodal discontinuous Galerkin scheme [37], and PyFR [64], which supports a range of flux reconstruction techniques.

How to Use This Document

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

2.5 How to Use This Document

In the next chapter, we will introduce the reader to various computer science tools and ideas upon which we rely heavily in our development and deployment of Nektar++. The remainder of this developer’s guide is then partitioned into three parts.

In Part I, we provide an overview of the various data structures and algorithms that live within the library subdirectory. We think of the library as containing all the basic building blocks of the Nektar++ code, as discussed above. In this part of the developer’s guide, we have dedicated a chapter to each subdirectory within the library. Each chapter contains three main sections. The first section of a chapter provides an overview of the mathematical concepts and terminology used within the chapter. This is not meant to be a detailed tutorial, but rather a reminder of basic concepts. The second section of a chapter provides a detailed description of the data structures introduced in that part of the library. The third section of a chapter is dedicated to explaining the algorithmic aspects of the library. Instead of going function by function or method by method, we have decided to structure this section in the style of “frequently asked questions" (FAQ). Based upon our long experience with students, postdocs and collaborators, we have distilled down a collection of questions that we will use (from the pedagogical perspective) to allow us to drill down into key algorithmic aspects of the library.

In Part II, we provide an overview of the many solvers implemented on top of the Nektar++ library. Each chapter is dedicated to a different solver, and correspondingly may take on a slightly different style of presentation to match the depth of mathematical, computer science and engineering knowledge to understand the basic data structures and algorithms represented there.

In Part III, we provide an overview of many of the utilities often used in our simulation pipeline – from things that aid the user in the preprocessing steps such as setting up parameter files and meshing to the postprocessing step of field conversion for visualization.

How you engage with this material is based upon your goals. If your goal is to:

Introduce a new solver, then you would want to jump to Part II and see how other solvers have been built. Then, depending on what library functions are needed, you may need to step back into the library parts of this manual to understand how to use some of our basic library functionality.

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

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

x2.4 Other Software Implementations and Frameworks x2.5 How to Use This Document

Summary of Development Tools and Best Practices

[next] [prev] [prev-tail] [tail] [up] 3.1 Summary of Development Tools and Best Practices

3.1.1 General Resources

Mailing List Please make sure that you sign up to the Nektar++ mailing list. This can be done at the following URL: https://mailman.ic.ac.uk/mailman/listinfo/nektar-users While the name suggests it is for users of Nektar++, this list is also used to ask questions about Nektar++ development.

Blog

The Nektar++ blog (https://www.nektar.info/cat/community) provides a broad range of posts on topics such as compiling the code on specific machines, to discussions of Nektar++ in specific application areas, to recently published papers which have made use of the code.

Contributing to this resource is a valuable way to support the project, as well as promoting your own work.

Annual Workshop

In 2015, we held the first Nektar++ workshop, which was a great success and followed by a similar event in 2016. It is now an annual event and allows first-hand access to the core Nektar++ development team as well as a range of other Nektar++ users.

3.1.2 Version Control (git)

The Nektar++ code is managed in a distributed version control system called git. To obtain the code directly from the repository and add new content to the repository requires the git command-line tool (or a suitable GUI equivalent). This is available for Linux, Mac OSX and Windows.

If you plan to work with the Nektar++ community, you will need to have a reasonable understanding of the git software. While it is beyond the scope of this document to discuss how to use git, it is important for someone new to git to spend time understanding how the tool works. For this purpose, we highly recommend familiarizing yourself with it using any of the many online resources (such as https://git-scm.com).

Anonymous Access

Nektar++ may be cloned anonymously using the following command:

git clone http://gitlab.nektar.info/clone/nektar/nektar.git nektar++

Your local copy of the code can be updated to include the latest changes in the repository using the command:

git pull

The Anonymous Access approach will provide you with a complete copy of the code, but you will be unable to push changes or new developments back to the repository. For this, you should use Collaborative Access.

Collaborative Access

Once you are familiar with git, have introduced yourself to the development community (see the mailing list information above), and are ready to become a contributing developer, you will need to register an account on the Nektar++ Gitlab website:

https://gitlab.nektar.info.

To use use authenticated access, using your new account, you must first upload the public part of your SSH key to your Gitlab profile. Generating and managing SSH keys is beyond the scope of this document. However, on Linux and OSX, one generally finds these keys in the

\$HOME/.ssh directory.

To upload the key, visit: https://gitlab.nektar.info/profile, select SSH keys from the menu on the left and follow the instructions.

Registering for an account and uploading your SSH key need only be done once.

To clone Nektar++, use the git command:

git clone git@gitlab.nektar.info:nektar/nektar.git nektar++

Note the different URL to use authenticated access, in comparison to the anonymous access.

Managing and Contributing code

Code contribution then follows three basic steps:

Create an issue to describe the code updates you are making;

Branch the Nektar++ master branch and make your changes on that branch; and

Submit a merge request on the Nektar++ Gitlab website that your updates are ready to be merged into the master branch.

More details regarding the concepts mentioned above is found below. The information is also available in the CONTRIBUTING.md file in the root directory of the code, or online at: https://gitlab.nektar.info/nektar/nektar/blob/master/CONTRIBUTING.md

Issues - The initial step for those who wish to add code to the master repository is to create an issue ticket that describes the defect, bug, proposed additions, changes, updates, etc. This is done on the Gitlab website at:

https://gitlab.nektar.info/nektar/nektar/issues

Please ensure you provide sufficient detail when creating the issue to cover all of the following (as required): Describe what is being requested, why it is important / necessary, an initial list of files that may be effected, any potential problems the change/addition may cause, and any other information that will help the development team understand the request and alert others to your work to avoid duplication of effort.

Branches - The second step is to create a git branch in which to do the actual code development. In your local copy of the Nektar++ repository, run the command:

   gitxbranchx<branch-name>

replacing <branch-name> with a suitable name for your branch. The naming convention used for branches reflects the nature of the change and is composed of a prefix and descriptor, separated by a forward slash. Prefixes are one of:

feature: used for developments which constitute a new capability in the code;

fix: used for changes to fix a bug in the code;

tidy: used for changes which improve the quality or readability of the code but do not change its function;

ticket: can be used to reference changes which address a specific issue.

Please choose concise descriptors in all-lowercase, using hyphens to separate words. An

example <branch-name> might therefore be: feature/low-energy-preconditioner

Now make your new branch active by running the command:

   gitxcheckoutx<branch-name>

Confirm you are on the correct branch by running:

   gitxstatus

which should print something similar to

   Onxbranchx<branch-name>

At this point you are all set to make the required modifications to the code in your branch. As you modify your branch, you can use git to save and track your changes. The following examples show how you can add a file to the list of local files that are being tracked, display differences between the current file and the original file, and commit the file.

   gitxaddxlibrary/LibUtilities/BasicUtils/my_new_file.cppxgitxdiffx--cached

Note --cached is necessary because my\_new\_file.cpp was staged using the gitxadd command above. Note, before you add (stage) the file, you can just use

gitxdiff.

To actually create the commit from the staged files, run:

   gitxcommitx-mx"Added␣X,␣Y,␣Z..."

This commits the file to your local repository. The first line of the log message should be a concise summary of the changes. You can use subsequent lines to provide more details if needed. Use gitxlog to see the list of previous commits and their messages. These changes are still local to your computer. To push them up to the main Nektar++ repository, use the command

   gitxpushxorigin

3.1.3 Building (CMake)

Nektar++ uses the CMake tool to manage the build process for the three supported operating systems: Linux, Windows, and OSX. For detailed instructions on how to use CMake to build Nektar++, including a list of required software dependencies and CMake option flags, please refer to the Nektar++ User Guide section 1.3.

3.1.4 Testing (CTest)

Before you are ready to have your code merged into the Nektar++ trunk, you should make sure that it passes the built-in test suite - in addition to any new tests that you have added for your updates. To run the test suite, on the command line type:

    ctestx[-j#]

The -j# optional argument will run # tests in parallel taking advantage of multiple cores on your machine. It is highly recommended that your use all available cores to minimize the amount of time spent waiting for the tests to complete. There are currently several hundred built-in unit tests for Nektar++.

For more information on testing, see “Software Testing Approaches” below.

3.1.5 Merge Requests (Gitlab)

The final step in contributing your code to the Nektar master repository is to submit a merge request to the development team using the Nektar++ gitlab website:

https://gitlab.nektar.info/nektar/nektar/merge\_requests

Submitting a merge request will automatically trigger the continuous integration (CI) system, which will build and test your code on a range of platforms. You can monitor the progress of these tests from the merge request page. Selecting individual workers will take you to the buildbot website from which you can examine any failures which have occurred.

Documentation and Tutorials

[next] [prev] [prev-tail] [tail] [up] 3.2 Documentation and Tutorials Documentation for Nektar++ is provided in a number of forms:

User Guide (LaTeX, compiled to pdf or html)

Source code documentation (Doxygen compiled to html)

3.2.1 Dependencies To build the User Guide and Developer’s Guide, the following dependencies are required:

texlive-base

texlive-latex-extra texlive-science texlive-fonts-recommended texlive-pstricks

imagemagick 3.2.2 Compiling the User Guide To compile the User Guide:

Configure the Nektar++ build tree as normal. Run

   makexuser-guide-pdf

to make the PDF version, or run

   makexuser-guide-html

to make the HTML version. 3.2.3 Developers Guide To compile the Developer’s Guide:

Configure the Nektar++ build tree as normal. Run

         makexdevelopers-guide-pdf

to make the PDF version, or run

         makexdevelopers-guide-html

to make the HTML version.

3.2.4 Compiling the code documentation To build the Doxygen documentation, the following dependencies are required:

doxygen

graphviz

To compile the code documentation enable the NEKTAR_BUILD_DOC option in the ccmake configuration tool. You can then compile the HTML code documentation using:

    makexdoc

Compiling Tutorials

[next] [prev] [prev-tail] [tail] [up] 3.3 Compiling Tutorials If you are using a clone of the Nektar++ git repository, you can also download the source for the Nektar++ tutorials which is available as a git submodule.

From a Nektar++ working directory (e.g. $NEKPP):

         gitxsubmodulexinitxgitxsubmodulexupdatex--remote

From your build directory (e.g. $NEKPP/build), re-run cmake to update the build system to include the tutorials

         cmakex../

Compile each required tutorial, for example

         makexflow-stability-channel

Core Nektar++ Programming Concepts

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

3.4 Core Nektar++ Programming Concepts

This section highlights some of the programming features that are used extensively within Nektar++. While much of the code consists of standard C++ practices, in some of the core infrastructure there are several practices that may be only familiar to programmers who have developed code using more advanced C++ features. Below we give a short summary of these entities in order to provide a starting point when working with these features. We begin with more well known features and end with some advanced techniques. Note, it is not the purpose of the following sections to cover in detail each of these important concepts, but instead to give a brief overview of them such that the developer may look to other, more in-depth, sources if they require further guidance.

3.4.1 Namespaces

Many C++ software projects place their code in a namespace so as to avoid conflicts with other code when included in larger applications. It is important to note that Nektar++ uses a hierarchy of namespaces for most of the defined data structures. The top level namespace is always “Nektar”, with the second level usually corresponding to the name of the library to which the code belongs. For example:

      1 namespacexNektarx 2 {x 3 namespacexStdRegionsx 4 {x 5 xxxx...x 6 }x 7 }

With this in mind, when you see something like Nektar::SpatialDomains::..., you can usually assume that the second item (in this case SpatialDomains) is a namespace, and not a class.

Note: To make better use of the 80 character width, generally enforced across the Nektar++ source code, we choose not to indent the contents of namespace blocks.

3.4.2 C++ Standard Template Library (STL)

Nektar++ uses of the C++ STL extensively. This consists of common data structures and algorithms, such as map and vector, as well as many of the extensions once found in the Boost library that have become part of the C++ standard and are now used directly.

One of the most important of these features is the use of Shared Pointers (std::shared_ptr). Most developers are somewhat familiar with “smart pointers” (pointers used to track memory allocation and to automatically deallocate the memory when it is no longer being used) for data blocks that are shared by multiple objects. These smart pointers are used extensively in Nektar++ and one should be familiar with the dynamic_pointer_cast function and the concept of the weak_ptr. Dynamic casting allows for safely converting one type of variable into its base type (or vice versa). For example:

      1 std::shared_ptr<FilterCheckpoint>xsptrx=x 2 xxxxxxxxstd::dynamic_pointer_cast<FilterCheckpoint>(xm_filters[k]x);x 3 if(xsptrx!=xnullptrx)x 4 {x 5 xxx//xCastxsucceeded!x 6 }

The advantage of using the dynamic cast, in comparison to the C style cast, is that you can check the return value at run time to verify that the casting was valid. A weak_ptr is a pointer to shared data with the explicit contract that the weak pointer does not own the data (and thus will not be responsible for deallocating it). Weak pointers are used mostly for short-term access to shared data.

Another modern code utility used by Nektar++ to support shared pointers can be seen in Nektar++ classes which inherit from std::enable_shared_from_this. This allows a class member function to return a shared pointer to itself. Specifically, it makes available the function shared_from_this() which returns a shard pointer to the object in the given context.

While C++ shared pointers are a powerful resource, there are a number of intricacies that must be understood and followed when creating classes and using objects that will be managed by them. For those not familiar with the C++11 (or previously Boost) implementation, it is highly recommended that you study them in more detail then presented here.

3.4.3 typedefs

Like most other large codes, Nektar++ uses typedefs to create short names for new variable types. You will see examples of this throughout the code and taking a few minutes to look at the definitions will help make it easier to follow the code. In the following example, we create (and explicitly name) the type ExpansionSharedPtr to make the code that uses this type easier to follow. This is particularly true of nested STL data structures where repeated template declarations would make the code harder to follow. A couple of examples are shown below:

      1 typedefxstd::shared_ptr<Expansion>xExpansionSharedPtr;x 2 typedefxstd::shared_ptr<std::vector<std::pair<GeometrySharedPtr,xint>>>x 3 GeometryLinkSharedPtr;

If you are not familiar with the use of typedefs, you should take time to read about them (there are many short summaries available on the web).

3.4.4 Forward Declarations

There are two ways that an existing class type can be specified when declaring a new class in a header file. The existing class can either be declared in name only, or declared in its entirety, before being used. In the latter case, one typically includes the header file declaring the full class. If the new class declaration only references the existing class in the form of a pointer or reference then the entire class declaration is not needed and the compiler only needs an assurance that the class exists. For this case, we can use a forward declaration which tells the compiler the name of the existing class. However, if functions of the existing class are called (within the new header file) or the class is used by value, then the full declaration is needed.

Forward declaring a class is achieved as shown in the following example:

      1 classxLinearSystem;

This statement tells the compiler the class LinearSystem exists and, as long as we only make reference to it as a pointer (LinearSystem*xl) or by reference (constxLinearSystemx&xl), then the compiler does not require any further information.

An advantage to using forward declarations where possible is that the header file does not need to #include the entire existing class and any header files referenced within. This allows for a cleaner header files and faster compilation as the compiler can process (often significantly) fewer lines of code.

Note: The full class declaration is most likely needed in the new class implementation file (.cpp) as reference to the existing class’s members will presumably be made.

3.4.5 Templated Classes and Specialization

Most C++ developers are familiar with basic class templating. However, many have not needed to use explicit template specialization. This is the process of implementing customised behaviour for one or more of the specific instances of a template when the compiler will not be able to instantiate a generic version for the class, or when different code is needed based on different versions of the class. For example:

      1 template<typenamexDim,xtypenamexDataType>x 2 classxArray;x 3 x 4 template<typenamexDataType>x 5 classxArray<OneD,xconstxDataType>x{x 6 xx//xExplicitxcodingxofxclassxmethodsxandxvariablesxspecificxtox 7 xx//xthisxversionxofxArrayxarexfoundxhere.x 8 };

In the above example, on the first line the generic templated Array class is declared. There are two template parameters: the dimension and the element type. The second line shows an explicit template specialisation of the Array class for a one-dimensional (version of) Array. When explicitly specialising a class, the programmer will write code that is specific to the datatype used in specifying the class. This includes explicitly writing code for one, some, or all of the methods of the class.

It is important to understand template specialization when dealing with the Nektar++ core libraries so that the developer can determine which (specialized version of the) class is being used, and to know that when updating classes with varied specializations, that it may be required to update code in several places (ie, for each of the specializations).

3.4.6 Multiple Inheritance and the virtual Keyword

When diving into many Nektar++ classes, you will see the use of multiple inheritance (where a class inherits from more than one parent class). When the parent class does not inherit from other classes, then the inheritance is straightforward and should not cause any confusion. However, when a class has grandparents, many times that grandparent class is the same class but is inherited through multiple parents. To account for this, class inheritance should use the virtual keyword. This specifies that if a class has multiple grandparents (that happen to be the same class), that only one copy of the grandparent class members should actually be instantiated. For example:

      1 classxExpansion2Dx:xvirtualxpublicxExpansion,x 2 xxxxxxxxxxxxxxxxxxxxvirtualxpublicxStdRegions::StdExpansion2Dx 3 {...}

3.4.7 Virtual Functions and Inheritance

Within Nektar++, classes that inherit from a parent class and override one of the parent class methods, use the concept of virtual functions. The function is prefixed with a v_, such as v_Function(), as a naming convention. This is a visual reminder that the function overrides a parent class function. For example:

      1 NekDoublexTriExp::v_Integral(constxArray<OneD,xconstxNekDouble>x&inarray)}

3.4.8 Const keyword

While the const keyword is known to most C++ developers, it is used (as it should be) liberally in Nektar++ for functions, function parameters, returning pointers to class data, and variable constants within functions. It is easy to neglect using const to mark all cases where a variable should be considered constant. However, its use can substantially reduce accidental errors and allow for accelerated debugging. The const qualifier should be used wherever a variable does not change including 1) parameters passed to functions, 2) variables in functions (or classes) that do not change value during their lifetime, 3) on the return type of functions that return pointers to data that should not be changed, and

4) on methods that do not change data within the class. The compiler will then produce an error if we (accidentally) attempt to make a change which violates a

const.

3.4.9 Function pointers and bind

Function pointers (std::function) are similar to pointers to data, except that they point to functions - and thus allow a function to be invoked indirectly (in other words, without explicitly writing the function call (name) directly in code). This technique is used by Nektar++ in a number of places, with NekManager being a prime example. The NekManager class is used to create objects of a specific type during the execution of the program. When a NekManager is created (constructed), it is provided with a pointer to a function that will (later) be called to generate the objects to be managed when required. While the creation function that is provided to the NekManager takes a number of parameters, in many cases some of the values to those parameters will be fixed. To handle this situation, Nektar++ uses the std::bind(xfx) function, which creates a new function based on supplied original function f, but specifies that one or more parameters of f are fixed at the time that f is created and only those bound parameter values will be used when f is later invoked.

3.4.10 Memory Pools and NekArray

An Array is a thin wrapper around native arrays. Arrays provide all the functionality of native arrays, with the additional benefits of automatic use of the Nektar++ memory pool, automatic memory allocation and deallocation, bounds checking in debug mode, and easier to use multi-dimensional arrays.

Arrays are templated to allow compile-time customization of its dimensionality and data type.

Parameters:

Dim Must be a type with a static unsigned integer called Value that specifies the array’s dimensionality. For example

           1 structxTenDx{x 2 xxxxstaticxunsignedxintxValuex=x10;x 3 };

DataType The type of data to store in the array.

It is often useful to create a class member Array that is shared with users of the object without letting the users modify the array. To allow this behavior, Array<Dim,xDataType> inherits from Array<Dim,xconstxDataType>. The following example shows what is possible using this approach:

      1 classxSamplex{x 2 public:x 3 xxxxArray<OneD,xconstxdouble>&xgetData()xconstx{xreturnxm_data;x}x 4 xxxxvoidxgetData(Array<OneD,xconstxdouble>&xout)xconstx{xoutx=xm_data;x}x 5 x 6 private:x 7 xxxxArray<OneD,xdouble>xm_data;x 8 };

In this example, each instance of Sample contains an array. The getData method gives the user access to the array values, but does not allow modification of those values.

Design Patterns

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

3.5 Design Patterns

3.5.1 Template pattern

The template pattern is used frequently within Nektar++ to provide a common interface to a range of related classes. The base class declares common functionality or algorithms, in the form of public functions, deferring to protected virtual functions where a specific implementation is required for each derived class. This ensures code external to the class hierarchy sees a common interface.

In the top level parent class you will find the interface functions, such as Function(), declared as public members. In some cases, these may implement a generic algorithm common to all classes. In the limit that the function is entirely dependent on the derived class, it may call through to a virtual counterpart, usually named v_Function(). These functions are usually protected to allow them to be called directly by other classes in the inheritance hierarchy, without exposing them to external classes.

As an example of this, let us consider a triangle element (TriExp). The TriExp class (eventually) inherits from the StdExpansion class. The StdExpansion defines the Integral() function which is used to provide integration over an element. However, in this case, the implementation is

shape-specific. Therefore StdExpansion::Integral() calls the (in this case) TriExp::v_Integral()

function. We should also note that while TriExp::v_Integral() does setup work, it then makes use of its parent’s StdExpansion2D::v_Integral() function to calculate the final value. This is only possible if the v_Integral() function was declared as protected.

3.5.2 Abstract Factory Pattern

Nektar++ makes extensive use of the Factory Pattern. Factories are used to create (allocate) instances of classes using class-specific creator functions. More specifically, a factory will create a new object of some sub-class type but return a base class pointer to the new object. In general, there are two ways that a factory knows what specific type of object to generate: 1) The Factory’s build function (CreateInstance()) is passed a key that details what to build; or 2) The factory may have some intrinsic knowledge detailing what objects to create. The first case is almost exclusively used throughout Nektar++. The factory pattern provides the following benefits:

Encourages modularisation of code such that conceptually related algorithms are grouped together;

Structuring of code such that different implementations of the same concept are encapsulated and share a common interface;

Users of a factory-instantiated modules need only be concerned with the interface and not the details of underlying implementations;

Simplifies debugging since code relating to a specific implementation resides in a single class;

The code is naturally decoupled to reduce header-file dependencies and improves compile times;

Enables implementations (e.g. relating to third-party libraries) to be disabled through the build process (CMake) by not compiling a specific implementation, rather than scattering preprocessing statements throughout the code.

Implementing the factory pattern The NekFactory class implements the factory pattern in Nektar++. There are two distinct aspects to creating a factory-instantiated collection of classes: defining the public interface, and; registering specific implementations. Both of these tasks involve adding mostly standard boilerplate code.

It is assumed that we are writing a code which implements a particular concept or capability within the code, for which there are (potentially) multiple implementations. The reasons for multiple implementations may be low level, such as alternative algorithms for solving a linear system, or high level, such as selecting from a range of PDEs to solve. The NekFactory can be used in both cases and applied in exactly the same way.

Creating a concept (base class)

A base class must be defined which prescribes an implementation-independent interface. In Nektar++, the template method pattern (see Sectionx3.5.1 above) is used, requiring public interface functions to be defined which call through to protected virtual implementation methods. This is because the factory returns the newly created object via a base-class pointer and the objects will almost always be used via this base class pointer. Without a public interface in the base class, much of the benefits and generalisation of code offered by the factory pattern would be lost. The virtual functions will be overridden in the specific implementation classes. In the base class these virtual methods should normally be defined as pure virtual, since there is typically no implementation and we will never be explicitly instantiating this base class.

As an example we will create a factory for instantiating different implementations of some

concept MyConcept, defined in MyConcept.h and MyConcept.cpp.

First in MyConcept.h, we need to include the NekFactory header

      1 #includex<LibUtilities/BasicUtils/NekFactory.hpp>

The following code should then be included just before the base class declaration (within the same namespace as the class):

      1 classxMyConceptx 2 x 3 //xDatatypexforxthexMyConceptxfactoryx 4 typedefxLibUtilities::NekFactory<xstd::string,xMyConcept,x 5 xxxxxxxxxxxxParamType1,x 6 xxxxxxxxxxxxParamType2x>xMyConceptFactory;x 7 MyConceptFactory&xGetMyConceptFactory();

The template parameters to the NekFactory define the datatype of the key used to retrieve a particular implementation (usually a string, enum or custom class such as MyConceptKey), the base class datatype (in our case MyConcept and a list of zero or more parameters which are taken by the constructors of all implementations of the type MyConcept (in our case we have two). Note that all implementations must take the same parameter list in their constructors. Since we have not yet declared the base class type MyConcept, we have forward-declared it above the NekFactory type definition.

The normal definition of our base class then follows:

      1 classxMyConceptx 2 {x 3 xxxxpublic:x 4 xxxxxxxxMyConcept(ParamType1xp1,xParamType2xp2);x 5 xxxxxxxx...x 6 };

We must also define a shared pointer for our base class, which should be declared outside the base class declaration.

      1 typedefxboost::shared_ptr<MyConcept>xMyConceptShPtr;

Creating a specific implementation (derived class)

A new class, derived from the base class above, is defined for each specific implementations of a concept. It is these specific implementations which are instantiated by the factory.

In our example we will have an implementations called MyConceptImpl1 defined in MyConceptImpl1.h and MyConceptImpl1.cpp. In the header file we include the base class header file

      1 #includex<MyConcept.h>

We then define the derived class as one would normally:

      1 classxMyConceptImpl1x:xpublicxMyConceptx 2 {x 3 ...x 4 };

In order for the factory to work, it must know two things:

     that MyConceptImpl1 exists; and

how to create an instance of it.

To allow the factory to create instances of our class we define a creator function in our class, which may have arbitrary name, but is usually called create out of convention:

      1 ///xCreatesxanxinstancexofxthisxclassx 2 staticxMyConceptSharedPtrxcreate(x 3 xxxxxxxxxxxxParamType1xp1,x 4 xxxxxxxxxxxxParamType2xp2)x 5 {x 6 xxxxreturnxMemoryManager<MyConceptImpl1>::AllocateSharedPtr(p1,xp2);x 7 }

The example above the create function simply creates an instance of MyConceptImpl1 using the Nektar++ memory manager and the supplied parameters. It must be a static function because we are not operating on any existing instance and it should return a shared pointer to a base class object (rather than a MyConceptImpl1 shared pointer), since the point of the factory is that the calling code does not know about specific implementations. An advantage of having each class providing a creator function is that it allows for two-stage initialisation – for example, initialising base-class variables based on the derived type.

The final task is to register each of our implementations with the factory. This is done using the RegisterCreatorFunction member function of the NekFactory. However, we wish this to happen as early on as possible (so we can use the factory straight away) and without needing to explicitly call the function for every implementation at the beginning of our program (since this would again defeat the point of a factory)! One solution is to use the function to initialise a static variable: it will force the function to be executed prior to the start of the main() routine, and can be located within the class it is registering, satisfying our code decoupling requirements.

In MyConceptImpl1.h we define a static class member variable with the same datatype as the key used in our factory (in our case std::string)

      1 staticxstd::stringxclassName;

The above variable can be private since it is typically never actually used within the code. We

then initialise it in MyConceptImpl1.cpp
      1 stringxMyConceptImpl1::classNamex 2 xxxxxxxx=xGetMyConceptFactory().RegisterCreatorFunction(x 3 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"Impl1",x 4 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxMyConceptImpl1::create,x 5 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"Firstximplementationxofxmyxconcept.");

The first parameter specifies the value of the key which should be used to select this implementation. The second parameter is a function pointer to the static function which should be used by the factory to instantiate our class. The third parameter provides a description which can be printed when listing the available MyConcept implementations. A specific implementation can be registered with the factory multiple times if there are multiple keys which should instantiate an object of this class.

Instantiating classes

To create instances of MyConcept implementations elsewhere in the code, we must first include the ”base class” header file

      1 #includex<MyConcept.h>

Note we do not include the header files for the specific MyConcept implementations anywhere in the code (apart from MyConceptImpl1.cpp). If we modify the implementation, only the implementation itself requires recompilation and the executable relinking.

We create an instance by retrieving the MyConceptFactory and calling the CreateInstance member function of the factory, for example,

      1 ParamTypexp1x=x...;x 2 ParamTypexp2x=x...;x 3 MyConceptShPtrxpx=xGetMyConceptFactory().CreateInstance(x"Impl1",xp1,xp2x);

Note that the instance of the specific implementation is used through the pointer p, which is of type MyConceptShPtr, allowing the use of any of the public interface functions in the base class (and therefore the specific implementations behind them) to be called, but not directly any functions declared solely in a specific implementation.

Software Testing Approaches

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

3.6 Software Testing Approaches

3.6.1 Unit Tests

Unit testing, sometimes called “module testing” or “element testing”, is a software testing method by which individual “units” of source code are tested to determine whether they are fit for use [42]. Unit tests are added to Nektar++ through the CMake system, and implemented using the Boost test framework. As an example, the set of linear algebra unit tests is listed in this file:

.../library/UnitTests/LibUtilities/LinearAlgebra/CMakeLists.txt

and the actual tests are implemented in this file:

.../library/UnitTests/LibUtilities/LinearAlgebra/TestBandedMatrixOperations.cpp

To register a new test, you use BOOST_AUTO_TEST_CASE(xTestNamex), implement the unit test, and

test the result using BOOST_CHECK_CLOSE(...), BOOST_CHECK_EQUAL(...), etc. Unit tests are

invaluable in maintaining the integrity of the code base and for localizing, finding, and debugging errors entered into the code. It is important to remember a unit test should test very specific functionality of the code - in the best case, a single function should be tested per unit test.

While it is beyond the scope of this document to go into more detail on writing unit tests, a good summary of the Boost test system can be found here:

http://www.boost.org/doc/libs/1\_63\_0/libs/test/doc/html/ .

3.6.2 Integration, System and Regression Tests

Integration testing involves testing ecosystems of components and their interoperability. System testing tests complete applications and regression testing focuses on ensuring previously fixed bugs do not resurface. In Nektar++ all of these are often colloquially referred to as regression testing. It is not white-box in that it does not examine how the code arrives at a particular answer, but rather in a black-box fashion tests to see if code when operating on certain data yields the predicted response [42].

3.6.3 Continuous Integration

Nektar++ uses the buildbot continuous integration to perform testing of the code across multiple operating systems. Builds are automatically instigated when merge requests are opened and subsequently when the associated branches receive additional commits.

For more information, go to:

http://buildbot.nektar.info [next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] I Building-Blocks of Our Framework (Inside the Library)

[next] [prev] [prev-tail] [tail] [up] PartxIBuilding-Blocks of Our Framework (Inside the Library)

4 Inside the Library: LibUtilities

[next] [prev] [prev-tail] [tail] [up] Chapterx4Inside the Library: LibUtilities In this chapter, we walk the reader through the different components of the LibUtilities Directory. We have ordered them in alphabetical order by directory name, not by level of importance or relevance to the code. Since all of these items are considered foundational to Nektar++, they should all be considered equally important and relevant. Along the same lines – since all of these areas of the code represent the deepest members of the code hierarchy, these items should rarely be modified.

xx3.6.3 Continuous Integration

BasicConst

[next] [prev] [prev-tail] [tail] [up] 4.1 BasicConst This directory contains two important files for all of Nektar++: NektarUnivConsts.hpp and NektarUnivTypeDefs.hpp. The file NektarUnivConsts.hpp contains various default constants used within Nektar++ as seen here:

       x 1 xxxxxxxxstaticxconstxNekDoublexkNekMinResidInitx=x1e16;x 2 xxxxxxxxstaticxconstxNekDoublexkVertexTheSameDoublexx=x1.0e-8;x 3 xxxxxxxxstaticxconstxNekDoublexkGeomFactorsTolx=x1.0e-8;x 4 xxxxxxxxstaticxconstxNekDoublexkNekZeroTolx=x1.0e-12;x 5 xxxxxxxxstaticxconstxNekDoublexkGeomRightAngleTolx=x1e-14;x 6 xxxxxxxxstaticxconstxNekDoublexkNekSqrtTolx=x1.0e-16;x 7 xxxxxxxxstaticxconstxNekDoublexkNekIterativeTolx=x1e-09;x 8 xxxxxxxxstaticxconstxNekDoublexkNekSparseNonZeroTolx=x1e-16;x 9 x 10 xxxxxxxx//xTolerancesxforxmeshxgenerationxandxCADxhandlingx 11 xxxxxxxxstaticxconstxNekDoublexGeomTolx=x1E-2;x 12 xxxxxxxxstaticxconstxNekDoublexCoinTolx=x1E-6;x 13 xxxx}

The file NektarUnivTypeDefs.hpp contains the low level typedefs such as: NekDouble, NekInt, OneD, TwoD, ThreeD, FourD, and enumerations such as Direction (xDir, yDir and zDir) and OutputFormat.

BasicUtils

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

4.2 BasicUtils

This directory contains some of the lowest level basic computer science routines within Nektar++. The directory currently contains the following:

11 ArrayEqualityComparison.cpp 11  ArrayPolicies.hpp  11 CompressData.cpp
11  CompressData.h  11 ConsistentObjectAccess.hpp 11  CsvIO.cpp
11  CsvIO.h  11  Equation.cpp  11  Equation.h
11  ErrorUtil.hpp  11  FieldIO.cpp  11  FieldIO.h
11  FieldIOHdf5.cpp  11  FieldIOHdf5.h  11 FieldIOXml.cpp
11  FieldIOXml.h  11  FileSystem.cpp  11 FileSystem.h
11  H5.cpp  11  H5.h  11 HashUtils.hpp
11  Metis.hpp  11  NekFactory.hpp  11 NekManager.hpp
11  OperatorGenerators.hpp  11  ParseUtils.cpp  11  ParseUtils.h
11  Progressbar.hpp  11  PtsField.cpp  11  PtsField.h
11  PtsIO.cpp  11  PtsIO.h  11 RawType.hpp
11  Scotch.hpp  11  SessionReader.cpp  11 SessionReader.h
11  ShapeType.hpp  11  SharedArray.hpp  11  Tau.hpp
11  Thread.cpp  11  Thread.h  11 ThreadBoost.cpp
11  ThreadBoost.h  11  Timer.cpp  11  Timer.h
11  VDmath.hpp  11  VDmathArray.hpp  11  Vmath.cpp
11  Vmath.hpp  11  VmathArray.hpp  11  VtkUtil.hpp

11

We have used bold to denote (as examples) routines at our used throughout Nektar++. They are in this sense “fundamental”. Note that this list includes input/output routines (e.g. CompressData, FieldIo, H5), partitioning (e.g. Metis and Scotch) and Threading (e.g. Thread and ThreadBoost).

Communication

[next] [prev] [prev-tail] [tail] [up] 4.3 Communication This directory contains files related to our distributed memory communication model. In particular, this directory contains files that help encapsulate MPI (Message Passing Interface) routines, as well as the Gather-Scatter (GS) and Xxt routines of Paul Fisher (Argonne National Lab and UIUC).

11  Comm.cpp  11 CommMpi.cpp 11  GsLib.hpp
11  Comm.h  11 CommMpi.h 11 Transposition.cpp
11 CommDataType.cpp 11 CommSerial.cpp 11 Transposition.h
11 CommDataType.h 11 CommSerial.h 11  Xxt.hpp

11

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

FFT

[next] [prev] [prev-tail] [tail] [up] 4.4 FFT This directory contains files related to the use of Fast Fourier Transforms within Nektar++.

The two groups of files are as follows:

NekFFTW.h/NekFFTW.cpp : Wrapper around FFTW library; and

NektarFFT.h/NektarFFT.cpp : Fast Fourier Transform base class in Nektar++.

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

Foundations

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

4.5 Foundations

The two basic building blocks of all that is done in Nektar++ are the concepts of Points and of a Basis. The Point objects denote positions in space, either on compact domains (normally [−1, 1]^{d} where d is the dimension in a reference domain mapped to world-space) or periodic domains such as (0, 2π] (i.e. in the case of points used in Fourier expansions). The Basis objects denote functions (e.g. polynomials) evaluated at a given set of points.

We rely heavily on the concept of “managers”, and this concept has its start right here in Foundations. When we started our development of Nektar++, we realized that a pure encapsulation strategy for Expansions would consist of every expansion holding a pointer to the basis (of some sort) that was used to specify it, and then correspondingly the basis pointing to a collection of points at which it was evaluated. From the computer science perspective, this is all very reasonable except for the fact that there are many cases in spectral/hp element expansions in which multiple elements use the same reference space basis functions, and correspondingly those basis functions are often evaluated at the same set of points (quadrature points). This led us to the conclusion that if we were going to have a code that could run over large number of elements (at the time tens of millions) without excessive memory usage, we would have to try to avoid as much of this duplication as possible. This led us to introduce to Nektar++ two important computer science tools: smart pointers and managers (i.e. a factory pattern).

Smart pointers were originally introduced as an add-on from the Boost library and later directly incorporated into Nektar++ when they were natively handled by C++11.

With the tradition view of pointers, one has a variable that holds an address that points to a place in memory. Typically, in type-safe programming, the way this memory is to be interpreted is denoted by the type of pointer: a double pointer (i.e. double*) is a pointer to memory that should be interpreted as a double. The downside of this traditional view occurs when multiple pointers all point to the same place in memory. The question then arises: which pointer (or more properly which thread of control containing a pointer) is responsible for deleting the memory when it is no longer needed? Again, in the traditional view that each pointer points to unique memory and hence when the memory is not longer needed the memory can be released no longer is viable when multiple pointers all point to the same piece of memory. This probably was solved with the invention of the smart pointer, an object that is wrapped up to look like a pointer, and yet contains a reference counter within it. When a smart pointer is initiated and memory is originally assigned to it, the reference counter is incremented. As more and more references to the piece of memory are created, the increment counter continues to increase. As pointers go out of scope and are longer valid, the increment is decremented. Only when all the references to a particular piece of memory are removed can a free (i.e. delete) operation be accomplished. The introduction of smart pointers into C++ and into Nektar++ were particularly advantageous to us as they allowed us to have lots of different elements, expansions, etc. point to the same fundamental data structure in memory (for instance, a Basis object evaluated at a particular set of points) without having to worry about improper deallocation or memory leakage due to pointer mishandling.

The second important feature we exploited within the redesign of Nektar++ was the idea of managers (i.e. a factory pattern). You will see this use of the term manager throughout Nektar++. For us, a manager has the following characteristics: it expects as input a key, which gives the manager sufficient information to know the object for which it is receiving a request. If the manager already has one of those objects in storage, it passes back a (smart) pointer that that object. The manager itself continues to hold on to a pointer for the object (in case other requests are made of it for that object). If the manager does not have a copy of the requested object, then it has the ability to generate the object “on the fly” and store it for future use. This last component of the manager is the factory feature that we mentioned earlier – a manager has registered with it a collection of methods that know how to generate the various objects that might be requested of it.

Based upon these two concepts: two fundamental managers that exist at this level of the library are the PointsManager and the BasisManager. These are contained within ManagerAccess.h/cpp, and are defined as singletons and defined globally so that they can be accessed throughout Nektar++. Within these files you will find the RegisterCreator method calls that link particular Create methods with various keys. Only in rare cases will these every need to be modified, and only infrequently possibly added to (if additional Points or Basis information is given).

4.5.1 Points

One of the two primary data structures in this directory is Points. A points object consists of a PointsKey and then extra information needed to facilitate various operations at points. The basic layout of the data structure is shown in Figure 4.1. pict Figurex4.1: The basic Points class data object. It consists of a PointsKey (used by a PointsManager) and various other data members needed for point operations.

The PointsKey object has three basis data members: the number of points, an enum specifying the point distribution, and a optional scaling factor. In almost all of Nektar++, only the first two are needed to uniquely specify a collection of points. The latter variable was added as part of the extensions of Nektar++ to accommodate NekMesh. As stated above, if you pass a PointsKey to a PointsManager you will receive back a pointer to a Points object with its various data members populated, among them a PointsKey identifier for the object. As can be seen by the diagram, the Points object has a collection (array) of three pointers to point positions. These three array pointers allow Points objects to point to 1D, 2D or 3D point positions. Corresponding to the dimension and associated with each point position is a weight associated with integration. This was done to facilitate both non-tensor product and tensor product constructions. To demonstrate non-tensor product constructions, consider if you were to be dealing with a 2D point set that is natively 2D (that is, not based upon tensor product construction such as the Electrostatic Points), then the number of points in the PointsKey would be the size of each array associated with m_points, and only one of the m_weight arrays (associated with m_weight[0]) would be valid such at:

.  ∫_{Ω_{2D}} f(x_{0}, x_{1})dx_{0}dx_{1} ≈ ∑_{i=0}^{numpoints−1} ω_{i}f(z_{0,i}, z_{1,i})

where z_{0,i} and z_{1,i} denote the points associated with m_points[0][i] and m_points[1][i] respectively. The weights in this case are calculated by routines found in NodalUtils via forming a Vandermonde system, solving for Lagrange basis functions that pass through a particular set of points, and then integrating those basis functions to obtain weights. This operation is done sufficiently often that we will try to spell it out after the Fekete Point paragraph below.

If one is dealing with a tensor product constructed 2D nodal point set, one would have:

.  ∫_{Ω_{2D}} f(x_{0}, x_{1})dx_{0}dx_{1} ≈ ∑_{i=0}^{numpoints−1} [ ∑_{j=0}^{numpoints−1} ω_{i}ω_{j}f(z_{0,i}, z_{1,j})]

where z_{0,i} and z_{1,i} denote the points associated with m_points[0][i] and m_points[1][i] respectively and ω_{i} and ω_{j} are the weights associated with the two coordinate directions respectively. The tensor product construction above was implemented in Nektar++ assuming we would create derived nodal point sets for quadrilaterals and hexahedron; however, throughout much of Nektar++, we use the tensor product forms explicitly by specifying 1D point sets in the different directions as needed. This particular point will become more apparent when you get to reading the Chapter on StdRegions (Chapter 5). You will encounter, for instance, that a Standard Quadrilateral Expansion (StdQuadExp) requires two 1D Point objects denoting the two coordinates for integration. In this case, the number of points and number of weights match up per direction.

Electrostatic Points:

The Electrostatic (Nodal) Points on the triangle and tetrahedron are based upon the work of Hesthaven [38]. They are an attempt to find the minimizer of a potential energy function based upon an electrostatic point source analogy, and provide points that correspondingly have low Lebesgue constant. The point positions as given in the paper are stored in the files NodalTriElecData.h and NodalTetElecData.h. Hesthaven uses a condensed format for the nodal positions which rely on the symmetries of the points (if you count the points in the file, you will see that fewer points are given than are needed to span the polynomial space. For example, for degree one, a single point it given. This point along with its three-fold symmetries represent all three points needed to support the linear space). The routine CalculatePoints() expands the condensed format of the points to the full set based upon the one, three (a and b denoting two types of three-symmetries) and six symmetries. Note that the ordering of the nodes expanded from the file does not respect any particular geometric ordering, so we reorder the points using NodalPointReorder2d(), which reorganizes the points in vertex points, edge points, and then interior points (i.e. a geometric decomposition that aligns more naturally with the way we organize nodes/modes within Nektar++). In Figure 4.2 we show an example of this reordering done for the points needed for expressing a total degree three expansion over a triangle.

Figurex4.2: Example of how nodes are reordered to meet a geometric ordering. On the left, we show the nodal positions that would be used for building a total degree three polynomial space over a triangle ordered in canonical form. On the right, we show the reordering of those nodes to follow a vertex, edge and interior ordering. Note that in the case of a tetrahedron, the ordering would be vertex, edge, face and then interior.

Fekete Points:

The Fekete (Nodal) Points on a triangle are based upon the work of Taylor and Wingate [56,x57], and are an alternative nodal point distribution on both triangles and tetrahedra. As an alternative strategy to the one mentioned above, Taylor and Wingate explicitly attempt to find a point distribution that minimizes the Lebesgue constant. For their particular point set, the minimization if over the Lebesgue function itself (not the electrostatic potential function, which can be viewed as a proxy for this function). The point positions as given in the paper are stored in the file NodalTriFeketeData.h and NodalTetFeketeData.h. Taylor and Wingate use a condensed format for the nodal positions similar to that used by Hesthaven.

Building Lagrange Interpolating Polynomials:

One common operation done in both the case of the Electrostatic Points and the Fekete Points is to build the Lagrange interpolating polynomial associated with a point set. This is needed, for instance, for computing the integration weights associated with a particular collection of points (as mentioned above). This operation is done sufficiently often that we will try to spell it out here. Forming the Lagrange basis functions for a particular set of points, in particular in 2D and 3D, is non-trivial. The most common way of arriving at a set of actionable basis functions is to use a known basis (for which we are able to explicitly write down its form) and find the linear combination of those basis functions that yield the Lagrange basis. For example, let us assume we are dealing with the set of Electrostatic Points on a triangle. The number of points, numpoints, will specify the number of “modes” (coefficients) we expect to find. The total degree of the space is given by N = \frac{(M+1)(M+2)}{2} where N is the degree of the polynomials and M is the number of points. As an example: note that when M = 3, we get N = 1 – that is, with three points we can support a polynomial of total degree one over a triangle. Let ϕ_{i} denote a known non-interpolating basis function. Although not done in practice, imagine that this is, for instance, a monomial. In 1D, this would equate to ϕ_{i}(x) = x^{i}, and in multiple dimensions this would equate to ϕ_{i}(x, y) = x^{i_{1}} y^{i_{2}} with some index mapping function i = σ(i_{1}, i_{2}) that gives us an index ordering. The total number of basis functions to span our degree N space is given by M. We first form the Vandermonde matrix V as follows:

. V = [ ϕ_{0}(x_{0,0}, x_{1,0}) ϕ_{1}(x_{0,0}, x_{1,0}) … ϕ_{N}(x_{0,0}, x_{1,0}) ; ϕ_{0}(x_{0,1}, x_{1,1}) ϕ_{1}(x_{0,1}, x_{1,1}) … ϕ_{N}(x_{0,1}, x_{1,1}) ; ⋮ ; ϕ_{0}(x_{0,N} , x_{1,N}) ϕ_{1}(x_{0,N} , x_{1,N}) … ϕ_{N}(x_{0,N} , x_{1,N}) ]

where each row i denotes the evaluation of our N basis functions at a given point (x_{0,i}, x_{1,i}). Since we have selected the number of basis functions to be what can be uniquely resolved by our N points, this matrix is square and invertible. The conditioning of the matrix is based upon our basis choice; hence we know that monomials are not a good choice beyond approximately cubics so in general we use our internal modified basis for this system. With V now available, we can form the linear system:

. V c = b

where c is the vector of coefficients and where b = b_{i} is a binary vector where the i^{th} entry is set to one when finding the coefficients for the i^{th} Lagrange basis function associated with a particular set of points. For each Lagrange basis function definition we seek to find, we update the right-hand-side vector and solve the linear system for the set of coefficients that denote the combination of our known basis that yields a Lagrange basis function. Note that since this operation (of “inverting" this linear system) must be done for each Lagrange basis function, we can optimize our operations by accomplishing LU decomposition on V first, which can be done in O(N^{2}) operations, and then each solution for the coefficient vector can be done in O(N) operations through backsolves.

Finding Integration Weights:

When the points are related to Gaussian quadrature, the array m_weights will contain the appropriate Gaussian quadrature weights computed using point and weight routines found in Foundations (our polymath library functions). In the case of weights associated with points sets that do not lie at the zeros of Jacobi polynomials, we must compute the weights directly based upon Lagrange interpolation through those points.

Although this can be done in general for any point distribution, we will focus here on evenly-spaced points (note that when one applies this approach to Chebyshev points, one arrives at Clenshaw-Curtis quadrature [58]). Within Nektar++, given evenly-spaced points on the interval [−1, 1] – if only one point is given, then the weight is set to 2.0 denoting a midpoint integration rule; otherwise, the weights are given by the following expression:

. w_{j} = ∫_{−1}^{1} ℓ_{j}(ξ)dξ

where ℓ_{j}(ξ) is the Lagrange basis function defined over an evenly-spaced set of points ξ_{k} ∈ [−1, 1]. For 2D and 2D nodal point sets (such as Electrostatic and Fekete Points), this same procedure is used. We first form the Lagrange interpolating functions via the algorithm given above (using the Vandermonde system), and then for each quadrature weight compute the integral of that Lagrange basis function using some known quadrature rule (such as tensor product quadrature via Duffy transformation; this point will be more easily understood after reading Chapter 5)

4.5.2 Basis

The second of the two primary data structures in this directory is Basis. A basis object consists of a BasisKey and the extra information needed to facilitate various operations on the basis. The basis layout of the data structure is shown in Figure 4.3. pict Figurex4.3: The basic Basis class data object. It consists of a BasisKey (used by a BasisManager) and various other data members needed for point operations.

The Basis object consists of a BasisKey (which in turn holds a PointsKey), a pointer to the positions at which the basis is evaluated and then arrays which hold the evaluation of the basis (m_bdata) and the evaluation of the derivates of the basis (m_dbdata). The storage layout for these data structures are shown in Figure 4.3. The number of rows is given by the number of modes associated with a particular expansion, and the number of columns is given by the number of points. These arrays are stored as contiguous memory to help avoid excessive memory hopping. pict Figurex4.4: Diagram of the basic memory layout of the arrays associated with the basis and its derivatives evaluated at points.

Interpreter

[next] [prev] [prev-tail] [tail] [up] 4.6 Interpreter This directory contains two files, AnalyticExpressionEvaluator.hpp and AnalyticExpressionEvaluator.cpp, both of which provide parser and evaluator functionality of analytic expressions.

Kernel

[next] [prev] [prev-tail] [tail] [up] 4.7 Kernel This directory contains two files, kernel.h and kernel.cpp, both of which are related to the line-SIAC (L-SIAC) filtering [26,x41] of FEM solutions. SIAC filtering takes the form:

. u^{∗}(x) = ∫_{−∞}^{∞} K_{H} ( \frac{ξ − x}{H} ) u_{h}(ξ)dξ

where u_{h} denotes the FEM (continuous Galkerkin or discontinuous Galerkin) solution and K_{H} represents a kernel function:

. K_{H}(x) = \frac{1}{H} ∑_{γ=0}^{M} c_{γ}ψ ( \frac{x}{H} − ζ_{γ}) .

In this expression, the functions ψ(x) are B-Splines. The two kernel files in this directory provide all the tools needed to construct the B-Spline functions and the kernel function described above.

Linear Algebra

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

4.8 Linear Algebra

This directory contains files related to linear algebra operations. Files such as NekMatrix, BlockMatrix and NekLinearSystem are fundamental to many of the operations we do within Nektar++. These files represent our attempt to encapsulate linear algebra operations in a way that make sense from the programmer perspective but for which we do not loose efficiency. As can be seen in these routines and throughout the code, we rely heavily on BLAS (Basic Linear Algebra Subroutines).

11  Arpack.hpp  11  Blas.hpp  11  BlasArray.hpp
11  BlockMatrix.cpp  11  BlockMatrix.hpp  11 CanGetRawPtr.hpp
11 ExplicitInstantiation.h 11  Lapack.hpp  11  MatrixBase.cpp
11  MatrixBase.hpp  11  MatrixFuncs.cpp  11  MatrixFuncs.h
11 MatrixOperations.cpp 11  MatrixOperations.hpp  11 MatrixStorageType.h
11  MatrixType.h  11 MatrixVectorMultiplication.cpp 11 NekLinAlgAlgorithms.hpp
11  NekLinSys.hpp  11  NekMatrix.hpp  11  NekMatrixFwd.hpp
11  NekPoint.hpp  11  NekTypeDefs.hpp  11  NekVector.cpp
11  NekVector.hpp  11  NekVectorFwd.hpp  11 NistSparseDescriptors.hpp
11  PointerWrapper.h  11  ScaledMatrix.cpp  11  ScaledMatrix.hpp
11 SparseDiagBlkMatrix.cpp 11 SparseDiagBlkMatrix.hpp 11  SparseMatrix.cpp
11  SparseMatrix.hpp  11  SparseMatrixFwd.hpp  11  SparseUtils.cpp
11  SparseUtils.hpp  11  StandardMatrix.cpp  11 StandardMatrix.hpp
11 StorageSmvBsr.cpp 11  StorageSmvBsr.hpp  11  TransF77.hpp
11  blas.cpp  11                             11

11

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

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

4.9 Memory

This directory contains three files, NekMemoryManager.hpp, ThreadSpecificPools.hpp and ThreadSpecificPools.cpp. The strategy used within Nektar++ was to preallocate a pool of arrays that could be used for various operations and then released back to the pool. This idea came about through profiling of the code early on – noticing that the new/delete operation of lots of small arrays used for temporary calculations was greatly slowing down the code. Like with our manager idea, we decided to invest in having a memory pool object what preallocated blocks of memory that could be requested and then returned back to the pool.

If Nektar++ is compiled with NEKTAR_MEMORY_POOL_ENABLED, the MemoryManager allocates from thread specific memory pools for small objects. Large objects are managed with the system supplied new/delete. These memory pools provide faster allocation and deallocation of small objects (particularly useful for shared pointers which allocate many 4 byte objects).

All memory allocated from the memory manager must be returned to the memory manager. Calling delete on memory allocated from the manager will likely cause undefined behavior. A particularly subtle violation of this rule occurs when giving memory allocated from the manager to a shared pointer.

Polylib

[next] [prev] [prev-tail] [tail] [up] 4.10 Polylib This directory contains polylib.h and polylib.cpp. These files contain foundational routines used for computing various operations related to Jacobi polynomials. The following abbreviations are used throughout the file:

z – Set of collocation/quadrature points w – Set of quadrature weights D – Derivative matrix h – Lagrange Interpolant I – Interpolation matrix g – Gauss k – Kronrod gr – Gauss-Radau gl – Gauss-Lobatto j – Jacobi m – point at minus 1 in Radau rules p – point at plus 1 in Radau rules

Points and Weights: The following routines are used to compute points and weights:

zwgj – Compute Gauss-Jacobi points and weights zwgrjm – Compute Gauss-Radau-Jacobi points and weights (z = −1) zwgrjp – Compute Gauss-Radau-Jacobi points and weights (z = 1) zwglj – Compute Gauss-Lobatto-Jacobi points and weights zwgk – Compute Gauss-Kronrod-Jacobi points and weights zwrk – Compute Radau-Kronrod points and weights zwlk – Compute Lobatto-Kronrod points and weights JacZeros – Compute Gauss-Jacobi points and weights

Derivative Matrices: The following routines are used to compute derivative matrices:

Dgj – Compute Gauss-Jacobi derivative matrix Dgrjm – Compute Gauss-Radau-Jacobi derivative matrix (z = −1) Dgrjp – Compute Gauss-Radau-Jacobi derivative matrix (z = 1) Dglj – Compute Gauss-Lobatto-Jacobi derivative matrix

Lagrange Interpolants: The following routines are used to compute Lagrange interpolation matrices:

hgj – Compute Gauss-Jacobi Lagrange interpolants hgrjm – Compute Gauss-Radau-Jacobi Lagrange interpolants (z = −1) hgrjp – Compute Gauss-Radau-Jacobi Lagrange interpolants (z = 1) hglj – Compute Gauss-Lobatto-Jacobi Lagrange interpolants

Interpolation Operators: The following routines are used to compute various interpolation operators:

Imgj – Compute interpolation operator gj->m Imgrjm – Compute interpolation operator grj->m (z = −1)

Imgrjp – Compute interpolation operator grj->m (z = 1) Imglj – Compute interpolation operator glj->m

Polynomial Evaluation: The following routines are used to evaluate Jacobi polynomials.

jacobfd – Returns value and derivative of Jacobi polynomial at point z jacobd – Returns derivative of Jacobi polynomial at point z (valid at z = −1, 1)

Time Integration

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

4.11 Time Integration

This directory consists of the following files:

11 TimeIntegrationScheme.h 11 TimeIntegrationScheme.cpp
11 TimeIntegrationWrapper.h 11 TimeIntegrationWrapper.cpp

11

The original time-stepping classes (found in TimeIntegrationScheme) where implemented by Vos et al. [62] and a detailed explanation of the mathematics and computer science concepts are provided there. After the original implementation, Cantwell et al. extended Vos’ original work by adding the time-stepping routines into a factory pattern (found in TimeIntegrationWrapper).

4.11.1 General Linear Methods

General linear methods (GLMs) can be considered as the generalization of a broad range of different numerical methods for ordinary differential equations. They were introduced by Butcher and they provide a unified formulation for traditional methods such as the Runge-Kutta methods and the linear multi-step methods such as Adams-Bashforth. From an implementation point of view, this means that all these numerical methods can be abstracted in a similar way. As this allows a high level of generality, it is chosen in Nektar++ to cast all time integration schemes in the framework of general linear methods.

For background information about general linear methods, please consult [14].

4.11.2 Introduction

The standard initial value problem can written in the form

\frac{dy}{dt} = f (t, y), y(t_{0}) = y_{0}

where y is a vector containing the variable (or an array of array containing the variables).

In the formulation of general linear methods, it is more convenient to consider the ODE in autonomous form, i.e.

\frac{dŷ}{dt} = f̂(ŷ), ŷ(t_{0}) = ŷ_{0}.

4.11.3 Formulation Suppose the governing differential equation is given in autonomous form, the n^{th} step of the GLM comprising

r steps (as in a multi-step method)

s stages (as in a Runge-Kutta method)

is formulated as:

\begin{aligned}Y_{i} & = Δt ∑ _{j=0} ^{s−1} a_{ij} F_{j} + ∑ _{j=0} ^{r−1} u_{ij} ŷ_{j}^{[n−1]} ,  i = 0, 1, x, s − 1 \\ ŷ_{i}^{[n]} & = Δt ∑ _{j=0} ^{s−1} b_{ij} F_{j} + ∑ _{j=0} ^{r−1} v_{ij} ŷ_{j}^{[n−1]} ,  i = 0, 1, x, r − 1 \end{aligned}

where Y_{i} are referred to as the stage values and F_{j} as the stage derivatives. Both quantities are related by the differential equation:

F_{i} = f̂(Y_{i}).

The matrices A = [a_{ij}], U = [u_{ij}], B = [b_{ij}], V = [v_{ij}] are characteristic of a specific method. Each scheme can then be uniquely defined by a partitioned (s + r) × (s + r) matrix

[ \begin{matrix}A & U \\ B & V \end{matrix} ]

4.11.4 Matrix notation Adopting the notation:

ŷ^{[n−1]} = [ ŷ_{0}^{[n−1]} ; ŷ_{1}^{[n−1]} ; ⋮ ; ŷ_{r−1}^{[n−1]} ] , ŷ^{[n]} = [ ŷ_{0}^{[n]} ; ŷ_{1}^{[n]} ; ⋮ ; ŷ_{r−1}^{[n]} ] , Y = [ Y_{0} ; Y_{1} ; ⋮ ; Y_{s−1} ] , F = [ F_{0} ; F_{1} ; ⋮ ; F_{s−1} ]

the general linear method can be written more compactly in the following form:

[ \begin{matrix}Y \\ ŷ & [n] \end{matrix} ] = [ \begin{matrix}A ⊗ I_{N} & U ⊗ I_{N} \\ B ⊗ I_{N} & V ⊗ I_{N} \end{matrix} ] [ \begin{matrix}ΔtF \\ ŷ & [n−1] \end{matrix} ]

where I_{N} is the identity matrix of dimension N × N.

4.11.5 General Linear Methods in Nektar++

Although the GLM is essentially presented for ODE’s in its autonomous form, in Nektar++ it will be used to solve ODE’s formulated in non-autonomous form. Given the ODE,

\frac{dy}{dt} = f (t, y), y(t_{0}) = y_{0}

a single step of GLM can then be evaluated in the following way: input

y^{[n−1]} , i.e. the r subvectors comprising y_{i}^{[n−1]} - t^{[n−1]} , i.e. the equivalent of y^{[n−1]} for the time variable t.

step 1: The stage values Y_{i}, T_{i} and the stage derivatives F_{i} are calculated through the relations:

Y_{i} = Δt ∑ _{j=0} ^{s−1} a_{ij} F_{j} + ∑ _{j=0} ^{r−1} u_{ij} y_{j}^{[n−1]} , i = 0, 1, x, s − 1
T_{i}x = xΔt ∑ _{j=0} ^{s−1} a_{ij} + ∑ _{j=0} ^{r−1} u_{ij} t_{j}^{[n−1]} , i = 0, 1, x, s − 1
F_{i}x = xf(T_{i}, Y_{i}), i = 0, 1, x, s − 1

step 2: The approximation at the new time level y^{[n]} is calculated as a linear combination of the stage derivatives F_{i} and the input vector y^{[n−1]} . In addition, the time vector t^{[n]} is also updated

y_{i}^{[n]} = Δt ∑ _{j=0} ^{s−1} b_{ij} F_{j} + ∑ _{j=0} ^{r−1} v_{ij} y_{j}^{[n−1]} , i = 0, 1, x, r − 1
t_{i}^{[n]}x = xΔt ∑ _{j=0} ^{s−1} b_{ij} + ∑ _{j=0} ^{r−1} v_{ij} t_{j}^{[n−1]} , i = 0, 1, x, r − 1

output

y^{[n]}, i.e. the r subvectors comprising y_{i}^{[n]}. y_{0}^{[n]} corresponds to the actual approximation at the new time level.

t^{[n]} where t_{0}^{[n]} is equal to the new time level t + Δt.

For a detailed description of the formulation and a deeper insight of the numerical method see 4.11.6 Types of time Integration Schemes

Nektar++ contains various classes and methods which implement the concept of GLMs. This toolbox is capable of numerically solving the generalised ODE using a broad range of different time-stepping methods. We distinguish the following types of general linear methods:

Formally Explicit Methods: These types of methods are considered explicit from an ODE point of view. They are characterised by a lower triangular coefficient matrix A, ”i.e.” a_{ij} = 0 for j ≥ i. To avoid confusion, we make a further distinction:

direct explicit method: Only forward operators are required.

indirect explicit method: The inverse operator is required.

Diagonally Implicit Methods: Compared to explicit methods, the coefficient matrix A has now non-zero entries on the diagonal. This means that each stage value depend on the stage derivative at the same stage, requiring an implicit step. However, the calculation of the different stage values is still uncoupled. Best known are the DIRK schemes.

IMEX schemes: These schemes support the concept of being able to split right hand forcing term into an explicit and implicit component. This is useful in advection diffusion type problems where the advection is handled explicitly and the diffusion is handled implicit.

Fully Implicit Methods Methods: The coefficient matrix has a non-zero upper triangular part. The calculation of all stages values is fully coupled.

The aim in Nektar++ is to fully support the first three types of GLMs. Fully implicit methods are currently not implemented.

4.11.7 Usage

The goal of abstracting the concept of general linear methods is to provide users with a single interface for time-stepping, independent of the chosen method. The classes tree allows the user to numerically integrate ODE’s using high-order complex schemes, as if it was done using the Forward-Euler method. Switching between time-stepping schemes is as easy as changing a parameter in an input file.

In the already implemented solvers the time-integration schemes have been set up according to the nature of the equations. For example the incompressible Navier-Stokes equations solver allows the use of three different Implicit-Explicit time-schemes if solving the equations with a splitting-scheme. This is because this kind of scheme has an explicit and an implicit operator that combined solve the ODE’s system.

Once aware of the problem’s nature and implementation, the user can easily switch between some (depending on the problem) of the following time-integration schemes:

11 AdamsBashforthOrder1  11 Adams-Bashforth Forward multi-step scheme of order 1
11 AdamsBashforthOrder2  11 Adams-Bashforth Forward multi-step scheme of order 2
11 AdamsBashforthOrder3  11 Adams-Bashforth Forward multi-step scheme of order 3
11 AdamsMoultonOrder1  11 Adams-Moulton Forward multi-step scheme of order 1
11 AdamsMoultonOrder2  11 Adams-Moulton Forward multi-step scheme of order 2
11 BDFImplicitOrder1  11 BDF multi-step scheme of order 1 (implicit)
11 BDFImplicitOrder2  11 BDF multi-step scheme of order 2 (implicit)
11 ClassicalRungeKutta4  11 Runge-Kutta multi-stage scheme 4th order explicit (old name)
11 RungeKutta4  11 Runge-Kutta multi-stage scheme 4th order explicit (new name for ClassicalRungeKutta4)
11 RungeKutta3_SSP  11 Nonlinear SSP RungeKutta3 explicit
11 RungeKutta2_ImprovedEuler 11 Improved RungeKutta2 explicit (old name meaning Heun’s method)
11 RungeKutta2_SSP  11 Nonlinear  SSP  RungeKutta2  explicit  (surrogate  for RungeKutta2_ImprovedEuler)
11 ForwardEuler  11 Forward-Euler scheme
11 BackwardEuler  11 Backward Euler scheme
11 IMEXOrder1  11 IMEX 1st order scheme using Euler Backwards Euler Forwards
11 IMEXOrder2  11 IMEX 2nd order scheme using Backward Different Formula & Extrapolation
11 IMEXOrder3  11 IMEX 3rd order scheme using Backward Different Formula & Extrapolation
11 Midpoint  11 Midpoint method (old name)
11 RungeKutta2  11 Classical RungeKutta2 method (new name for Midpoint)
11 DIRKOrder2  11 Diagonally Implicit Runge-Kutta scheme of order 2
11 DIRKOrder3  11 Diagonally Implicit Runge-Kutta scheme of order 3
11 CNAB  11 Crank-Nicolson-Adams-Bashforth Order 2 (CNAB)
11 IMEXGear  11 IMEX Gear Order 2
11 MCNAB  11 Modified Crank-Nicolson-Adams-Bashforth Order 2 (MCNAB)
11 IMEXdirk_1_1_1  11 Forward-Backward Euler IMEX DIRK(1,1,1)
11 IMEXdirk_1_2_1  11 Forward-Backward Euler IMEX DIRK(1,2,1)
11 IMEXdirk_1_2_2  11 Implicit-Explicit Midpoint IMEX DIRK(1,2,2)
11 IMEXdirk_2_2_2  11 L-stable, two stage, second order IMEX DIRK(2,2,2)
11 IMEXdirk_2_3_2  11 L-stable, three stage, third order IMEX DIRK(3,4,3)
11 IMEXdirk_2_3_3  11 L-stable, two stage, third order IMEX DIRK(2,3,3)
11 IMEXdirk_3_4_3  11 L-stable, three stage, third order IMEX DIRK(3,4,3)
11 IMEXdirk_4_4_3  11 L-stable, four stage, third order IMEX DIRK(4,4,3)

11

Nektar++ input file for your problem will ask you just the string corresponding the time-stepping scheme you want to use (between quotation marks in the previous list), and few parameters to define your integration in time (time-step and number of steps or final time). For example:

   1 <SOLVERINFO>x 2 xxx<IxPROPERTY="EQTYPE"xVALUE="UnsteadyNavierStokes"x/>x 3 xxx<IxPROPERTY="SolverType"xVALUE="VelocityCorrectionScheme"x/>x 4 xxx<IxPROPERTY="AdvectionForm"xVALUE="Convective"x/>x 5 xxx<IxPROPERTY="Projection"xVALUE="Galerkin"x/>x 6 xxx<IxPROPERTY="TimeIntegrationMethod"xVALUE="IMEXOrder1"x/>x 7 x</SOLVERINFO>x 8 x 9 x<PARAMETERS>x 10 xxx<P>xTimeStepx=x0.001xxxxx</P>x 11 xxx<P>xNumStepsx=x1000xxxxxx</P>x 12 xxx<P>xKinvisx=x1xxxxxxxxxxx</P>x 13 </PARAMETERS>

4.11.8 Implementation of a time-dependent problem

In order to implement a new solver which takes advantage of the time-integration class in Nektar++, two main ingredients are required:

A main files in which the time-integration of you ODE’s system is initialized and performed.

A class representing the spatial discretization of your problem, which reduces your system of PDE’s to a system of ODE’s.

Your pseudo-main file, where you actually loop over the time steps, will look like

1 NekDoublextimestepxxxx=x0.1;x 2 NekDoublextimexxxxxxxx=x0.0;x 3 intxNumStepsxxxxxxxxxx=x1000;x 4 x 5 YourClassxsolver(Input);x 6 x 7 LibUtilities::TimeIntegrationMethodxxTIME_SCHEME;x 8 LibUtilities::TimeIntegrationSchemeOperatorsxODE;x 9 x 10 ODE.DefineOdeRhs(&YourClass::YourExplicitOperatorFunction,solver);x 11 ODE.DefineProjection(&YourClass::YourProjectionFunction,solver);x 12 ODE.DefineImplicitSolve(&YourClass::YourImplicitOperatorFunction,solver);x 13 x 14 Array<OneD,xLibUtilities::TimeIntegrationSchemeSharedPtr>xIntScheme;x 15 x 16 LibUtilities::TimeIntegrationSolutionSharedPtrxode_solution;x 17 x 18 Array<OneD,xArray<OneD,xNekDouble>x>xU;x 19 x 20 TIME_SCHEMEx=xLibUtilities::eForwardEuler;x 21 intxnumMultiSteps=1;x 22 IntSchemex=xArray<OneD,LibUtilities::TimeIntegrationSchemeSharedPtr>(numMultiSteps);x 23 LibUtilities::TimeIntegrationSchemeKeyxIntKey(TIME_SCHEME);x 24 IntScheme[0]x=xLibUtilities::TimeIntegrationSchemeManager()[IntKey];x 25 ode_solutionx=xIntScheme[0]->InitializeScheme(timestep,U,time,ODE);x 26 x 27 for(intxnx=x0;xnx<xNumSteps;x++n)x{x 28 xxxxUx=xIntScheme[0]->TimeIntegrate(timestep,ode_solution,ODE);x 29 }

We can distinguish three different sections in the code above

Definitions

1 NekDoublextimestepxxxx=x0.1;x 2 NekDoublextimexxxxxxxx=x0.0;x 3 intxNumStepsxxxxxxxxxx=x1000;x 4 x 5 YourClassxequation(Input);x 6 x 7 LibUtilities::TimeIntegrationMethodxxTIME_SCHEME;x 8 LibUtilities::TimeIntegrationSchemeOperatorsxODE;x 9 x 10 ODE.DefineOdeRhs(&YourClass::YourExplicitOperatorFunction,equation);x 11 ODE.DefineProjection(&YourClass::YourProjectionFunction,xequation);x 12 ODE.DefineImplicitSolve(&YourClass::YourImplicitOperatorFunction,xequation);

In this section you define the basic parameters (like time-step, initial time, etc.) and the time-integration objects. The operators are not all required, it depends on the nature of your problem and on the type of time integration schemes you want to use. In this case, the problem has been set up to work just with Forward-Euler, then for sure you will not need the implicit operator. An object named equation has been initialized, is an object of type YourClass, where your spatial discretization and the functions which actually represent your operators are implemented. An example of this class will be shown later in this page.

Initialisations

      1 Array<OneD,LibUtilities::TimeIntegrationSchemeSharedPtr>xIntScheme;x 2 x 3 LibUtilities::TimeIntegrationSolutionSharedPtrxode_solution;x 4 x 5 Array<OneD,xArray<OneD,xNekDouble>x>xU;x 6 x 7 TIME_SCHEMEx=xLibUtilities::eForwardEuler;x 8 intxnumMultiSteps=1;x 9 IntSchemex=xArray<OneD,LibUtilities::TimeIntegrationSchemeSharedPtr>(numMultiSteps);x 10 LibUtilities::TimeIntegrationSchemeKeyxIntKey(TIME_SCHEME);x 11 IntScheme[0]x=xLibUtilities::TimeIntegrationSchemeManager()[IntKey];x 12 ode_solutionx=xIntScheme[0]->InitializeScheme(timestep,U,time,ODE);

The second part consists in the scheme initialization. In this example we set up just Forward-Euler, but we can set up more then one time-integration scheme and quickly switch between them from the input file. Forward-Euler does not require any other scheme for the start-up procedure. High order multi-step schemes may need lower-order schemes for the start up.

Integration

      1 for(intxnx=x0;xnx<xNumSteps;x++n)x{x 2 xxxxUx=xIntScheme[0]->TimeIntegrate(timestep,ode_solution,ODE);x 3 }

The last step is the typical time-loop, where you iterate in time to get your new solution at each time-level. The solution at time t^{n+1} is stored into vector U (you need to properly initialize this vector). U is an Array of Arrays, where the first dimension corresponds to the number of variables (eg. u,v,w) and the second dimension corresponds to the variables size (e.g. the number of modes or the number of physical points).

The variable ODE is an object which contains the methods. A class representing a PDE equation (or a system of equations) must have a series of functions representing the implicit/explicit part of the method, which represents the reduction of the PDE’s to a system of ODE’s. The spatial discretization and the definition of this method should be

implemented in YourClass. &YourClass::YourExplicitOperatorFunction is a

functor, i.e. a pointer to a function where the method is implemented. equation is a pointer to the object, i.e. the class, where the function/method is implemented. Here a pseudo-example of the .h file of your hypothetical class representing the set of equations. The implementation of the functions is meant to be in the related .cpp file.

1 classxYourCalssx 2 {x 3 public:x 4 xxxxYourClass(INPUT);x 5 x 6 xxxx~YourClass(void);x 7 x 8 xxxxvoidxYourExplicitOperatorFunction(x 9 xxxxxxxxconstxArray<OneD,xArray<OneD,xNekDouble>x>x&xinarray,x 10 xxxxxxxxxxxxxxArray<OneD,xArray<OneD,xNekDouble>x>x&xoutarray,x 11 xxxxxxxxconstxNekDoublextime);x 12 x 13 xxxxvoidxYourProjectionFunction(x 14 xxxxxxxxconstxArray<OneD,xArray<OneD,xNekDouble>x>x&xinarray,x 15 xxxxxxxxxxxxxxArray<OneD,xArray<OneD,xNekDouble>x>x&xoutarray,xconstx 16 xxxxxxxxxxxxxxNekDoublextime);x 17 x 18 xxxxvoidxYourImplicitOperatorFunction(x 19 xxxxxxxxconstxArray<OneD,xArray<OneD,xNekDouble>x>x&xinarray,x 20 xxxxxxxxxxxxxxArray<OneD,xArray<OneD,xNekDouble>x>x&xoutarray,xconstx 21 xxxxxxxxxxxxxxNekDoublextime,x 22 xxxxxxxxconstxNekDoublexlambda);x 23 x 24 xxxxvoidxInternalMethod1(...);x 25 x 26 xxxxNekDoublexinternalvariale1;x 27 x 28 protected:x 29 xxxx...x 30 x 31 private:x 32 xxxx...x 33 };

4.11.9 Strongly imposed essential boundary conditions

Dirichlet boundary conditions can be strongly imposed by lifting the known Dirichlet solution. This is equivalent to decompose the approximate solution y into an known part, y^{D}, which satisfies the Dirichlet boundary conditions, and an unknown part, y^{H}, which is zero on the Dirichlet boundaries, i.e.

y = y^{D} + y^{H}

In a Finite Element discretisation, this corresponds to splitting the solution vector of coefficients y into the known Dirichlet degrees of freedom y^{D} and the unknown homogeneous degrees of freedom y^{H}. Ordering the known coefficients first, this corresponds to:

y = [ \begin{matrix}y & D \\ y & H \end{matrix} ]

The generalised formulation of the general linear method (i.e. the introduction of a left hand side operator) allows for an easier treatment of these types of boundary conditions. To better appreciate this, consider the equation for the stage values for an explicit general linear method where both the left and right hand side operator are linear operators, i.e. they can be represented by a matrix.

M Y_{i} = Δt ∑ _{j=0} ^{i−1} a_{ij} LY_{j} + ∑ _{j=0} ^{r−1} u_{ij} y_{j}^{[n−1]} , i = 0, 1, x, s − 1

In case of a lifted known solution, this can be written as:

\begin{matrix}[ \begin{matrix}M ^{DD} & M ^{DH} \\ M ^{HD} & M ^{HH} \end{matrix} ] [ \begin{matrix}Y & D i \\ Y & H i \end{matrix} ] = Δt ∑ _{j=0} ^{i−1} a_{ij} [ \begin{matrix}L^{DD} & L^{DH} \\ L^{HD} & L^{HH} \end{matrix} ] [ \begin{matrix}Y & D j \\ Y & H j \end{matrix} ] + ∑ _{j=0} ^{r−1} u_{ij} [ \begin{matrix}y & {D[n−1] \atop j} \\ y & {H[n−1] \atop j} \end{matrix} ] , \\ i = 0, 1, x, s − 1 \end{matrix}

In order to calculate the stage values correctly, the explicit operator should now be implemented to do the following:

[ \begin{matrix}b & D \\ b & H \end{matrix} ] = [ \begin{matrix}L^{DD} & L^{DH} \\ L^{HD} & L^{HH} \end{matrix} ] [ \begin{matrix}y & D \\ y & H \end{matrix} ]

Note that only the homogeneous part b^{H} will be used to calculate the stage values. This means essentially that only the bottom part of the operation above, i.e. L^{HD} y^{D} + L^{HH} y^{H} is required. However, sometimes it might be more convenient to use/implement routines for the explicit operator that also calculate b^{D}.

An implicit method should solve the system:

( [ \begin{matrix}M ^{DD} & M ^{DH} \\ M ^{HD} & M ^{HH} \end{matrix} ] − λ [ \begin{matrix}L^{DD} & L^{DH} \\ L^{HD} & L^{HH} \end{matrix} ]) [ \begin{matrix}y & D \\ y & H \end{matrix} ] = [ \begin{matrix}H ^{DD} & H ^{DH} \\ H ^{HD} & H ^{HH} \end{matrix} ] [ \begin{matrix}y & D \\ y & H \end{matrix} ] = [ \begin{matrix}b & D \\ b & H \end{matrix} ]

for the unknown vector y. This can be done in three steps:

Set the known solution y^{D}

Calculate the modified right hand side term b^{H} − H ^{HD} y^{D}

Solve the system below for the unknown y^{H}, i.e. H ^{HH} y^{H} = b^{H} − H ^{HD} y^{D}

4.11.10 How to add a new time-stepping method To add a new time integration scheme, follow the steps below:

Choose a name for the method and add it to the TimeIntegrationMethod enum list.

Populate the switch statement in the TimeIntegrationScheme constructor with the coefficients of the new method.

Use ( or modify) the function InitializeScheme to select (or implement) a proper initialisation strategy for the method.

Add documentation for the method (especially indicating what the auxiliary parameters of the input and output vectors of the multi-step method represent)

4.11.11 Examples of already implemented time stepping schemes

Here we show some examples time-stepping schemes implemented in Nektar++, to give an idea of what is required to add one of them.

Forward Euler

[ \begin{matrix}A & U \\ B & V \end{matrix} ] = [ \begin{matrix}0 & 1 \\ 1 & 1 \end{matrix} ] , y^{[n]} = [ y_{0}^{[n]} ] = [ m (t^{n}, y^{n}) ]

Backward Euler

[ \begin{matrix}A & U \\ B & V \end{matrix} ] = [ \begin{matrix}1 & 1 \\ 1 & 1 \end{matrix} ] , y^{[n]} = [ y_{0}^{[n]} ] = [ m (t^{n}, y^{n}) ]

2^{nd} order Adams-Bashforth

[ \begin{matrix}A & U \\ B & V \end{matrix} ] = [ \begin{matrix}0 & 1 & 0 \\ \frac{3}{2} & 1 & \frac{−1}{2} \\ 1 & 0 & 0 \end{matrix} ] , y^{[n]} = [ \begin{matrix}y & {[n] \atop 0} \\ y & {[n] \atop 1} \end{matrix} ] = [ \begin{matrix}m & (t^{n}, y^{n}) \\ Δtl(t^{n−1} , y^{n−1} ) \end{matrix} ]

1^{st} order IMEX Euler-Backward/ Euler-Forward

[ \begin{matrix}A^{IM} & A^{EM} & U \\ B^{IM} & B^{EM} & V \end{matrix} ] = [ \begin{matrix}[ 1 ] & [ 0 ] & [ 1 1 ] \\ [ \begin{matrix}1 \\ 0 \end{matrix} ] & [ \begin{matrix}0 \\ 1 \end{matrix} ] & [ \begin{matrix}1 & 1 \\ 0 & 0 \end{matrix} ] \end{matrix} ] with y^{[n]} = [ \begin{matrix}y & {[n] \atop 0} \\ y & {[n] \atop 1} \end{matrix} ] = [ \begin{matrix}m & (t^{n}, y^{n}) \\ Δtl (t^{n}, y^{n}) \end{matrix} ]

2^{nd} order IMEX Backward Different Formula & Extrapolation

[ \begin{matrix}A^{IM} & A^{EM} & U \\ B^{IM} & B^{EM} & V \end{matrix} ] = [ \begin{matrix}[ \frac{2}{3} ] & [ 0 ] & [ \frac{4}{3} −\frac{1}{3} \frac{4}{3} −\frac{2}{3} ] \\ [ \frac{2}{3} 0 0 0 ] & [ \begin{matrix}0 \\ 0 \\ 1 \\ 0 \end{matrix} ] & [ \begin{matrix}\frac{4}{3} & −\frac{1}{3} & \frac{4}{3} & −\frac{2}{3} \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{matrix} ] \end{matrix} ]

with

y^{[n]} = [ \begin{matrix}y & {[n] \atop 0} \\ y & {[n] \atop 1} \\ y & {[n] \atop 2} \\ y & {[n] \atop 3} \end{matrix} ] = [ \begin{matrix}m & (t^{n}, y^{n}) \\ m & (t^{n−1} , y^{n−1} ) \\ Δtl (t^{n}, y^{n}) \\ Δtl (t^{n−1} , y^{n−1} ) \end{matrix} ]

Note: The first two rows are normalised so the coefficient on y_{n}^{[n+1]} is one. In the standard formulation it is 3/2.

3rdorder IMEX Backward Different Formula & Extrapolation

[ \begin{matrix}A^{IM} & A^{EM} & U \\ B^{IM} & B^{EM} & V \end{matrix} ] = [ \begin{matrix}[ \frac{6}{11} ] & [ 0 ] & [ \frac{18}{11} −\frac{9}{11} \frac{2}{11} \frac{18}{11} −\frac{18}{11} \frac{6}{11} ] \\ [ \frac{6}{11} 0 0 0 0 0 ] & [ \begin{matrix}0 \\ 0 \\ 0 \\ 1 \\ 0 \\ 0 \end{matrix} ] & [ \begin{matrix}\frac{18}{11} & −\frac{9}{11} & \frac{2}{11} & \frac{18}{11} & −\frac{18}{11} & \frac{6}{11} \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \end{matrix} ] \end{matrix} ]

with

y^{[n]} = [ \begin{matrix}y & {[n] \atop 0} \\ y & {[n] \atop 1} \\ y & {[n] \atop 2} \\ y & {[n] \atop 3} \\ y & {[n] \atop 4} \\ y & {[n] \atop 5} \end{matrix} ] = [ \begin{matrix}m & (t^{n}, y^{n}) \\ m & (t^{n−1} , y^{n−1} ) \\ m & (t^{n−2} , y^{n−2} ) \\ Δtl (t^{n}, y^{n}) \\ Δtl (t^{n−1} , y^{n−1} ) \\ Δtl (t^{n−2} , y^{n−2} ) \end{matrix} ]

Note: The first two rows are normalised so the coefficient on y_{n}^{[n+1]} is one. In the standard formulation it is 11/6.

2^{nd} order Adams-Moulton

[ \begin{matrix}A & U \\ B & V \end{matrix} ] = [ \begin{matrix}\frac{1}{2} & 1 & \frac{1}{2} \\ \frac{1}{2} & 1 & \frac{1}{2} \\ 1 & 0 & 0 \end{matrix} ] , y^{[n]} = [ \begin{matrix}y & {[n] \atop 0} \\ y & {[n] \atop 1} \end{matrix} ] = [ \begin{matrix}m & (t^{n}, y^{n}) \\ Δtl(t^{n}, y^{n}) \end{matrix} ]

Midpoint Method

[ \begin{matrix}A & U \\ B & V \end{matrix} ] = [ \begin{matrix}0 & 0 & 1 \\ \frac{1}{2} & 0 & 1 \\ 0 & 1 & 1 \end{matrix} ] , y^{[n]} = [ y_{0}^{[n]} ] = [ m (t^{n}, y^{n}) ]

RK4: the standard fourth-order Runge-Kutta scheme

[ \begin{matrix}A & U \\ B & V \end{matrix} ] = [ \begin{matrix}0 & 0 & 0 & 0 & 1 \\ \frac{1}{2} & 0 & 0 & 0 & 1 \\ 0 & \frac{1}{2} & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 1 \\ \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} & 1 \end{matrix} ] , y^{[n]} = [ y_{0}^{[n]} ] = [ m (t^{n}, y^{n}) ]

2^{nd} order Diagonally Implicit Runge-Kutta (DIRK)

[ \begin{matrix}A & U \\ B & V \end{matrix} ] = [ \begin{matrix}λ & 0 & 1 \\ (1 − λ) & λ & 1 \\ (1 − λ) & λ & 1 \end{matrix} ] with λ = \frac{2 − \sqrt{2}}{2} , y^{[n]} = [ y_{0}^{[n]} ] = [ m (t^{n}, y^{n}) ]

3^{rd} order Diagonally Implicit Runge-Kutta (DIRK)

[ \begin{matrix}A & U \\ B & V \end{matrix} ] = [ \begin{matrix}λ & 0 & 0 & 1 \\ \frac{1}{2} (1 − λ) & λ & 0 & 1 \\ \frac{1}{4} (−6λ^{2} + 16λ − 1) & \frac{1}{4} (6λ^{2} − 20λ + 5) & λ & 1 \\ \frac{1}{4} (−6λ^{2} + 16λ − 1) & \frac{1}{4} (6λ^{2} − 20λ + 5) & λ & 1 \end{matrix} ]

with

λ = 0.4358665215, y^{[n]} = [ y_{0}^{[n]} ] = [ m (t^{n}, y^{n}) ]

3^{rd} order L-stable, three stage IMEX DIRK(3,4,3)

[ [ \begin{matrix}0 & 0 & 0 & 0 \\ 0 & λ & 0 & 0 \\ 0 & \frac{1}{2} (1 − λ) & λ & 0 \\ 0 & \frac{1}{4} (−6λ^{2} + 16λ − 1) & \frac{1}{4} (6λ^{2} − 20λ + 5) & λ \end{matrix} ] [ \begin{matrix}0 & 0 & 0 & 0 \\ λ & 0 & 0 & 0 \\ 0.3212788860 & 0.3966543747 & 0 & 0 \\ −0.105858296 & 0.5529291479 & 0.5529291479 & 0 \end{matrix} ] [ \begin{matrix}1 \\ 1 \\ 1 \\ 1 \end{matrix} ]_{[ 0 \frac{1}{4} (−6λ^{2} + 16λ − 1)  \frac{1}{4} (6λ^{2} − 20λ + 5)  λ ] [ 0 \frac{1}{4} (−6λ^{2} + 16λ − 1)  \frac{1}{4} (6λ^{2} − 20λ + 5)  λ ] [ 1 ]} ]

with λ = 0.4358665215, y^{[n]} = [ y_{0}^{[n]} ] = [ m (t^{n}, y^{n}) ] [next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] 5 Inside the Library: StdRegions

[ \begin{matrix}A^{IM} & A^{EM} & U \\ B^{IM} & B^{EM} & V \end{matrix} ] =

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

Chapterx5Inside the Library: StdRegions

In this chapter, we walk the reader through the different components of the StdRegions Directory. We begin with a discussion of the mathematical fundamentals, for which we use the book by Karniadakis and Sherwin [44] as our principle reference. We then provide the reader with an overview of the primary data structures introduced within the StdRegions Directory (often done through C++ objects), and then present the major algorithms – expressed as either object methods or functions – employed over these data structures.

The Fundamentals Behind StdRegions

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

5.1 The Fundamentals Behind StdRegions

The idea behind standard regions can be traced back to one of the fundamental principles we learn in calculus: differentiation and integrating on arbitrary regions are often difficult; however, differentiation and integration have been worked out on canonical (straight-sided or planar-sided), right-angled, coordinate aligned domains. In the case of Nektar++, we do not need to work on truly arbitrary domains, but rather on seven fundamental domains: segments in 1D, triangles and quadrilaterals in 2D, and tetrahedra, hexahedra, prisms and pyramids in 3D. Since Nektar++ deals with polynomial methods, the natural domain over which reference elements should be build is [−1, 1] as this is the interval over which Jacobi Polynomials are defined and over which Gaussian quadrature is defined [19]. We show in Figure 5.1 a pictorial representation of the standard segment, standard quadrilateral (often shorthanded quad) and standard triangle. pict Figurex5.1: Example of the 1D and 2D reference space elements (segment, quad and triangle).

The standard quad and standard hexahedra (shorthanded ’hex’) are geometrically tensor-product constructions defined on [−1, 1]^{d} for d = 2 and d = 3 respectively. The standard triangle is constructed by taking ξ_{1} ∈ [−1, 1] and taking ξ_{2} ≤ −ξ_{1}. The standard tetrahedra (shorthanded ’tet’) is built upon the standard triangle and has all four faces being triangles, with the two triangles along the coordinate directions looking like the standard triangle. The standard prism consists of a standard triangle along the ξ_{1} − −ξ_{2} plane extruded into the third direction (yielding three quadrilateral faces). The standard pyramid consists of a standard quadrilateral at the base with four triangular faces reaching up to its top vertex We show this pictorially in Figure 5.2.

Figurex5.2: Hexahedron to tetrahedron transformation showing how to go from an element on [−1, 1]^{3} (i.e. a standard hexahedron) to the three other element types which are subsets on [−1, 1]^{3}. Image taken from [44].

Regardless of the particular element type, we use the fact one can build polynomial spaces over these geometric objects. In the case of triangles and tetrahedra, these are polynomial spaces in the mathematical sense, i.e. P(k) spaces. In the case of quadrilaterals and hexahedra, these are bi– and tri-polynomial spaces, i.e. Q(k). In the case of prisms and pyramids, the spaces are more complicated and are a mixture in some sense of these spaces, but they are indeed polynomial in form. This allows us to define differentiation exactly and to approximate integrals exactly up to machine precision.

In what is to follow, we describe (1) how you can build differentiation and integration operators over standard regions by mapping them to a domain over which you can exploit index separability; and then (2) how you can build differentiation and integration operators natively over various standard region shapes. The former allows one to benefit from tensor-product operators, and is at present the primary focus of all Nektar++ operators. The latter as implemented for nodal element types is currently an is-a class definition extended from the (tensor-product-based) standard element definitions.

5.1.1 Reference Element Transformations That Facilitate Separability

Differentiation and integration over the standard elements within Nektar++, in general, always try to map things back to tensor-product (e.g. tensor-contraction indexing). This allows one to transform operators that would normally have iterators that go from i = 0, x, N^{d} where d is the dimension of the shape to operators constructed with d · N operations (i.e., the product of one-dimensional operations).

As an example (to help the reader gain an intuitive understanding), consider the integration of a function f(x⃗) over a hexahedral element. If we had a quadrature rule in 3D that needed N^{d} points to integrate this function exactly, we would express the operation as:

. ∫_{E} f(x⃗)dx⃗ ≈ ∑_{i=0}^{N^{d}} ω_{i}f(z⃗_{i})

where E denotes our d-dimensional element and the set Q = {z⃗_{i}, ω_{i}} denotes the set of points and quadrature weights over the element E. Now, suppose that both our function f(x⃗) and our Q were separable: that is, f(x⃗) can be written as f_{1}(x_{1}) · f_{2}(x_{2}) · f_{3}(x_{3}) where x⃗ = (x_{1}, x_{2}, x_{3})^{T} and Q can we written in terms of 1D quadrature: z⃗_{i} = z_{i_{1}} ^{(1)} · z_{i_{2}} ^{(2)} · z_{i_{3}} ^{(3)} with the index handled by an index map σ defined by i = σ(i_{1}, i_{2}, i_{3}) = i_{1} + N · i_{2} + N^{2} · i_{3}, and similarly for the weights ω_{i}. We can then re-write the integral above as follows:

. \begin{aligned}∫_{E} f(x⃗)dx⃗ & = & ∫_{E} f(x_{1}) · f(x_{2}) · f(x_{3})dx_{1}dx_{2}dx_{3} \\ & ≈ & ∑_{i=0}^{N^{d}} ω_{i}f(z⃗_{i}) \\ & = & ∑_{i_{1}=0}^{N} ∑_{i_{2}=0}^{N} ∑_{i_{3}=0}^{N} ω_{σ(i_{1},i_{2},i_{3})} [f(z_{i_{1}}^{(1)}) · f(z_{i_{2}}^{(2)}) · f(z_{i_{3}}^{(3)})] \\ & = & [ ∑_{i_{1}=0}^{N} ω_{i_{1}} f(z_{i_{1}}^{(1)})] · [ ∑_{i_{2}=0}^{N} ω_{i_{2}} f(z_{i_{2}}^{(2)})] · [ ∑_{i_{3}=0}^{N} ω_{i_{3}} f(z_{i_{3}}^{(3)})] \end{aligned}

As you can see, when the functions are separable and the quadrature (or collocating points) can be written in separable form, we can take O(N^{d}) operators and transform them into O(d · N)) operators. The above discussion is mainly focusing on the mathematical transformations needed to accomplish this; in subsequent sections, we will point out the memory layout and index ordering (i.e. now σ is ordered and implemented) to gain maximum performance.

The discussion of tensor-product operations on quadrilaterals and hexahedra may seem quite natural as the element construction is done with tensor-products of the segment; however, how do you create these types of operations on triangles (and anything involving triangles such as tetrahedra, prisms and pyramids)? The main mathematical building block we use that allows such operators is a transformation due to Duffy [28], which was subsequently used by Dubiner [27] in the context of finite elements and was extended to general polyhedral types by Ainsworth [7,x6].

An image denoting the Duffy transformation is shown in Figure 5.3. On the left is an example of the right-sided standard triangle with coordinate system (ξ_{1}, ξ_{2}), and on the right is the separable tensor-product domain with coordinates (η_{1}, η_{2}). We often refer to this transformation as “collapsed coordinates" as it appears as the “collapsing" of a quadrilateral domain to a triangle. Note that the edge along η_{1} = 1 on the quadrilateral collapses down to the single point (ξ_{1}, ξ_{2}) = (1, 1). This, in general, does not cause issues with our integration over volumes as this is a point of measure zero. However, special case will need be taken when doing edge integrals. We will make it a point to highlight those places in which special care is needed. pict Figurex5.3: Duffy mapping between a right-sided triangle and a square domain.

With the 2D Duffy transformation in place, we can actually create all four 3D element types from the starting point of a hexahedra (our base tensor-product shape). One collapsing operation yields a prism. The second collapsing operation yields a pyramid. Lastly, a final collapsing operation yields a tetrahedra. A visual representation of these operations are shown in Figure 5.4. pict Figurex5.4: Duffy mapping application to generate all the 3D element types. Image taken from [44].

This is meant merely to be a summary of tensor-product operations, enough to allow us to discuss the basic data types and algorithms within Nektar++. If more details are needed, we encourage the reader to consult Karniadakis and Sherwin [44] and the references therein.

5.1.2 Reference Elements On Primitive Geometric Types

Although the primary backbone of Nektar++ is tensor-product operations across all element types through the use of collapsed coordinates, we have implemented two commonly used non-tensorial basis set definitions. These are based upon Lagrange polynomial construction from a nodal (collocating) point set. As derived element types (in the C + + sense), we have implemented NodalStdTriangle (and Tet) based upon the electrostatic points of Hesthaven [38] and the Fekete points of Taylor and Wingate [56,x57].

Although these types of elements, at first glance, might seem to be “sub-optimal" (in terms of operations), there are various reasons why people might choice to use them. For instance, the point set you use for defining the collocation points dictates the interpolative projection operator (and correspondingly the properties of that operator) – an important issue when evaluating boundary conditions, initial conditions and for some non-linear operator evaluations. Within the Nektar++ team, we started an examination of this within the paper by Kirby and Sherwin [45]. There has since been various papers in the literature (beyond the scope of this developer’s guide) that give the pros and cons of tensor-product (separable) and non-tensor-product element types.

The Fundamental Data Structures within StdRegions

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

5.2 The Fundamental Data Structures within StdRegions

In almost all object-oriented languages (which includes C + +), there exists the concepts of class attributes and object attributes. Class attributes are those attributes shared by all object instances (both immediate and derived objects) of a particular class definition, and object attributes (sometimes called data members) are those attributes whose values vary from object to object and hence help to characterize (or make unique) a particular object. In C + +, object attributes are specified a header file containing class declarations; within a class declaration, attributes are grouped by their accessibility: public attributes, protected attributes and private attributes. A detailed discussion of the nuances of these categories are beyond the scope of this guide; we refer the interested reader to the following books for further details: [55,x51]. For our purposes, the main thing to appreciate is that categories dictate access patters within the inheritance hierarchy and to the “outside” world (i.e. access from outside the object). We have summarized the relationships between the categories and their accessibility in Tables 5.1, 5.2 and 5.3

1 1These tables are based upon information provided at http://www.programiz.com/cpp-programming/public-protected-private-inheritance, accessed 6 April 2018. .

Tablex5.1: Accessibility in Public Inheritance

11  Accessibility  11 private variables 11 protected variables 11 public variables
11  Accessibility from own class?  11  yes  11  yes  11  yes
11 Accessibility from derived class? 11  no  11  yes  11  yes
11 Accessibility from 2nd derived class? 11  no  11  yes  11  yes

Tablex5.2: Accessibility in Protected Inheritance

11  Accessibility  11 private variables 11 protected variables 11 public variables
11  Accessibility from own class?  11  yes  11  yes  11  yes
11 Accessibility from derived class? 11  no  11  yes  11  yes
11                                 11                11                  11 (inherited as
11                                 11                11                  11 protected variable)
11 Accessibility from 2nd derived class? 11  no  11  yes  11  yes

Tablex5.3: Accessibility in Private Inheritance

11  Accessibility  11 private variables 11 protected variables 11 public variables
11  Accessibility from own class?  11  yes  11  yes  11  yes
11 Accessibility from derived class? 11  no  11  yes  11  yes
11                                 11                11 (inherited as 11 (inherited as
11                                 11                11 private variable) 11 private variable)
11 Accessibility from 2nd derived class? 11  no  11  no  11  no

Within the StdRegions directory of the library, there exists a class inheritance hierarchy designed to try to encourage re-use of core algorithms (while simultaneously trying to minimize duplication of code). We present this class hierarchy in Figure 5.5. pict Figurex5.5: Class hierarchy derived from StdExpansion, the base class of the StdRegions Directory.

As is seen in Figure 5.5, the StdRegions hierarchy consists of three levels: the base level from which all StdRegion objects are derived is StdExpansion. This object is then specialized by dimension, yielding StdExpansion0D, StdExpansion1D, StdExpansion2D and StdExpansion3D. The dimension-specific objects are then specialized based upon shape.

The object attributes (variables) at various levels of the hierarchy can be understood in light of Figure 5.6. At its core, an expansion is a means of representing a function over a canonically-defined region of space evaluated at a collection of point positions. The various data members hold information to allow all these basic building blocks to be specified.

PIC

Figurex5.6: Diagram to help understand the various data members (object attributes) contained within StdRegions and how they connect with the mathematical representation presented earlier.

The various private, protected and public data members contained within StdRegions are provided in the subsequent sections.

5.2.1 Variables at the Level of StdExpansion

Private: There are private methods but no private data members within StdExpansion.

Protected:

Array of Basis Shared Pointers: m_base

Integer element id: m_elmt_id

Integer total number of coefficients in the expansion: m_ncoeffs

     Matrix Manager: m_stdMatrixManager
     Matrix Manager: m_stdStaticCondMatrixManager

IndexKeyMap Matrix Manager: m_IndexMapManager Public: There are public methods but no public data members within StdExpansion.

5.2.2 Variables at the Level of StdExpansion$D for various Dimensions

Private: There are private methods but no private data members within StdExpansion$D.

Protected:

0D and 1D: std::map<int, NormalVector> m_vertexNormals

2D: Currently does not have any data structure. It should probably have

     m_edgeNormals

3D: std::map<int, NormalVector> m_faceNormals

3D: std::map<int, bool> m_negatedNormals

Public: There are public methods but no public data members within StdExpansion$D.

5.2.3 Variables at the Level of Shape-Specific StdExpansions

Private:

Protected:

Public:

5.2.4 General Layout of the Basis Functions in Memory

5.2.5 General Layout

Basis functions are stored in a 1D array indexed by both mode and quadrature point. The fast index runs over quadrature points while the slow index runs over modes. This was done to match the memory access pattern of the inner product, which is the most frequently computed kernel for most solvers.

Bases are built from the tensor product of three different equation types (subsequently called Type 1, Type II and Type III respectively):

. ϕ_{p}(z) = \{ \begin{matrix}\frac{1−z}{2} & p = 0 \\ \frac{1+z}{2} & p = 1 \\ ( \frac{1−z}{2} ) ( \frac{1+z}{2} ) P_{p−1}^{1,1} (z) & 2 ≤ p < P \end{matrix}

(5.1)

. ϕ_{pq}(z) = \{ \begin{matrix}ϕ_{q}(z) & p = 0 & 0 ≤ q < P \\ ( \frac{1−z}{2} ) p & 1 ≤ p < P & q = 0 \\ ( \frac{1−z}{2} )^{p} ( \frac{1+z}{2} ) P_{q−1}^{2p−1,1} (z) & 1 ≤ p < P, & 1 ≤ q < P − p \end{matrix}

.ϕ_{pqr}(z) = \{ \begin{matrix}ϕ_{qr} & p = 0 & 0 ≤ q < P & 0 < r < P − q \\ ( \frac{1−z}{2} ) p+q & 1 ≤ p < P, & 0 ≤ q < P − p, & r = 0 \\ ( \frac{1−z}{2} )^{p+q} ( \frac{1+z}{2} ) P_{r−1}^{2p+2q−1,r} (z) & 1 ≤ p < P, & 0 ≤ q < P − p, & 1 ≤ r < P − p − qr. \end{matrix}

Here, P is the polynomial order of the basis and P_{p}^{α,β} are the pth order jacobi polynomial.

A Note Concerning Adjustments For C_{0} Continuity

Before going further it is worth reviewing the spatial shape of each node. The term \frac{1+z}{2} is an increasing function which is equal to zero at z = −1 and equal to one at z = 1. Similarly, \frac{1−z}{2} is a decreasing function which is equal to one at z = −1 and equal to zero at z = +1. These two functions are thus non-zero at one of each of the boundaries. If we need to maintain C_{0} continuity with surrounding elements (as we do in the continuous Galerkin method), then these local modes must be assembled together with the correct local modes in adjacent elements to create a continuous, global mode. For instance \frac{1+z}{2} in the left element would be continuous with \frac{1−z}{2} in the right element. The union of these two modes under assembly form a single “hat” function. By contrast, functions of the form

. _{\frac{1 − z}{2} \frac{1 + z}{2}}

are zero at both end points z = ±1. As a result, they are trivially continuous with any other function which is also equal to zero on the boundary. These “bubble” functions may be treated entirely locally and thus are used to construct the interior modes of a basis. Only bases with p > 1 have interior modes.

All of this holds separately in one dimension. Higher dimensional bases are constructed via the tensor product of 1D basis functions. As a result, we end up with a greater number of possibilities in terms of continuity. When the tensor product is taken between two bubble functions from different bases, the result is still a bubble function. When the tensor product is taken between a hat function and a bubble function we get “edge” modes (in both 2D and 3D). These are non-zero along one edge of the standard domain and zero along the remaining edges. When the tensor product is taken between two hat functions, they form a “vertex” mode which reaches its maximum at one of the vertices and is non-zero on two edges. The 3D bases are constructed similarly.

Based upon this convention, the 1D basis set consists of vertex modes and bubble modes. The 2D basis function set consists of vertex modes, edge modes and bubble modes. The 3D basis set contains vertex modes, edge modes, face modes and bubble modes.

5.2.6 2D Geometries

Quadrilateral Element Memory Organization

Quads have the simplest memory organization. The quadrilateral basis is composed of the tensor product of two Type I functions ϕ_{p}(ξ_{0,i})ϕ_{q}(ξ_{1,j}). This function would then be indexed as

      1 basis0[p*nq0x+xi]x*xbasis1[q*nq1x+xj]

where nq<b> is the number of quadrature points for the bth basis. Unlike certain mode orderings (e.g. Karniadakis and Sherwin [44]), the two hat functions are accessed as the first and second modes in memory with interior modes placed afterward. Thus,

      1 basis[i],xbasis[nqx+xi]

correspond to \frac{1−z}{2} and \frac{1+z}{2} respectively.

Triangle Element Memory Organization

Due to the use of collapsed coordinates, triangular element bases are formed via the tensor product of one basis function of Type I, and one of Type II, i.e. ϕ_{p}(η_{0,i}ϕ_{p}q(η1, j)). Since ϕ_{p} is also a Type I function, its memory ordering is identical to that used for quads. The second function is complicated by the mixing of ξ_{0} and ξ_{1} in the construction of η_{1}.

In particular, this means that the basis function has two modal indices, p and q. While p can run all the way to P, The number of q modes depends on the value of the p index q index such that 0 ≤ q < P − p. Thus, for p = 0, the q index can run all the way up to P. When p = 1, it runs up to P − 1 and so on. Memory is laid out in the same way starting with p = 0. To access all values in order, we write

      1 modex=x0x 2 forxpxinxP:x 3 xxxxforxqxinxPx-xp:x 4 xxxxxxxxout[mode*nqx+xq]x=xbasis0[p*nq]*basis1[mode*nqx+xq]x 5 xxxxmodex+=xP-p

Notice the use of the extra “mode” variable. Since the maximum value of q changes with p, basis1 is not simply a linearized matrix and instead has a triangular structure which necessitates keeping track of our current memory location.

The collapsed coordinate system introduces one extra subtlety. The mode

. ϕ_{1}(η_{1})ϕ_{1}(η_{2})

represents the top right vertex in the standard basis. However, when we move to the standard element basis, we are dealing with a triangle which only has three vertices. During the transformation, the top right vertex collapses into the top left vertex. If we naively construct an operators by iterating through all of our modes, the contribution from this vertex to mode Φ_{01} will not be included. To deal with this, we add its contribution as a correction when computing a kernel. The correction is Φ_{01} = ϕ_{0}ϕ_{01} + ϕ_{1}ϕ_{10} for a triangle.

5.2.7 3D Geometries Hexahedral Element Memory Organization

The hexahedral element does not differ much from the quadrilateral as it is the simply the product of three Type I functions.

. Φ_{pqr} = ϕ_{p}(ξ_{0})ϕ_{q}(ξ_{1})ϕ_{r}(ξ_{2}).

Prismatic Element Memory Organization

Cross sections of a triangular prism yield either a quad or a triangle based chosen direction. The basis, therefore, looks like a combination of the two different 2D geometries.

. Φ_{pqr} = ϕ_{p}(η_{0})ϕ_{q}(ξ_{1})ϕ_{pr}(η_{1}).

Taking ϕ_{p}ϕ_{pr} on its own produces a triangular face while taking ϕ_{p}ϕ_{q} on its own produces a quadrilateral face. When the three basis functions are combined into a single array (as in the inner product kernel), modes are accessed in the order p,q,r with r being the fastest index and p the slowest. The access pattern for the prism thus looks like

      1 mode_pqrx=x0x 2 mode_prx=x0x 3 forxpxinxP:x 4 xxxxforxqxinxQ:x 5 xxxxxxxxforxrxinxPx-xp:x 6 xxxxxxxxxxxxout[mode_pqr*nqx+xr]x=xbasis0[p*nq]*basis1[q*nq]*basis2[mode_prx+xr]x 7 xxxxxxxxmode_pqrx+=xPx-xpx 8 xxxxmode_prx+=xPx-xp

As with the triangle, we have to deal with complications due to collapsed coordinates. This time, the singular vertex from the triangle occurs along an entire edge of the prism. Our correction must be added to a collection of modes indexed by q . Φ_{0q1}+ = ϕ_{1}ϕ_{q}ϕ_{10}.

Tetrahedral Element Memory Organization

The tetrahedral element is the most complicated of the element constructions. It cannot simply be formed as the composition of multiple triangles since η_{2} is constructed by mixing three coordinate directions. We thus need to introduce our first Type III function.

. Φ_{pqr}(η_{0}, η_{1}, η_{2}) = ϕ_{p}(η_{0})ϕ_{pq}(η_{1})ϕ_{pqr}(η_{2}).

The r index is constrained by both p and q indices. It runs from P − p − q to 1 in a similar manner to the Type II function. Our typical access pattern is thus

      1 mode_pqrx=x0x 2 mode_pqx=x0x 3 forxpxinxP:x 4 xxxxforxqxinxPx-xp:x 5 xxxxxxxxforxrxinxPx-xpx-xr:x 6 xxxxxxxxxxxxout[mode_pqr*nqx+xr]x=xbasis0[p*nq]*basis1[mode_pqx+xq]*basis2[mode_pqrx+xr]x 7 xxxxxxxxmode_pqrx+=x(Px-xpx-xr)x 8 xxxxmode_pqx+=x(Px-xp)

The tetrahedral element also has to add a correction due to collapsed coordinates. Similar to the prism, the correction must be applied to an entire edge indexed by r

. Φ_{01r}+ = ϕ_{1}ϕ_{1}ϕ_{11r}.

Pyramidic Element Memory Organization

Like the tetrahedral element, a pyramid contains a collapsed coordinate direction which mixes three standard coordinates from the standard region. Unlike the tetrahedra, the collapse only occurs along one axis. Thus it is constructed from two Type I functions and one Type III function

. Φ_{pqr} = ϕ_{p}(η_{1})ϕ_{q}(η_{2})ϕ_{pqr}(η_{3}).

The product ϕ_{p}ϕ_{q} looks like the a quad construction which reflects the quad which serves as the base of the pyramid. A typical memory access looks like

      1 mode_pqrx=x0x 2 forxpxinxP:x 3 xxxxforxqxinxPx-xp:x 4 xxxxxxxxforxrxinxPx-xpx-xr:x 5 xxxxxxxxxxxxout[mode_pqr*nqx+xr]x=xbasis0[p*nq]*basis1[q*nq]*basis2[mode_pqr*nqx+xr]x 6 xxxxxxxxmode_pqrx+=x(Px-xpx-xr)

The Fundamental Algorithms within StdRegions

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

5.3 The Fundamental Algorithms within StdRegions

As stated in the introduction, this section of this guide is structured in question-answer form. This is not meant to capture every possible question asked of us on the Nektar++ users list; however, this set of (ever-growing) questions are meant to capture the “big ideas” that developers want to know about how StdRegions work and how they can be used.

In this section, we will through question and answer format try to cover the following basic algorithmic concepts that are found within the StdRegions part of the library:

xx

With the big ideas in place, let us now start into our questions.

[next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] 6 Inside the Library: SpatialDomains

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

Chapterx6Inside the Library: SpatialDomains

In this chapter, we walk the reader through the different components of the SpatialDomains Directory. We begin with a discussion of the mathematical fundamentals, for which we use the book by Karniadakis and Sherwin [44] as our principle reference. We then provide the reader with an overview of the primary data structures introduced within the SpatialDomains Directory (often done through C++ objects), and then present the major algorithms – expressed as either object methods or functions – employed over these data structures.

The SpatialDomains Directory and its corresponding class definitions serve two principal purposes:

To hold the elemental geometric information (i.e. vertex information, curve information and reference-to-world mapping information); and

To facility reading in and writing out geometry-related information.

When designing Nektar++, developing a class hierarchy for StdRegions (those fundamental domains over which we define integration and differentiation) and LocalRegions (i.e. elements in world-space) was fairly straightforward following [44]. For instance, a triangle in world-space is-a standard triangle. The first question that arose was where to store geometric information, as information within the LocalRegions element or as information encapsulated from the element so that multiple Expansions could all point to the same geometric information. The decision we made was to store geometric information – that is, the vertex information in world-space that defines an element and the edge and face curvature information – in its own data structure that could be shared by multiple Expansions (functions) over the same domain (element) in world-space. Hence SpatialDomains started as the directory containing Geometry and GeomFactors class definitions to meet the first item listed above. A LocalRegion is-a StdRegion and has-a SpatialDomain (i.e. Geometry and GeomFactors).

We then realized that in order to jump-start the process of constructing elements and combining them together into MultiRegions (collections of elements that represent a (sub)-domain of interest), we needed devise a light-weight data structure into which we could load geometric information from our geometry file and from which we could then construct Expansions (with their mappings, etc.). The light-weight data structure we devised was MeshGraph, and it was meant to meet the second item listed above.

xx5.2.4 General Layout of the Basis Functions in Memory

xx5.2.5 General Layout

xx5.2.6 2D Geometries

xx5.2.7 3D Geometries

x5.3 The Fundamental Algorithms within StdRegions

The Fundamentals Behind SpatialDomains

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

6.1 The Fundamentals Behind SpatialDomains

As mentioned in our discussions of the fundamentals of StdRegions (i.e. Section 5.1), one of the most fundamental tools from calculus that we regularly employ is the idea of mapping from general domains to canonical domains. General domains are the regions in world space over which we want to solve engineering problems, and thus want to be able to take derivatives and compute integrals. But as we learned in calculus, it is often non-trivial to accomplish differentiation or integration over these regions. We resort to the mapping arbitrary domains back to canonical domains over which we can define various operations. We introduced our canonical domains, which we call standard regions, in Chapter 5. We refer to our world space regions as local regions, which we will present in Chapter 7. How these two are connected are via SpatialDomains.

As will be further discussed in Section 6.2, there are two fundamental purposes served by SpatialDomains: (1) holding basic geometric information (e.g. vertex values and curvature information) and (2) holding geometric factors information. The former information relates to the geometric way we map standard regions to local regions. The latter information relates to how we use this map to allow us to accomplish differentiation and integration of functions in world space via operations on standard regions with associated map (geometric) information. In this section, we will highlight the important mathematical principles that are relevant to this section. We will first discuss the mapping itself: vertex and curvature information and how it is used. We will then discuss how geometric factor information is computed. We break this down into two subsections following the convention of the code. We will first discussed what is labeled in the code as Regular, and denotes mappings between elements of the same dimension (i.e. standard region triangles to 2D triangles in world space). Although our notation will be slightly different, we will use [44] as our guide; we refer the interested reader there for a more in-depth discussion of these topics. We will then discuss what is labeled in the code as Deformed, and denotes the mappings between standard region elements and their world space variants in a higher embedded dimension (i.e. standard region triangles to triangles lying on a surface embedded in a 3D space representing a manifold). Since the manifold work within Nektar++ was introduced as an area of research, we will use [18] as our guide. The notation therein is slightly different than that of [44] because of the necessity to use broader coordinate system transformation principles (e.g. covariance and contravariance of vectors, etc.). We will abbreviate the detailed mathematical derivations here, but encourage the interested reader to review [18] and references therein as needed. For those unfamiliar with covariant and contravariant spaces, we encourage the reader to review [13].

6.1.1 Vertex and Curvature Mapping Information

When we load in a mesh into Nektar++, elements are often described in world space based upon their vertex positions. In traditional FEM formats, this can be as simple as a list of d-dimensional vertex coordinates, followed by a list of element definitions: each row holding four integers (in the case of tetrahedra) denoting the four vertices in the vertex list that comprise an element. Nektar++ uses an HTML-based geometry file with a more rich definition of the basic geometric information that just described; we encourage developers and users to review our User Guide(s) for the organization and conventions used within our geometry files. For the purposes of this section, the important pieces of information are as follows. Let us assume that for each element, we have through our MeshGraph data structure (described in the next section) access to the vertex positions of an element. In general, each vertex v⃗_{j} is a n-tuple of dimension d denoting the dimension in which the points are specified in world space. For example, when considering a quadrilateral on the 2D plane, our vertex points each contain two coordinates denotes the x– and y–coordinates. As shown in Figure 6.1, we can express the mapping between a world space quadrilateral and the standard region quadrilateral on [−1, 1] × [−1, 1] via a bi-linear mapping function χ(ξ⃗). pict Figurex6.1: Reference space to world space mapping of a 2D quadrilateral to a straight-sided (2D) quadrilateral in the plane via a bi-linear (i.e., Q(1)) mapping.

In the case of segments, this mapping function χ(ξ⃗) is a function of one variable and is merely an affine mapping. In two dimensions, the mapping of a straight-sided triangle is a linear mapping (i.e., P(1) in the language of traditional finite elements – a total degree k = 1 space) and the mapping of a straight-sided quadrilateral is a bi-linear mapping (i.e, Q(1) in the language of traditional finite elements – a degree k = 1 map along each coordinate direction combined through tensor-product). In three dimensions, the mapping of a planar-sided tetrahedron is also a linear mapping, the mapping of a planar-sided hexahedron is a tri-linear mapping, and the prism and pyramid are mathematically somewhere in-between these two canonical types as given in [44]. The key point is that in the case of straight-sided (or planar-sided) elements, the mapping between reference space and world space can be deduced solely based upon the vertex positions. Furthermore in these cases, as denoted in Figure 6.1, the form of the mapping function is solely determined by type (shape) of the element. If only planar-sided elements are used, pre-computation involving the mapping functions can be done so that when vertex value information is available, all the data structures can be finalized.

As presented in ??, there are many components of Nektar++ that capitalize on the geometric nature of the basis functions we use. We often speak in terms of vertex modes, edge modes, face modes and internal modes – i.e., the coefficients that provide the weighting of vertex basis functions, edge basis functions, etc. It is beyond the scope of this developer guide to go into all the mathematical details of their definitions, etc. However, we do want to point out a few common developer-level features that are important. In the case of straight-sided (planar-sided) elements, the aforementioned mapping functions can be fully described by vertex basis functions. The real benefit of this approach (of connecting the mapping representation with a geometric basis) is seen when moving to curved elements.

Consider Figure 6.2 in which we modify the example given above to accommodate on curved edge. From the mathematical perspective, we know that the inclusion of this (quadratic) edge will require our mapping function to now be in Q(2). If we were not to use the fact that our basis is geometric in nature, we would be forced to form a Vandermonde system for a set of coefficients used to combine the tensor-product quadratic functions (nine basis functions in all), and use the five pieces of information available to us (the four vertex values and the one point c⃗_{0} that informs the curve on edge 1. As shown in Figure 6.2, we would expect that this updated (to accommodate a curved edge) mapping function to consist of the bi-linear mapping function with an additional term C_{1}(ξ_{1}, ξ_{2}) that encompasses the new curvature information. pict Figurex6.2: Reference space to world space mapping of a 2D quadrilateral to a curve (2D) quadrilateral in the plane via a bi-quadratic (i.e., Q(2)) mapping.

The basis we use, following [44], allows us to precisely specify C_{1}(ξ_{1}, ξ_{2}) using the edge basis function associated with edge 1, and to use the point value c⃗_{0} to specify the coefficient to be used. In the figure, we assume that the form of the function is collocating, but in practice it need not be so.

In practice, edge (and face) information can be given either as a set of point positions in world space that correspond to a particular point distribution in the reference element (i.e., evenly-spaced points or GLL points) or modal information corresponding to the geometric basis we use internally. Our geometric file formats assume the former – that curve information is provided to us as physical values at specified positions from which we infer (calculate) the modal values and store these values within SpatialDomain data structures.

Assuming we now have, for each element, a way of specifying the mapping function χ(ξ⃗), we can now move to how we compute the geometric factors for regular as well as for deformed mappings. We use the term “regular mappings" to refer to mappings from reference space elements to world space elements with the same dimension. For instance, a collection of triangles that lie in a plane can be represented by vertex coordinates that an only 2-tuples and the mappings between the reference space triangle the and world space triangle does not need to consider a larger embedded dimension. This was a fundamental assumption of the original Nektar code: segments lived on the 1D line, triangles and quadrilaterals lived on the 2D plane and that hexahedra, tetrahedra, etc., lived in the 3D volume. When redesigning Nektar++, we purposefully enabled geometric entities to live in a high-dimensional embedding space, different from their parameterized dimension. For instance, in Nektar++, it is possible to define segment expansions (functions that live over one dimensional parameterized curves) that are embedded in 3D. The same is true for quadrilaterals and triangles – although their parameterized dimension is two, both may live in a higher dimensional embedding space and thus represent a manifold of co-dimension one in that space. For mappings that maintain the co-dimension is the opposite of the dimension (i.e., quadrilateral in a plane represented by vertices with two coordinates mapped to a two parameter reference quadrilateral), we keep to the mapping conventions originally outlined in [44] and denote these mapping operations in the code by the enumerated value Regular; for mappings in which the co-dimension is greater than zero, we following the modified convention outlined in

[18] and denote these mapping operations in the code by the enumerated value Deformed.

6.1.2 Regular Mappings: Geometric Factors

Following [44], we assume that the vertex positions of an element in world space are given by x⃗ = x_{i} and that our reference space coordinates are given by ξ⃗ = ξ_{j}, there i and j run from zero to dim − 1 where dim is the parameter dimension of the element (e.g. for a triangle, dim = 2).

The Jacobian (matrix) is correspondingly defined as:

. J = \frac{∂x_{i}}{∂ξ_{j}} .

(6.1)

Note that this matrix is always square, but also note that it is not always constant across an element. Only in special cases such as the linear mapping of triangles and tetrahedra does the Jacobian matrix reduce to a constant (matrix) over the entire element. In general, this matrix can be evaluated at any point over the element for which it is constructed.

There are two (high-level) times in which this information is needed: when computing derivatives and when computing integrals. When computing derivatives, we employ the chain rule for differentiation, which in Einstein notation is given by the following expression:

. \frac{∂}{∂x_{i}} = \frac{∂ξ_{j}}{∂x_{i}} \frac{∂}{∂ξ_{j}} .

Note that this expression requires the reciprocal of the expression above – that is, J^{−1} . The polynomial mappings we use in Nektar++ are defined in terms of their forward mappings (reference to world). If the determinant of the mapping is non-zero, the inverse of the mapping exists but is not available analytically (except is special cases). As a consequence, we in general limit the places at which we compute the inverse of the Jacobian. Typically, the quadrature point positions are the places at which you need these values (since it as at these points we take physical space derivates and then use integration rules to construct weak form operators). Thus, procedurally, we do the following:

For a given element, compute the Jacobian matrix using the expression given in Equation 6.2 for each quadrature point position on the element (or for any points positions within an element at which it is needed).

Explicitly form the inverse of the matrix at each point position.

Because we accomplish the inversion of the Jacobian matrix at particular positions, this introduces an approximation to this computation. Although the inverse Jacobian matrix can be computed exactly at each point – when we then correspondingly use this information in the inner product over an element – we are in effect assuming that we are using a polynomial interpolative projection of this operator. The approximation error is introduced at the point of our quadrature approximation. Although the forward mapping is polynomial and hence we could find a polynomial integration rule and number of points/weights to integrate our bilinear forms exactly, our use of the inverse mapping in our bilinear forms means that we can only approximate our integrals. From our experiments, the impact of his is negligible in most cases, and only becomes in a concern in highly curved geometries. In such cases, over-integration might be requires to minimize the errors introduced due this approximation.

The other place at which we need the Jacobian matrix is to compute its determinant to be used in integration. The determinant of the Jacobian matrix (sometimes also called “the Jacobian" of the mapping) provides us the scaling of the metric terms used in integration. In all our computations, we assume that the determinant of the Jacobian matrix is strictly positive. In the area of mesh generation, the value of the determinant is used to estimate how good or bad the quality of the mapping is – in effect, if you have reasonable elements in your mesh. Negative Jacobian elements are inadmissible but even elements with small Jacobian determinants might cause issues. At the level of the library and solvers, we assume that these issues have been addressed by the user prior to attempting to run Nektar++ solvers and interpret their results.

6.1.3 Deformed Mappings: Geometric Factors

Following [18], we again assume that the vertex positions of an element in world space are given by x⃗ = x_{i} and that our reference space coordinates are given by ξ⃗ = ξ_{j}, there i runs from zero to dim − 1 where dim is the embedding dimension and j runs from zero to M − 1 where M is the parameter dimension of the element (e.g.for a triangle on a manifold embedded in 3D with vertex values in 3D, dim = 3 whilst M = 2).

The Jacobian (matrix) is correspondingly defined as:

. J = J_{j}^{i} = \frac{∂x_{i}}{∂ξ_{j}} .

(6.2)

In this case, we use the notational conventions of [8] which delineate covariant and contravariant forms. In general, this matrix is not square, also also note that it is not always constant across an element. In general, this matrix can be evaluated at any point over the element for which it is constructed. The metric tensor related to this transformation can be formed as:

. g = JJ^{T}

and the Jacobian factor associated with this mapping is then given by:

. g = detg.

Because various mappings are necessary when dealing with covariant and contravariant vectors, we have encapsulated all these routines into the directory GlobalMapping (see Chapter 10). At this stage, we do not implement general Piola transformations [50] that further respect H(div) and H(curl) constraints on these mappings as would be needed in solid mechanics or electromagnetics; however, there is nothing inherent within the Nektar++ framework that would preclude someone from adding these additional features as necessary.

The Fundamental Data Structures within SpatialDomains

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

6.2 The Fundamental Data Structures within SpatialDomains

As mentioned earlier, in almost all object-oriented languages (which includes C + +), there exists the concepts of class attributes and object attributes. For a summary of attributes and access patterns, please review Section 5.2. Within the SpatialDomains directory of the library, there exists a class inheritance hierarchy designed to try to encourage re-use of core algorithms (while simultaneously trying to minimize duplication of code). We present this class hierarchy in Figure 6.3. pict Figurex6.3: Class hierarchy derived from Geometry, the base class of the SpatialDomains Directory.

At its core, the items contained within SpatialDomains are meant to represent the mapping of StdRegion information into world-space. The various attributes contained herein related to this geometric (mesh, curvature and mapping) information. The various private, protected and public data members contained within StdRegions are provided in the subsequent sections.

6.2.1 Variables at the Level of Geometry

Private:

Protected:

Public:

6.2.2 Variables at the Level of Geometry$D for various Dimensions

Private:

Protected:

Public: 6.2.3 Variables at the Level of Shape-Specific Geometry Information Private: Protected: Public: 6.2.4 Reference to World-Space Mapping Geometry GeomFactors 6.2.5 MeshGraph The Fundamental Algorithms within SpatialDomains

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

6.3 The Fundamental Algorithms within SpatialDomains

As stated in the introduction, this section of this guide is structured in question-answer form. This is not meant to capture every possible question asked of us on the Nektar++ users list; however, this set of (ever-growing) questions are meant to capture the “big ideas” that developers want to know about how SpatialDomains work and how they can be used.

In this section, we will through question and answer format try to cover the following basic algorithmic concepts that are found within the SpatialDomains part of the library:

xx

With the big ideas in place, let us now start into our questions.

[next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] 7 Inside the Library: LocalRegions

[next] [prev] [prev-tail] [tail] [up] Chapterx7Inside the Library: LocalRegions In this chapter, we walk the reader through the different components of the LocalRegions Directory. We begin with a discussion of the mathematical fundamentals, for which we use the book by Karniadakis and Sherwin [44] as our principle reference. We then provide the reader with an overview of the primary data structures introduced within the LocalRegions Directory (often done through C++ objects), and then present the major algorithms – expressed as either object methods or functions – employed over these data structures.

xx6.2.1 Variables at the Level of Geometry xx6.2.2 Variables at the Level of Geometry$D for various Dimensions xx6.2.3 Variables at the Level of Shape-Specific Geometry Information xx6.2.4 Reference to World-Space Mapping xx6.2.5 MeshGraph x6.3 The Fundamental Algorithms within SpatialDomains

The Fundamentals Behind LocalRegions

[next] [prev] [prev-tail] [tail] [up] 7.1 The Fundamentals Behind LocalRegions The idea behind local regions is strongly connected to that of standard regions, but from the top-down perspective. As an example: in standard regions, we only had to consider one type of triangle, the one that is straight-sided, right-angled, and whose principle horizontal and vertical sides where aligned with the coordinate axes. Of course, meshes of elements consist of elements that may are may not be right-angled, planar-sided, etc. The starting point for us is the question of how to build build basis functions that exist of a world-space element – that is, an element whose vertex positions lie in the engineering (PDE) coordinate system of interest. Such an expansion is a local region. In Nektar++, each local region is-a standard region and has-a spatial domain data structure. The local region inherits common expansion methods from its standard region parent, and it uses its spatial domain information to specialize its operators to its local coordinate system.

The Fundamental Data Structures within LocalRegions

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

7.2 The Fundamental Data Structures within LocalRegions

As mentioned earlier, in almost all object-oriented languages (which includes C + +), there exists the concepts of class attributes and object attributes. For a summary of attributes and access patterns, please review Section 5.2. Within the LocalRegions directory of the library, there exists a class inheritance hierarchy designed to try to encourage re-use of core algorithms (while simultaneously trying to minimize duplication of code). We present this class hierarchy in Figure 7.1. pict Figurex7.1: Class hierarchy derived from Expansion, the base class of the LocalRegions Directory.

As is seen in Figure ??, the LocalRegions hierarchy consists of three levels: the base level from which all LocalRegion objects are derived is Expansion. This object is then specialized by dimension, yielding Expansion0D, Expansion1D, Expansion2D and Expansion3D. The dimension-specific objects are then specialized based upon shape.

The object attributes (variables) at various levels of the hierarchy can be understood in light of Figure 5.6. At its core, an expansion is a means of representing a function over a world-space region evaluated at a collection of point positions. The various data members hold information to allow all these basic building blocks to be specified. Many of the attributes are inherited from StdRegions as they are not unique to LocalRegions; however, each LocalRegion Expansion is uniquely defined based upon its geometric factors (which it stores via SpatialDomain information).

PIC

Figurex7.2: Diagram to help understand the various data members (object attributes) contained within LocalRegions and how they connect with the mathematical representation presented earlier. Recall that a LocalRegion is-a StdRegion and has-a SpatialDomain.

The various private, protected and public data members contained within LocalRegions are provided in the subsequent sections.

7.2.1 Variables at the Level of Expansion

Private: There are private methods but no private data members within Expansion.

Protected: As discussed above, the primary data in LocalRegions that distinguishes it from StdExpansions is the has-a relationship with SpatialDomains, given by the following:

SpatialDomains::GeometrySharedPtr m_geom

SpatialDomains::GeomFactorsSharedPtr m_metricinfo

     MetricMap m_metrics

Public: There are public methods but no public data members within Expansion.

7.2.2 Variables at the Level of Expansion$D for various Dimensions Private: Protected: Public:

7.2.3 Variables at the Level of Shape-Specific Expansions Private: Protected: Public:

The Fundamental Algorithms within LocalRegions

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

7.3 The Fundamental Algorithms within LocalRegions

As stated in the introduction, this section of this guide is structured in question-answer form. This is not meant to capture every possible question asked of us on the Nektar++ users list; however, this set of (ever-growing) questions are meant to capture the “big ideas” that developers want to know about how LocalRegions work and how they can be used.

In this section, we will through question and answer format try to cover the following basic algorithmic concepts that are found within the LocalRegions part of the library:

xx

With the big ideas in place, let us now start into our questions.

[next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] 8 Inside the Library: Collections

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

Chapterx8Inside the Library: Collections

In this chapter, we walk the reader through the different components of the Collections Directory. We begin with a discussion of the mathematical fundamentals, for which we use the book by Karniadakis and Sherwin [44] as our principle reference. We then provide the reader with an overview of the primary data structures introduced within the Collections Directory (often done through C++ objects), and then present the major algorithms – expressed as either object methods or functions – employed over these data structures.

The Fundamentals Behind Collections

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

8.1 The Fundamentals Behind Collections

The concept of “collections” is that there are many operations that we (conceptually) accomplish on an element-by-element basis that can be aggregated together into a single operation. The basic idea is that whenever you see a piece of code that consists of a loop in which, in the body of the loop, you are accomplishing the same operation (e.g. function evaluation) based upon an input variable tied to the loop index – you may be able to make the operation a “collection”. For instance, consider if you were given the following snippet of code in Matlab notation:

xxxxxdoxjx=x1:N,
xxxxxxxxxxxy(j)x=xdot(a,x(:,j))
xxxxxxend

where y is a row vector of length N, x is an M × N matrix and a is a row vector of length M. Note that in this code snipped, the vector a remains constant; only the vector x involved in the dot product is changing based upon the index (i.e. selecting a column of the array). From linear algebra, you might recognize this as the way we accomplish the product of a row vector with a matrix – we march through the matrix accomplishing dot products. Although these two things are equivalent mathematically, they are not equivalent computationally from the perspective of computational efficiency. Most linear algebra operations within Nektar++ are done using BLAS – Basic Linear Algebra Subroutines – a collection of specialized routines for accomplishing Level 1 (single loop), Level 2 (double nested loop) and Level 3 (triple nested loop) operations. A dot product is an example of a Level 1 BLAS operation; A matrix-vector multiplication is an example of a Level 2 BLAS operation; and finally a matrix-matrix multiplication is an example of a Level 3 BLAS operation. The general rule of thumb is that as you move up the levels, the operations can be made more efficient (using various algorithmic reorganization and memory layout strategies). Hence when you see a loop over Level 1 BLAS calls (like in the above piece of code), you should always ask yourself if this could be transformed into a Level 2 BLAS call. Since the vector a is not changing, this piece of code should be replaced with an appropriate vector-matrix multiply.

As seen in StdRegions and LocalRegions, we often conceptually thing of many of our basic building block operations as being elemental. We have elemental basis matrices, local stiffness matrices, etc. Although it is correct to implement operations as an iterator over elements in which you accomplish an action like evaluation (done by taking a basis evaluation matrix times a coefficient vector), this operation can be made more efficient by ganging all these elemental operations together into a collection.

As mentioned earlier, one of the caveats of this approach is that you must assume some level of consistency of the operation. For instance, in the case of physical evaluations, you must assume that a collection consists of elements that are all of the same time, have the same basis number and ordering, and are evaluated at the same set of points – otherwise the operation cannot be expressed as a single (basis) matrix multiplied by a matrix whose columns consist of the modes of all the elements who have joined the collective operation.

As will be seen in the data structure and algorithms sections of this discussion, these assumptions lead us to two fundamental types of collections: collections that live at the StdRegions level (i.e. collection operations that act on a group of elements as represented in reference space) and collections that live at the LocalRegions level (i.e. collection operations that act on a group of elements as represented in world space).

The Fundamental Data Structures within Collections

[next] [prev] [prev-tail] [tail] [up] 8.2 The Fundamental Data Structures within Collections The Fundamental Algorithms within Collections

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

8.3 The Fundamental Algorithms within Collections

As stated in the introduction, this section of this guide is structured in question-answer form. This is not meant to capture every possible question asked of us on the Nektar++ users list; however, this set of (ever-growing) questions are meant to capture the “big ideas” that developers want to know about how Collections work and how they can be used.

In this section, we will through question and answer format try to cover the following basic algorithmic concepts that are found within the Collections part of the library:

xx

With the big ideas in place, let us now start into our questions.

[next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] 9 Inside the Library: MultiRegions

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

Chapterx9Inside the Library: MultiRegions

In this chapter, we walk the reader through the different components of the MultiRegions Directory. We begin with a discussion of the mathematical fundamentals, for which we use the book by Karniadakis and Sherwin [44] as our principle reference. We then provide the reader with an overview of the primary data structures introduced within the MultiRegions Directory (often done through C++ objects), and then present the major algorithms – expressed as either object methods or functions – employed over these data structures.

xx7.2.1 Variables at the Level of Expansion

xx7.2.2 Variables at the Level of Expansion$D for various Dimensions

xx7.2.3 Variables at the Level of Shape-Specific Expansions

x7.3 The Fundamental Algorithms within LocalRegions

The Fundamentals Behind MultiRegions

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

9.1 The Fundamentals Behind MultiRegions

Up to now in the library, the various data structures and methods associated with standard regions, spatial domains and local regions are not specifically dictated by any particular numerical method. In fact, at this stage, they can all be viewed in light of approximation theory. With local regions in place, we have a region in world space over which we can represent an expansion (i.e. linear combination of basis functions) and form its derivatives and its integral. It is at the level of MultiRegions that we now combine two fundamental concepts: the idea of a collection of elements together to form a “global" expansion and the idea of how these (local) elements communicate (in the sense of how does one form approximations of a PDE of these collections of elements). Hence MultiRegions is important because it gives us a way of dealing with general tessellation and also because it is the first place at which we can connect to a specific numerical PDE approximation methods of choice (i.e. continuous Galerkin FEM methods, discontinuous Galerkin finite volume methods, etc.).

Because MultiRegions is both about grouping of elements in space (to form a domain) and about solving PDEs over these domains, you will find two primary collections of routines contained within the MultiRegions directory. You will find things related to collections of elements: ExpList (Expansion List), DisContField (discontinuous field) and ConField (continuous field); and you will find objects related to the linear systems formed based upon the particular numerical method once selects (i.e. GlobalLinSys, which stands for Global Linear System).

At present, the Nektar++ framework supports three types of numerical PDE discretizations for conservation laws:

Discontinuous Galerkin Methods: These weak-form (variational) methods do not require element continuity, but do put restrictions on the flux of information between elements. In general, these methods can be thought of as being in the class of finite volume (FV) methods. One feature of these methods that is often exploited computationally is that many operations can be considered as elemental. See [20,x39] and references therein for a more complete summary.

Continuous Galkerin Methods: These weak-form (variational) methods require at least C^{0} continuity. Mathematically, there have been extensions to higher levels of continuity, e.g. Isogeometric Analysis [21], these are not implemented in Nektar++ and would require further constraints on our SpatialDomain representations than we currently accommodate. In general, these methods can be thought of as being in the class of finite element methods (FEM). Although these methods are technically (mathematically) formulated as global methods, their elemental construction and compact basis types allow for many local operations. Many (but not all) of the linear system routines that are contained with the MultiRegions directory are focussed on this discretization type. See [54,x25,x44] and references therein for a more complete summary.

Flux Reconstruction Methods: These strong-form methods do not require element continuity, but like dG methods they impose restrictions on the flux of information between elements. In general, these methods can be though of as being in the class of generalized finite difference (FD) or collocating methods. See [46,x63] and references therein for a more complete summary.

Assembly Map pict Figurex9.1: Diagram to help explain assembly. UPDATE. pict Figurex9.2: Diagram to help explain assembly. UPDATE.

The Fundamental Data Structures within MultiRegions

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

9.2 The Fundamental Data Structures within MultiRegions

As mentioned earlier, in almost all object-oriented languages (which includes C + +), there exists the concepts of class attributes and object attributes. For a summary of attributes and access patterns, please review Section 5.2. Within the MultiRegions directory of the library, there exists a class inheritance hierarchy designed to try to encourage re-use of core algorithms (while simultaneously trying to minimize duplication of code). We present this class hierarchy in Figure 9.3. pict Figurex9.3: Class hierarchy derived from ExpList, the base class of the MultiRegions Directory.

At its core, the items contained within MultiRegions are meant to represent sets of LocalRegions. In the most abstract sense, an ExpList is merely a set of LocalRegion objects that may or may not have any relevance to each other. We then, as we do in other parts of the library, specialize on the dimension of the objects these sets will contain. At the subsequent levels of the hierarchy, we now involve information about we want to treat a collection of elements when evaluating them as a field. We consider a field to be the representation of a function over a (sub–)domain. A domain consists of a collection of elements that are connected (that is, they form a contiguous region in space). We think of the region of space as being tiled by elements over which expansions are build. If we consider a MultiRegion field to be a collection of these expansions with no constraints on their continuity, we arrive at the DisContField family of class definitions (which vary by dimension). For the currently available solvers within Nektar++, this level of field is used for the solution of PDEs via the discontinuous Galerkin (dG) and Flux-Reconstruction (FR) methods. If we consider a function as a collection of expansion in which we require continuity (in Nektar++, only C^{0} continuity), then we employ a further derived class called ContField (which again varies by dimension). For the currently available solvers within Nektar++, this level of field is used for the solution of PDEs via the continuous Galerkin (cG) method. Since many of the operations at the continuous field level do not rely upon the continuity of the field, we have structured the continuous MultiRegion object as with an is-a relationship with the DisContFields.

The various private, protected and public data members contained within MultiRegions are provided in the subsequent sections.

9.2.1 Variables at the Level of ExpList

Private:

Protected:

Public: 9.2.2 Variables at the Level of ExpList$D for various Dimensions Private: Protected: Public: 9.2.3 Variables at the Level of Discontinuous Field Expansions Private: Protected:

Public: 9.2.4 Variables at the Level of Continuous Field Expansions Private: Protected:

The Fundamental Algorithms within MultiRegions

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

9.3 The Fundamental Algorithms within MultiRegions

As stated in the introduction, this section of this guide is structured in question-answer form. This is not meant to capture every possible question asked of us on the Nektar++ users list; however, this set of (ever-growing) questions are meant to capture the “big ideas” that developers want to know about how MultiRegions work and how they can be used.

In this section, we will through question and answer format try to cover the following basic algorithmic concepts that are found within the MultiRegions part of the library:

xx

With the big ideas in place, let us now start into our questions.

Preconditioners

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

9.4 Preconditioners

Most of the solvers in Nektar++, including the incompressible Navier-Stokes equations, rely on the solution of a Helmholtz equation,

. ∇^{2}u(x) + λu(x) = f(x),

(9.1)

an elliptic boundary value problem, at every time-step, where u is defined on a domain Ω of N_{el} non-overlapping elements. In this section, we outline the preconditioners which are implemented in Nektar++. Whilst some of the preconditioners are generic, many are especially designed for the modified basis only.

9.4.1 Mathematical formulation

The standard spectral/hp approach to discretise (9.1) starts with an expansion in terms of the elemental modes:

. u^{δ}(x) = ∑_{n=0}^{N_{dof} −1} û_{n}Φ_{n}(x) = ∑_{e=1}^{N_{el}} ∑_{n=0}^{N_{m}^{e} −1} û_{n}^{e}ϕ_{n}^{e}(x)

(9.2)

where N_{el} is the number of elements, N_{m}^{e} is the number of local expansion modes within the element Ω^{e}, ϕ_{n}^{e}(x) is the n^{th} local expansion mode within the element Ω^{e}, x_{n}^{e} is the n^{th} local expansion coefficient within the element Ω^{e}. Approximating our solution byx(9.2), we adopt a Galerkin discretisation of equationx(9.1) where for an appropriate test space V^{δ} we find an approximate solution u^{δ} ∈ V^{δ} such that

. L (v, u) = ∫_{Ω} ∇v^{δ} · ∇u^{δ} + λv^{δ}u^{δ}dx = ∫_{Ω} v^{δ}fdx ∀v^{δ} ∈ V^{δ}

This can be formulated in matrix terms as

. Hû = f

where H represents the Helmholtz matrix, û are the unknown global coefficients and f the inner product the expansion basis with the forcing function.

C^{0} formulation

We first consider the C^{0} (i.e. continuous Galerkin) formulation. The spectral/hp expansion basis is obtained by considering interior modes, which have support in the interior of the element, separately from boundary modes which are non-zero on the boundary of the element. We align the boundary modes across the interface of the elements to obtain a continuous global solution. The boundary modes can be further decomposed into vertex, edge and face modes, defined as follows:

vertex modes have support on a single vertex and the three adjacent edges and faces as well as the interior of the element;

edge modes have support on a single edge and two adjacent faces as well as the interior of the element;

face modes have support on a single face and the interior of the element.

When the discretisation is continuous, this strong coupling between vertices, edges and faces leads to a matrix of high condition number κ. Our aim is to reduce this condition number by applying specialised preconditioners. Utilising the above mentioned decomposition, we can write the matrix equation as:

. [ \begin{matrix}H_{bb} & H_{bi} \\ H_{ib} & H_{ii} \end{matrix} ] [ û_{b} ; û_{i} ] = [ f̂_{b} ; f̂_{i} ]

where the subscripts b and i denote the boundary and interior degrees of freedom respectively. This system then can be statically condensed allowing us to solve for the boundary and interior degrees of freedom in a decoupled manor. The statically condensed matrix is given by

. [ \begin{matrix}H_{bb}− H_{bi}H_{ii}^{−1} H_{ib} & 0 \\ H_{ib} & H_{ii} \end{matrix} ] [ û_{b} ; û_{i} ] = [ f̂_{b} − H_{bi}H_{ii}^{−1} f̂_{i} ; f̂_{i} ]

This is highly advantageous since by definition of our interior expansion this vanishes on the boundary, and so H_{ii} is block diagonal and thus can be easily inverted. The above sub-structuring has reduced our problem to solving the boundary problem:

. S_{1}û = f̂_{1}

where S_{1} = H_{bb}− H_{bi}H_{ii}^{−1} H_{ib} and f̂_{1} = f̂_{b}− H_{bi}H_{ii}^{−1} f̂_{i}. Although this new system typically has better convergence properties (i.e lower κ), the system is still ill-conditioned, leading to a convergence rate of the conjugate gradient (CG) routine that is prohibitively slow. For this reason we need to precondition S_{1}. To do this we solve an equivalent system of the form:

. M^{−1} (S_{1}û − f̂_{1}) = 0

where the preconditioning matrix M is such that κ (M^{−1} S_{1}) is less than κ (S_{1}) and speeds up the convergence rate. Within the conjugate gradient routine the same preconditioner M is applied to the residual vector r̂_{k+1} of the CG routine every iteration:

.  ẑ_{k+1} = M^{−1} r̂_{k+1} .

HDG formulation

When utilising a hybridizable discontinuous Galerkin formulation, we perform a static condensation approach but in a discontinuous framework, which for brevity we omit here. However, we still obtain a matrix equation of the form . Λû = f̂.

where Λ represents an operator which projects the solution of each face back onto the three-dimensional element or edge onto the two-dimensional element. In this setting then, f̂ consists of degrees of freedom for each egde (in 2D) or face (in 3D). The overall system does not, therefore, results in a weaker coupling between degrees of freedom, but at the expense of a larger matrix system.

9.4.2 Preconditioners

Within the Nektar++ framework a number of preconditioners are available to speed up the convergence rate of the conjugate gradient routine. The table below summarises each method, the dimensions of elements which are supported, and also the discretisation type support which can either be continuous (CG) or discontinuous (hybridizable DG).

11 Name  11 Dimensions 11 Discretisations

11 Null                              11 All  11 All
11 Diagonal                          11 All  11 All
11 FullLinearSpace                   11 2/3D  11 CG
11 LowEnergyBlock                    11 3D  11 CG
11 Block                             11 2/3D  11 All
11 FullLinearSpaceWithDiagonal       11 All  11 CG
11 FullLinearSpaceWithLowEnergyBlock 11 2/3D  11 CG
11 FullLinearSpaceWithBlock          11 2/3D  11 CG

11

The default is the Diagonal preconditioner. The above preconditioners are specified through the Preconditioner option of the SOLVERINFO section in the session file. For example, to enable FullLinearSpace one can use:

   1 xx<IxPROPERTY="Preconditioner"xVALUE="FullLinearSpace"x/>

Alternatively one can have more control over different preconditioners for each solution field by using the GlobalSysSoln section. For more details, consult the user guide. The following sections specify the details for each method.

Diagonal

Diagonal (or Jacobi) preconditioning is amongst the simplest preconditioning strategies. In this scheme one takes the global matrix H = (h_{ij}) and computes the diagonal terms h_{ii}. The preconditioner is then formed as a diagonal matrix M^{−1} = (h_{ii}^{−1} ).

Linear space

The linear space (or coarse space) of the matrix system is that containing degrees of freedom corresponding only to the vertex modes in the high-order system. Preconditioning of this space is achieved by forming the matrix corresponding to the coarse space and inverting it, so that

. M^{−1} = (S_{1}^{−1} )_{vv}

Since the mesh associated with higher order methods is relatively coarse compared with traditional finite element discretisations, the linear space can usually be directly inverted without memory issues. However such a methodology can be prohibitive on large parallel systems, due to a bottleneck in communication.

In Nektar++ the inversion of the linear space present is handled using the XX^{T} library. XX^{T} is a parallel direct solver for problems of the form Ax̂ = b̂ based around a sparse factorisation of the inverse of A. To precondition utilising this methodology the linear sub-space is gathered from the expansion and the preconditioned residual within the CG routine is determined by solving

. (S_{1})_{vv}ẑ = r̂

The preconditioned residual ẑ is then scattered back to the respective location in the global degrees of freedom.

Block

Block preconditioning of the C^{0} continuous system is defined by the following:

. M^{−1} = [ \begin{matrix}(S_{1}^{−1} )_{vv} & 0 & 0 \\ 0 & (S_{1}^{−1} )_{eb} & 0 \\ 0 & 0 & (S_{1}^{−1} )_{ef} \end{matrix} ]

where diag[(S_{1})_{vv}] is the diagonal of the vertex modes, (S_{1})_{eb} and (S_{1})_{fb} are block diagonal matrices corresponding to coupling of an edge (or face) with itself i.e ignoring the coupling to other edges and faces. This preconditioner is best suited for two dimensional problems.

In the HDG system, we take the block corresponding to each face and invert it. Each of these inverse blocks then forms one of the diagonal components of the block matrix M^{−1} .

9.4.3 Low energy

Low energy basis preconditioning follows the methodology proposed by Sherwin & Casarin. In this method a new basis is numerically constructed from the original basis which allows the Schur complement matrix to be preconditioned using a block preconditioner. The method is outlined briefly in the following.

Elementally the local approximation u^{δ} can be expressed as different expansions lying in the same discrete space V^{δ}

. u^{δ}(x) = ∑_{i}^{dim(V^{δ})} û_{1i}ϕ_{1i}(x) = ∑_{i}^{dim(V^{δ})} û_{2i}ϕ_{2j} (x)

Since both expansions lie in the same space it’s possible to express one basis in terms of the other via a transformation, i.e.

. ϕ_{2} = Cϕ_{1} ⟹ û_{1} = C^{T}û_{2}

Applying this to the Helmholtz operator it is possible to show that,

. H_{2} = CH_{1}C^{T}

For sub-structured matrices (S) the transformation matrix (C) becomes:

. C = [ \begin{matrix}R & 0 \\ 0 & I \end{matrix} ]

Hence the transformation in terms of the Schur complement matrices is: . S_{2} = RS_{1}R^{T} Typically the choice of expansion basis ϕ_{1} can lead to a Helmholtz matrix that has undesirable properties i.e poor condition number. By choosing a suitable transformation matrix C it is possible to construct a new basis, numerically, that is amenable to block diagonal preconditioning.

. S_{1} = [ \begin{matrix}S_{vv} & S_{ve} & S_{vf} \\ S_{ve}^{T} & S_{ee} & S_{ef} \\ S_{vf}^{T} & S_{ef}^{T} & S_{ff} \end{matrix} ] = [ \begin{matrix}S_{vv} & S_{v,ef} \\ S_{v,ef}^{T} & S_{ef,ef} \end{matrix} ]

Applying the transformation S_{2} = RS_{1}R^{T} leads to the following matrix

. S_{2} = [ \begin{matrix}S_{vv} + R_{v}S_{v,ef}^{T} + S_{v,ef} R_{v}^{T} + R_{v}S_{ef,ef} R_{v}^{T} & [S_{v,ef} + R_{v}S_{ef,ef} ]A^{T} \\ A[S_{v,ef}^{T} + S_{ef,ef} R_{v}^{T}] & AS_{ef,ef} A^{T} \end{matrix} ]

where AS_{ef,ef} A^{T} is given by

. AS_{ef,ef} A^{T} = [ \begin{matrix}S_{ee} + R_{ef} S_{ef}^{T} + S_{ef} R_{ef}^{T} + R_{ef} S_{ff} R_{ef}^{T} & S_{ef} + R_{ef} S_{ff} \\ S_{ef}^{T} + S_{ff} R_{ef}^{T} & S_{ff} \end{matrix} ]

To orthogonalise the vertex-edge and vertex-face modes, it can be seen from the above that . R_{ef}^{T} = −S_{ff}^{−1} S_{ef}^{T} and for the edge-face modes:

. R_{v}^{T} = −S_{ef,ef}^{−1} S_{v,ef}^{T}

Here it is important to consider the form of the expansion basis since the presence of S_{ff} ^{−1} will lead to a new basis which has support on all other faces; this is problematic when creating a C^{0} continuous global basis. To circumvent this problem when forming the new basis, the decoupling is only performed between a specific edge and the two adjacent faces in a symmetric standard region. Since the decoupling is performed in a rotationally symmetric standard region the basis does not take into account the Jacobian mapping between the local element and global coordinates, hence the final expansion will not be completely orthogonal.

The low energy basis creates a Schur complement matrix that although it is not completely orthogonal can be spectrally approximated by its block diagonal contribution. The final form of the preconditioner is:

. M^{−1} = [ \begin{matrix}diag[(S_{2})_{vv}] & 0 & 0 \\ 0 & (S_{2})_{eb} & 0 \\ 0 & 0 & (S_{2})_{fb} \end{matrix} ]^{−1}

where diag[(S_{2})_{vv}] is the diagonal of the vertex modes, (S_{2})_{eb} and (S_{2})_{fb} are block diagonal matrices corresponding to coupling of an edge (or face) with itself i.e ignoring the coupling to other edges and faces.

[next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] 10 Inside the Library: GlobalMapping

[next] [prev] [prev-tail] [tail] [up] Chapterx10Inside the Library: GlobalMapping In this chapter, we walk the reader through the different components of the GlobalMapping Directory. We begin with a discussion of the mathematical fundamentals, for which we use research article by Cantwell et al. [18] and the book [8] as our principle references. We then provide the reader with an overview of the primary data structures introduced within the GlobalMapping Directory (often done through C++ objects), and then present the major algorithms – expressed as either object methods or functions – employed over these data structures.

xx9.2.2 Variables at the Level of ExpList$D for various Dimensions xx9.2.3 Variables at the Level of Discontinuous Field Expansions xx9.2.4 Variables at the Level of Continuous Field Expansions x9.3 The Fundamental Algorithms within MultiRegions x9.4 Preconditioners xx9.4.1 Mathematical formulation xx9.4.2 Preconditioners xx9.4.3 Low energy

The Fundamentals Behind GlobalMapping

[next] [prev] [prev-tail] [tail] [up] 10.1 The Fundamentals Behind GlobalMapping Based upon the appendix in [18], we outline a rigorous derivation of the Laplace-Beltrami operator. We use the convention that indices appearing once in the upper position and once in the lower position are considered dummy indices and are implicitly summed over their range, while non-repeated indices are considered free to take any value. Derivatives are also denoted using the lower-index comma notation, for example g_{ij,k}. Covariant vectors such as the gradient are those which, under a change of coordinate system, change under the same transformation in order to maintain coordinate system invariance. In contrast, contravariant vectors, such as velocity, should remain fixed under a coordinate transformation requiring that their components change under the inverse of the transformation to maintain invariance. With this in mind, we now construct the fundamental differential operators we require for a 2-dimensional manifold embedded in a 3-dimensional space. In order to express these operators in curvilinear coordinates we start by assuming that we have a smooth surface parametrization given by

x(ξ^{1}, ξ^{2}) := (x^{1}(ξ^{1}, ξ^{2}), x^{2}(ξ^{1}, ξ^{2}), x^{3}(ξ^{1}, ξ^{2})).

Next we will define the Jacobian of x as the tensor

J_{i}^{j} = \frac{∂x^{j}}{∂ξ^{i}}

where J_{i}^{j} can be viewed as a covariant surface vector (by fixing the upper index) or as a contravariant space vector (by fixing the lower index). The surface metric tensor g_{ij} can be defined in terms of the J_{i}^{j} as

g_{ij} = ∑ _{k=1} ^{3}J_{i}^{k}J_{j}^{k}.  (10.1)

which can be considered to transform a contravariant quantity to a covariant quantity. Similarly the conjugate tensor g^{ij}, which does the reverse transformation, is given by

g^{11} = g_{22}/g, g^{12} = g^{21} = −g_{12}/g, g^{22} = g_{11}/g,  (10.2)

where g is the determinant of g_{ij}. The metric tensor and its conjugate satisfy the condition

δ_{i}^{j} = g_{ik}g^{jk} = \{ \begin{matrix}1, & if i = j, \\ 0, & if i ≠ j. \end{matrix}

To construct the divergence operator we will also need the derivative of g with respect to components of the metric, g_{ij}. We know that g is invertible and from linear algebra we have that the inverse of the metric (10.2) satisfies g^{−1} = \frac{1}{g} g̃^{⊤} , where g̃ is the cofactor matrix of g. Therefore g̃^{⊤} = gg^{−1} , or in components g̃^{ij} = g(g^{−1} )_{ji}. Using Jacobi’s formula for the derivative of a matrix determinant with respect its entries, and since g is invertible, the derivative of the metric determinant is

\frac{∂g}{∂g_{ij}} = tr( g̃^{⊤} \frac{∂g}{∂g_{ij}} ) = g̃^{ij} = g(g^{−1} )_{ji} = gg^{ij} .  (10.3)

10.1.1 Divergence operator

The partial derivative of a tensor with respect to a manifold coordinate system is itself not a tensor. In order to obtain a tensor, one has to use covariant derivative, defined below. The covariant derivative of a contravariant vector is given by

∇_{k}a^{i} = a_{,k}^{i} + a^{j}Γ_{jk}^{i}.  (10.4)

where Γ_{jk}^{i} are Christoffel Symbols of the second kind. The Christoffel symbols of the first kind are defined by

Γ_{ijk} = \frac{1}{2} [g_{kj,i} + g_{ik,j} − g_{ij,k}] .

Here we note that Γ_{ijk} is symmetric in the first two indices. To obtain the Christoffel symbols of the second kind we formally raise the last index using the conjugate tensor,

Γ_{ij} ^{l} = Γ_{ijk} g^{kl}  (10.5)

which retains the symmetry in the lower two indices. We can now express the derivatives of the metric tensor in terms of the Christoffel symbols as

g_{ij,k} = Γ_{ikj} + Γ_{jki} = g_{lj} Γ_{ik}^{l} + g_{li}Γ_{jk}^{l}.

We now define the divergence operator on the manifold, ∇ · X = ∇_{k}X^{k}. Consider first the derivative of the determinant of the metric tensor g with respect to the components of some local coordinates system ξ^{1}, ξ^{2}. We apply the chain rule, making use of the derivative of the metric tensor with respect to components of the metric (10.3) and the relationship (10.5), to get

\frac{∂g}{∂ξ^{k}} = \frac{∂g}{∂g_{ij}} \frac{∂g_{ij}}{∂ξ^{k}} = gg^{ij} g_{ij,k} = gg^{ij} (Γ_{ikj} + Γ_{jki}) = g(Γ_{ik}^{i} + Γ_{jk}^{j}) = 2gΓ_{ik}^{i}.  (10.6)

We can therefore express the Christoffel symbol Γ_{ik}^{i} in terms of this derivative as

Γ_{ik}^{i} = \frac{1}{2g} \frac{∂g}{∂ξ^{k}} = \frac{1}{\sqrt{g}} \frac{∂\sqrt{g}}{∂ξ^{k}}.  (10.7)

Finally, by substituting for Γ_{ik}^{i} in the expression for the divergence operator

∇_{k}X^{k}  = X_{,k}^{k} + X^{i}Γ_{ki}^{k}

= X_{,k}^{k} + X^{k}Γ_{ik}^{i}

= X_{,k}^{k} + X^{k} \frac{1}{\sqrt{g}}(\sqrt{g})_{,k}

we can deduce a formula for divergence of a contravariant vector as

∇ · X = ∇_{k}X^{k} = \begin{matrix}(X^{k}\sqrt{g}) & ,k \\ √ & \overline{g} \end{matrix}  (10.8)

10.1.2 Laplacian operator

The covariant derivative (gradient) of a scalar on the manifold is identical to the partial derivative, ∇_{k}ϕ = ϕ_{,k}. To derive the Laplacian operator we need the contravariant form of the covariant gradient above which can be found by raising the index using the metric tensor, giving

∇^{k}ϕ = g^{kj}ϕ_{,j},  (10.9)

and substituting (10.9) for X^{k} in (10.8) to get the Laplacian operator on the manifold as

Δ_{M}ϕ = \frac{1}{\sqrt{g}} ( \sqrt{g}g^{ij} ϕ_{,j} ) _{,i}.  (10.10)

10.1.3 Anisotropic diffusion

We now extend the above operator to allow for anisotropic diffusion in the domain by deriving an expression for the surface conductivity from the ambient conductivity. The gradient of a surface function scaled by the ambient conductivity tensor ˜∇^{p} is given by

˜∇^{p}f = g^{mp}J_{m}^{l}σ_{kl}J_{j}^{k}g^{ij} \frac{∂f}{∂x^{i}} .  (10.11)

The surface gradient is mapped to the ambient space through the Jacobian J_{j}^{k}, scaled by the ambient conductivity, and mapped back to the surface through J_{m}^{l}. The anisotropic Laplacian operator is given by

˜∇^{2}f = ∇_{k} ˜∇^{k}f = ∇_{k}σ̃^{ij} ∇_{j}f  (10.12)

Therefore the surface conductivity tensor can be computed using the Jacobian tensor and the inverse metric as

σ̃ = g^{−1} JσJ^{⊤} g^{−1} .  (10.13)

10.1.4 Anisotropic Laplacian operator

Anisotropic diffusion is important in many applications. In the ambient Euclidean space, this can be represented by a diffusivity tensor σ in the Laplacian operator as

Δ_{M} = ∇ · σ∇.

On our manifold, we seek the generalisation of 10.10, in the form

Δ̃_{M}ϕ = ∇_{j} σ̃_{i}^{j}∇^{i}ϕ.

where the σ̃_{i}^{j} are entries in the surface diffusivity. For a contravariant surface vector a^{j} we can find the associated space vector A^{i} as A^{i} = J_{j}^{i}a^{j}. Similarly if A_{i} is a covariant space vector, then a_{j} = J_{j}^{i}A_{i} is a covariant surface vector. Using these we can construct the anisotropic diffusivity tensor σ̃ on the manifold by constraining the ambient diffusivity tensor σ to the surface. The contravariant surface gradient ∇^{i}ϕ is mapped to the corresponding space vector, which lies in the tangent plane to the surface. This is scaled by the ambient diffusivity and then projected back to a covariant surface vector. Finally, we use the conjugate metric to convert back to a contravariant form. The resulting surface Laplacian is

Δ̃_{M}ϕ = ∇_{m}g^{lm}J_{l}^{k}σ_{jk}J_{i}^{j}∇^{i}ϕ.

Following on from this we deduce that

σ̃_{i}^{j} = g^{jm} J_{m}^{l}σ_{lk}J_{i}^{k}.

It can be seen that in the case of isotropic diffusion that σ̃_{i}^{j} = δ_{j}^{i} ⇔ σ = I,

σ̃_{i}^{j} = g^{im}J_{m}^{k}σ_{lk}J_{j}^{l} = g^{im}J_{m}^{k}J_{j}^{k} = g^{im}g_{mj} = δ_{j}^{i}.

The Fundamental Data Structures within GlobalMapping

[next] [prev] [prev-tail] [tail] [up] 10.2 The Fundamental Data Structures within GlobalMapping The Fundamental Algorithms within GlobalMapping

[next] [prev] [prev-tail] [tail] [up] 10.3 The Fundamental Algorithms within GlobalMapping [next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] 11 Inside the Library: FieldUtils

[next] [prev] [prev-tail] [tail] [up] Chapterx11Inside the Library: FieldUtils In this chapter, we walk the reader through the different components of the FieldUtils Directory. We begin with a discussion of the mathematical fundamentals, for which we use the book by Karniadakis and Sherwin [44] as our principle reference. We then provide the reader with an overview of the primary data structures introduced within the FieldUtils Directory (often done through C++ objects), and then present the major algorithms – expressed as either object methods or functions – employed over these data structures.

xx9.2.4 Variables at the Level of Continuous Field Expansions x9.3 The Fundamental Algorithms within MultiRegions x9.4 Preconditioners xx9.4.1 Mathematical formulation xx9.4.2 Preconditioners xx9.4.3 Low energy

The Fundamentals Behind FieldUtils

[next] [prev] [prev-tail] [tail] [up] 11.1 The Fundamentals Behind FieldUtils The Fundamental Data Structures within FieldUtils

[next] [prev] [prev-tail] [tail] [up] 11.2 The Fundamental Data Structures within FieldUtils The Fundamental Algorithms within FieldUtils

[next] [prev] [prev-tail] [tail] [up] 11.3 The Fundamental Algorithms within FieldUtils [next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] 12 Inside the Library: SolverUtils

[next] [prev] [prev-tail] [tail] [up] Chapterx12Inside the Library: SolverUtils In this chapter, we walk the reader through the different components of the SolverUtils Directory. We begin with a discussion of the mathematical fundamentals, for which we use the book by Karniadakis and Sherwin [44] as our principle reference. We then provide the reader with an overview of the primary data structures introduced within the SolverUtils Directory (often done through C++ objects), and then present the major algorithms – expressed as either object methods or functions – employed over these data structures.

x9.3 The Fundamental Algorithms within MultiRegions x9.4 Preconditioners xx9.4.1 Mathematical formulation xx9.4.2 Preconditioners xx9.4.3 Low energy

The Fundamentals Behind SolverUtils

[next] [prev] [prev-tail] [tail] [up] 12.1 The Fundamentals Behind SolverUtils The Fundamental Data Structures within SolverUtils

[next] [prev] [prev-tail] [tail] [up] 12.2 The Fundamental Data Structures within SolverUtils The Fundamental Algorithms within SolverUtils

[next] [prev] [prev-tail] [tail] [up] 12.3 The Fundamental Algorithms within SolverUtils [next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] II Solvers

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

PartxIISolvers

In this part of the developer’s guide, we walk you through development aspects of the various solvers that are available within publicly-available Nektar++repository. For each of these solvers, user guides have been developed that help you see how to use these for various engineering applications. We encourage you to use these solvers in your science and engineering work flows. If you find that the current solvers do not meet your needs and consequently you end up building a new solver based upon our framework, we encourage you to consider submitting this to us for incorporation in the publicly-available master branch.

10 Inside the Library: GlobalMapping x10.1 The Fundamentals Behind GlobalMapping x10.2 The Fundamental Data Structures within GlobalMapping x10.3 The Fundamental Algorithms within GlobalMapping 11 Inside the Library: FieldUtils x11.1 The Fundamentals Behind FieldUtils x11.2 The Fundamental Data Structures within FieldUtils x11.3 The Fundamental Algorithms within FieldUtils 12 Inside the Library: SolverUtils x12.1 The Fundamentals Behind SolverUtils x12.2 The Fundamental Data Structures within SolverUtils x12.3 The Fundamental Algorithms within SolverUtils

13 ADRSolver: Solving the Advection-Reaction-Diffusion Equation

[next] [prev] [prev-tail] [tail] [up] Chapterx13ADRSolver: Solving the Advection-Reaction-Diffusion Equation In this chapter, we walk the reader through our Advection-Reaction-Diffusion Solver (ADRSolver).

[next] [prev] [prev-tail] [front] [up] 14 IncNavierStokesSolver: Solving the Incompressible Navier-Stokes Equations

[next] [prev] [prev-tail] [tail] [up] Chapterx14IncNavierStokesSolver: Solving the Incompressible Navier-Stokes Equations In this chapter, we walk the reader through our 2D and 3D incompressible Navier-Stokes Solver (IncNavierStokesSolver).

[next] [prev] [prev-tail] [front] [up] 15 CompressibleFlowSolver: Solving the Compressible Navier-Stokes Equations

[next] [prev] [prev-tail] [tail] [up] Chapterx15CompressibleFlowSolver: Solving the Compressible Navier-Stokes Equations In this chapter, we walk the reader through our 2D and 3D compressible Navier-Stokes Solver (CompressibleFlowSolver).

xx9.4.3 Low energy

Fundamental Theories of CompressibleFlowSolver

[next] [prev] [prev-tail] [tail] [up] 15.1 Fundamental Theories of CompressibleFlowSolver Data Structure of CompressibleFlowSolver

[next] [prev] [prev-tail] [tail] [up] 15.2 Data Structure of CompressibleFlowSolver Figurex15.1: CompressibleFlowSolver DataStructure

CompressibleFlowSolverSystem EulerCFE NavierStokesCFE RiemannSolver CFSBndConditions m bndConds DiffusionLDGNS m diffusion ArtificialDiffusion m artificialdiffusion m artificialdiffusion AdvectionSystem Advection m advObject m riemnn Diffusion m diffusion forcing m forcing UnsteadySystem EquationSystem TimeIntegrationScheme m intScheme Inheritance Date Footprint Factory Solver Library Flow Chart of CompressibleFlowSolver

[next] [prev] [prev-tail] [tail] [up] 15.3 Flow Chart of CompressibleFlowSolver Figurex15.2: CompressibleFlowSolver Main Flow Chart

CompressibleFlowSolver.cpp

Initialize Objects eg. DriverStandard::v InitObject

See Figure 15.3 Execute Solver beginning from driver:

eg. DriverStandard::v Execute

Initial conditions: m equ[0]- >DoInitialise

See Figure 15.4 Main execution of solver: m equ[0]- >DoSolve

See Figure 15.5 and Figure 15.6 Output brief information:

m equ[0]- >Output Figurex15.3: CompressibleFlowSolver InitObject PIC ^{CompressibleFlowSolver::drv− >GetDriverFactory Initialize the driver eg. DriverStandard::v InitObject eg. NavierStokesCFE::v InitObject CompressibleFlowSystem::v InitObject m diffusion− >InitObject AdvectionSystem::v InitObject Specific CFS advection-related object eg. RiemannSolver::v InitObject Specific CFS boundary-related object m bndConds UnsteadySystem::v InitObject EquationSystem::v InitObject TimeIntegrationSystem::v InitObject If unsteady BoundaryConditions::v InitObject} Figurex15.4: CompressibleFlowSolver Initialize Conditions PIC ^{CompressibleFlowSolver::drv− >Execute Do Initialization, eg. DriverStandard::v Execute: m equ[0]− > DoInitialise UnsteadySystem::v DoInitialise EquationSystem::SetBoundaryConditions CompressibleFlowSystem::v SetInitialConditions EquationSystem::v SetInitialConditions} Figurex15.5: CompressibleFlowSolver Execute Advection PIC ^{CompressibleFlowSolver::drv− >Execute CompressibleFlowSystem::DoOdeRhs CompressibleFlowSystem::DoAdvection Advection::Advect(){v Advect} eg. AdvectionWeakDG::v Advect Get advection flux through CompressibleFlowSystem::GetFluxVector: m fluxVector fields[i]− >IProductWRTDerivBase Solve flux at element BC: m riemann− >Solve fields[i]− >AddTraceIntegral} Figurex15.6: CompressibleFlowSolver Execute Diffusion PIC ^{CompressibleFlowSolver::drv− >Execute CompressibleFlowSystem::DoOdeRhs CompressibleFlowSystem::DoDiffusion eg. NavierStokesCFE::v DoDiffusion() {m diffusion− > Diffuse} Diffusion::Diffuse(){v Diffuse} eg. DiffusionLDGNS::v Diffuse Calculate numerical flux for first order derivitives v NumericalFluxO1 Calculate numerical flux for second order derivitives v NumericalFluxO2} [next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] III Utilities

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

PartxIIIUtilities

In this part of the developer’s guide, we walk you through development aspects of the various utilities that are available within publicly-available Nektar++repository. We encourage you to incorporate these utilities into your science and engineering work flows. If you find that the current utilities do not meet your needs and consequently you end up building a new utility based upon our framework, we encourage you to consider submitting this to us for incorporation in the publicly-available master branch.

10 Inside the Library: GlobalMapping x10.1 The Fundamentals Behind GlobalMapping x10.2 The Fundamental Data Structures within GlobalMapping x10.3 The Fundamental Algorithms within GlobalMapping 11 Inside the Library: FieldUtils x11.1 The Fundamentals Behind FieldUtils x11.2 The Fundamental Data Structures within FieldUtils x11.3 The Fundamental Algorithms within FieldUtils 12 Inside the Library: SolverUtils x12.1 The Fundamentals Behind SolverUtils x12.2 The Fundamental Data Structures within SolverUtils x12.3 The Fundamental Algorithms within SolverUtils

16 FieldConvert

[next] [prev] [prev-tail] [tail] [up] Chapterx16FieldConvert In this chapter, we walk the reader through FieldConvert.

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

[next] [prev] [prev-tail] [tail] [up] Chapterx17NekMesh In this chapter, we walk the reader through NekMesh.

[next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] IV NekPy: Python interface to Nektar++

[next] [prev] [prev-tail] [tail] [up] PartxIVNekPy: Python interface to Nektar++ x10.2 The Fundamental Data Structures within GlobalMapping x10.3 The Fundamental Algorithms within GlobalMapping 11 Inside the Library: FieldUtils x11.1 The Fundamentals Behind FieldUtils x11.2 The Fundamental Data Structures within FieldUtils x11.3 The Fundamental Algorithms within FieldUtils 12 Inside the Library: SolverUtils x12.1 The Fundamentals Behind SolverUtils x12.2 The Fundamental Data Structures within SolverUtils x12.3 The Fundamental Algorithms within SolverUtils

18 Introduction

[next] [prev] [prev-tail] [tail] [up] Chapterx18Introduction This part of the guide contains the information on using and developing NekPy Python wrappers for the Nektar++ spectral/hp element framework. As a disclaimer, these wrappings are experimental and incomplete. You should not rely on their current structure and API remaining unchanged. Currently, representative classes from the LibUtilities, StdRegions, SpatialDomains, LocalRegions and MultiRegions libraries have been wrapped in order to show the proof-of-concept.

xx10.1.3 Anisotropic diffusion xx10.1.4 Anisotropic Laplacian operator x10.2 The Fundamental Data Structures within GlobalMapping x10.3 The Fundamental Algorithms within GlobalMapping

Features and functionality

[next] [prev] [prev-tail] [tail] [up] 18.1 Features and functionality NekPy uses the Boost.Python library to provide a set of high-quality, hand-written Python bindings for selected functions and classes in Nektar++. It is worth noting that Python (CPython, the standard Python implementation written in C, in particular) includes C API and that everything in Python is strictly speaking a C structure called PyObject. Hence, defining a new class, method etc. in Python is in reality creating a new PyObject structure. Boost.Python is essentially a wrapper for Python C API which conveniently exports C++ classes and methods into PyObjects. At compilation time a dynamic library is created which is then imported to Python, as shown in Figure 18.1.

PIC Figurex18.1: A schematic diagram of how C++ code is converted into Python with Boost.Python [49] A typical snippet could look something like:

Listingx18.1: NekPy sample

snippet 1 fromxNekPy.LibUtilitiesximportxPointsKey,xPointsType,xBasisKey,xBasisTypex 2 fromxNekPy.StdRegionsximportxStdQuadExpx 3 importxnumpyxasxnpx 4 numModesx=x8x 5 numPtsxxx=x9x 6 ptsKeyxxx=xPointsKey(numPts,xPointsType.GaussLobattoLegendre)x 7 basisKeyx=xBasisKey(BasisType.Modified_A,xnumModes,xptsKey)x 8 quadExpxx=xStdQuadExp(basisKey,xbasisKey)x 9 x,xyxxxxx=xquadExp.GetCoords()x 10 fxxxxxxxx=xnp.sin(x)x*xnp.cos(y)x 11 projxxxxx=xquadExp.FwdTrans(fx)

NekPy uses the Boost.NumPy library, contained in Boost 1.63+, to automatically convert C++ Array<OneD, > objects to and from the commonly-used numpy.ndarray object, which makes the integration more seamless between Python and C++.

[next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] 19 Installing NekPy

[next] [prev] [prev-tail] [tail] [up] Chapterx19Installing NekPy NekPy has the following list of requirements:

Boost with Python support

Nektar++ master branch compiled from source (i.e. not from packages)

Python 2.7+ (note that examples rely on Python 2.7)

NumPy

Most of these can be installed using package managers on various operating systems, as we describe below. We also have a requirement on the Boost.NumPy package, which is available in Boost 1.63 or later. If this isn’t found on your system, it will be automatically downloaded and compiled.

x10.3 The Fundamental Algorithms within GlobalMapping

Compiling and installing Nektar++

[next] [prev] [prev-tail] [tail] [up] 19.1 Compiling and installing Nektar++ Nektar++ should be compiled as per the user guide instructions and installed into a directory which we will refer to as $NEKDIR. By default this is the dist directory inside the Nektar++ build directory.

Note that Nektar++ must, at a minimum, be compiled with NEKTAR_BUILD_LIBRARY,

NEKTAR_BUILD_UTILITIES , NEKTAR_BUILD_SOLVERS and NEKTAR_BUILD_PYTHON. This will

automatically download and install Boost.NumPy if required. Note that all solvers may be disabled as long as the NEKTAR_BUILD_SOLVERS option is set.

19.1.1 macOS

Homebrew

Users of Homebrew should make sure their installation is up-to-date with brew upgrade. Then run

      1 brewxinstallxpythonxboost-python

To install the NumPy package, use the pip package manager:

      1 pipxinstallxnumpy

MacPorts Users of MacPorts should sure their installation is up-to-date with sudo port selfupdate &&

sudo port upgrade outdated. Then run
      1 sudoxportxinstallxpython27xpy27-numpyx 2 sudoxportxselectx--setxpythonxpython27

19.1.2 Linux: Ubuntu/Debian Users of Debian and Ubuntu Linux systems should sure their installation is up-to-date with

sudo apt-get update && sudo apt-get upgrade

      1 sudoxapt-getxinstallxlibboost-python-devxpython-numpy

19.1.3 Compiling the wrappers Run the following command in $NEKDIR/build directory to install the Python package for the current user:

      1 makexnekpy-install-user

Alternatively, the following command can be used to install the package for all users:

1 makexnekpy-install-system

Using the bindings

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

19.2 Using the bindings

By default, the bindings will install into the dist directory, along with a number of examples that are stored in the $NEKDIR/library/Demos/Python directory. To test your installation, you can for example run one of these (e.g. python Basis.py) or launch an interactive session:

      1 $xcdxbuildsx 2 $xpythonx 3 Pythonx2.7.13x(default,xAprxx4x2017,x08:47:57)x 4 [GCCx4.2.1xCompatiblexApplexLLVMx8.1.0x(clang-802.0.38)]xonxdarwinx 5 Typex"help",x"copyright",x"credits"xorx"license"xforxmorexinformation.x 6 >>>xfromxNekPy.LibUtilitiesximportxPointsKey,xPointsTypex 7 >>>xPointsKey(10,xPointsType.GaussLobattoLegendre)x 8 <NekPy.LibUtilities._LibUtilities.PointsKeyxobjectxatx0x11005c310>

19.2.1 Examples

A number of examples of the wrappers can be found in the $NEKDIR/library/Demos/Python directory, along with a sample mesh newsquare_2x2.xml:

SessionReader.py is the simplest example and shows how to construct a session

     reader object. Run it as python SessionReader.py mesh.xml.

Basis.py shows functionality of basic LibUtilities points and basis classes. Run

     this as python Basis.py.

StdProject.py shows how to use some of the StdRegions wrappers and duplicates the functionality of Basis.py using the StdExpansion class. Run this as python

StdProject.py.

MeshGraph.py loads a mesh and prints out some basic properties of its

quadrilateral elements. Run it as python MeshGraph.py newsquare_2x2.xml.

If you want to modify the source files, it’s advisable to edit them in the $NEKDIR/library/Demos/Python directory and re-run make install, otherwise local changes will be overwritten by the next make install.

[next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] 20 Package structure

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

Chapterx20Package structure

The NekPy wrapper is designed to mimic the library structure of Nektar++, with directories for the LibUtilities, SpatialDomains and StdRegions libraries. This is a deliberate design decision, so that classes and definitions in Nektar++ can be easily located inside NekPy.

There are also some other directories and files:

     LibUtilities/Python/NekPyConfig.hpp is a convenience header that all .cpp

files should import. It sets appropriate namespaces for boost::python and boost::python::numpy, depending on whether the Boost.NumPy library was compiled or is included in Boost,

cmake/python contains templates for init.py and setup.py files which every Python package should contain,

cmake/ThirdPartyPython.cmake is a CMake configuration file which searches for Boost.Python and prepares make targets for installing NekPy.

Figure 20.1 shows the location of Python wrapper files within Nektar++ structure. Every sub-module of Nektar++ hosts an additional folder for Python files and the structure of the Python folder mimics the structure of the sub-module itself, as shown on the left of the figure with LibUtilities sub-module. Individual .cpp files, such as Expansion.cpp contain the wrappers for classes and methods from corresponding Nektar++ library files whereas .cpp files named after sub-modules (e.g. LibUtilities.cpp) contain BOOST_PYTHON_MODULE definitions which create package modules

(e.g. NekPy.LibUtilities).

PIC Figurex20.1: The location of Python wrapper files within Nektar++ structure.

[next] [prev] [prev-tail] [front] [up] 21 NekPy wrapping guide

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

Chapterx21NekPy wrapping guide

This section attempts to outline some of the basic principles of the NekPy wrapper, which relies entirely on the excellent Boost.Python library. An extensive documentation is therefore beyond the scope of this document, but we highlight aspects that are important for the NekPy wrappers.

In general, note that when things go wrong with Boost.Python, it’ll be indicated either by an extensive compiler error, or a runtime error in the Python interpreter when you try to use your wrapper. Judicious use of Google is therefore recommended to track down these issues!

You may also find the following resources useful:

The Boost.Python tutorial [22],

The Boost.Python entry on the Python wiki [2],

Examples on GitHub [5],

The OpenStreetGraph cookbook [4] and rationale for using manual wrapping [3] which served as a starting point for this project.

To demonstrate how to wrap classes, we’ll refer to a number of existing parts of the code below.

Defining a library

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

21.1 Defining a library First consider LibUtilities. An abbreviated version of the base file, LibUtilities.cpp has the following structure:

Listingx21.1: Defining a library with

            Boost.Python 1 #includex<LibUtilities/Python/NekPyConfig.hpp>x 2 x 3 voidxexport_Basis();x 4 voidxexport_SessionReader();x 5 x 6 BOOST_PYTHON_MODULE(_LibUtilities)x 7 {x 8 xxxx//xInitialisexBoost.NumPy.x 9 xxxxnp::initialize();x 10 x 11 xxxx//xExportxclasses.x 12 xxxxexport_Basis();x 13 xxxxexport_SessionReader();x 14 }

The BOOST_PYTHON_MODULE(name) macro allows us to define a Python module inside C++. Note that in this case, the leading underscore in the name (i.e. _LibUtilities) is deliberate. To define the contents of the module, we call a number of functions that are prefixed by export_, which will define one or more Python classes that live in this module. These Python classes correspond with our Nektar++ classes. We adopt this approach simply so that we can split up the different classes into different files, because it is not possible to call BOOST_PYTHON_MODULE more than once. These functions are defined in appropriately named files, for example export_Basis() lives

in the file LibUtilities/Python/Foundations/Basis.cpp. This corresponds to
the Nektar++ file LibUtilities/Foundations/Basis.cpp and the classes defined

therein.

Basic class wrapping

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

21.2 Basic class wrapping

As a very basic example of wrapping a class, let’s consider the SessionReader wrapper.

Listingx21.2: Basic class wrapping with

Boost.Python 1 voidxexport_SessionReader()x 2 {x 3 xxxxpy::class_<SessionReader,x 4 xxxxxxxxxxxstd::shared_ptr<SessionReader>,x 5 xxxxxxxxxxxboost::noncopyable>(x 6 xxxxxxxxxxxxxxx"SessionReader",xpy::no_init)x 7 x 8 xxxxxxxx.def("CreateInstance",xSessionReader_CreateInstance)x 9 xxxxxxxx.staticmethod("CreateInstance")x 10 x 11 xxxxxxxx.def("GetSessionName",x&SessionReader::GetSessionName,x 12 xxxxxxxxxxxxxpy::return_value_policy<py::copy_const_reference>())x 13 x 14 xxxxxxxx.def("Finalise",x&SessionReader::Finalise)x 15 xxxxxxxx;x 16 }

21.2.1 py::class_<>

This Boost.Python object defines a Python class in C++. It is templated, and in this case we have the following template arguments:

SessionReader is the class that will be wrapped std::shared_ptr<SessionReader> indicates that this object should be stored inside a shared (or smart) pointer, which we frequently use throughout the library, as can be seen by the frequent use of SessionReaderSharedPtr

boost::noncopyable indicates that Boost.Python shouldn’t try to automatically wrap the copy constructor of SessionReader. We add this here because of compiler errors due to subclasses used inside SessionReader, but generally, this should be used for abstract classes which can’t be copied.

We then have two arguments:

"SessionReader" is the name of the class in Python.

py::no_init indicates this object has no publically-accessible initialiser. This is because for SessionReader, we define a factory-type function called

CreateInstance instead.

21.2.2 Wrapping member functions

We then call the .def function on the class_<>, which allows us to define member functions on our class. This is equivalent to def-ing a function in Python. .def has two required parameters, and one optional parameter:

The function name as a string, e.g. "GetSessionName"

A function pointer that defines the C++ function that will be called

An optional return policy, which we need to define when the C++ function returns a reference.

Boost.Python is very smart and can convert many Python objects to their equivalent C++ function arguments, and C++ return types of the function to their respective Python object. Many times therefore, one only needs to define the .def() call.

However, there are some instances where we need to do some additional conversion, mask some C++ complexity from the Python interface, or deal with functions that return references. We describe ways to deal with this below.

Thin wrappers

Instead of defining a function pointer to a member of the C++ class, we can define a function pointer to a separate function that defines some extra functionality. This is called a thin wrapper.

As an example, consider the CreateInstance function. In C++ we pass this function the command line arguments in the usual argc, argv format. In Python, command line arguments are defined as a list of strings inside sys.argv. However, Boost.Python does not know how to convert this list to argc, argv, so we need some additional code.

Listingx21.3:

Thin wrapper

          example 1 SessionReaderSharedPtrxSessionReader_CreateInstance(py::listx&ns)x 2 {x 3 xxxx//x...xsomexcodexherexthatxconvertsxaxPythonxlistxtoxthexstandardx 4 xxxx//xc/c++x(intxargc,xcharx**argv)xformatxforxcommandxlinexarguments.x 5 xxxx//xThenxusexthisxtoxconstructxaxSessionReaderxandxreturnxit.x 6 xxxxSessionReaderSharedPtrxsrx=xSessionReader::CreateInstance(argc,xargv);x 7 xxxxreturnxsr;x 8 }

In Python, we can then simply call session = SessionReader.CreateInstance(sys.argv).

Dealing with references

When dealing with functions in C++ that return references, e.g. const NekDouble &GetFactor() we need to supply an additional argument to .def(), since Python immutable types such as strings and integers cannot be passed by reference. For a full list of options, consult the Boost.Python guide. However a good rule of thumb is to use copy_const_reference as highlighted above, which will create a copy of the const reference and return this.

Dealing with Array<OneD, >
The LibUtilities/Python/BasicUtils/SharedArray.cpp file contains a number of

functions that allow for the automatic conversion of Nektar++ Array<OneD, > storage to and from NumPy ndarray objects. This means that you can wrap functions that take these as parameters and return arrays very easily. However bear in mind the following caveats:

Any NumPy ndarray created from an Array<OneD, > (and vice versa) will share their memory. Although this avoids expensive memory copies, it means that changing the C++ array changes the contents of the NumPy array (and vice versa).

Many functions in Nektar++ return Arrays through argument parameters. In Python this is a very unnatural way to write functions. For example:

           1 #xThisxisxgoodx 2 x,xy,xzx=xexp.GetCoords()x 3 #xThisxisxbadx 4 x,xy,xzx=xnp.zeros(10),xnp.zeros(10),xnp.zeros(10)x 5 exp.GetCoords(x,y,z)

Use thin wrappers to overcome this problem. For examples of how to do this, particularly in returning tuples, consult the StdRegions/StdExpansion.cpp wrapper which contains numerous examples.

TwoD and ThreeD arrays are not supported.

More information on the memory management and how the memory is shared can be found in Sectionx23.

21.2.3 Inheritance

Nektar++ makes heavy use of inheritance, which can be translated to Python quite easily using Boost.Python. For a good example of how to do this, you can examine the StdRegions wrapper for StdExpansion and its elements such as StdQuadExp. In a cut-down form, these look like the following:

Listingx21.4: Inheritance with

      Boost.Python 1 voidxexport_StdExpansion()x 2 {x 3 xxxxpy::class_<StdExpansion,x 4 xxxxxxxxxxxxxxxstd::shared_ptr<StdExpansion>,x 5 xxxxxxxxxxxxxxxboost::noncopyable>(x 6 xxxxxxxxxxxxxxxxxxx"StdExpansion",xpy::no_init);x 7 }x 8 voidxexport_StdQuadExp()x 9 {x 10 xxxxpy::class_<StdQuadExp,xpy::bases<StdExpansion>,x 11 xxxxxxxxxxxxxxxstd::shared_ptr<StdQuadExp>x>(x 12 xxxxxxxxxxxxxxxxxxx"StdQuadExp",xpy::init<constxLibUtilities::BasisKey&,x 13 xxxxxxxxxxxxxxxxxxxconstxLibUtilities::BasisKey&>());x 14 }

Note the following:

StdExpansion is an abstract class, so it has no initialiser and is non-copyable.

We use py::bases<StdExpansion> in the definition of StdQuadExp to define

its parent class. This does not necessarily need to include the full hierarchy of C++ inheritance: in StdRegions the inheritance graph for StdQuadExp looks like

StdExpansion -> StdExpansion2D -> StdQuadExp. In the above wrapper, we

omit the StdExpansion2D call entirely.

py::init<> is used to show how to wrap a C++ constructor. This can accept any arguments for which you have either written explicit wrappers or Boost.Python already knows how to convert.

21.2.4 Wrapping enums

Most Nektar++ enumerators come in the form:

      1 enumxMyEnumx{x 2 xxxxeItemOne,x 3 xxxxeItemTwo,x 4 xxxxSIZE_MyEnumx 5 };x 6 staticxconstxcharx*MyEnumMap[]x=x{x 7 xxxx"ItemOne"x 8 xxxx"ItemTwo"x 9 xxxx"ItemThree"x 10 };

To wrap this, you can use the NEKPY_WRAP_ENUM macro defined in NekPyConfig.hpp, which in this case can be used as NEKPY_WRAP_ENUM(MyEnum, MyEnumMap). Note that if instead of const char * the map is defined as a const std::string, you can use the

NEKPY_WRAP_ENUM_STRING macro.

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

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

Chapterx22Documentation

The NekPy package certainly had to be documented in order to provide an easily accessible information about the wrapped classes to both users and developers. Ideally, the documentation should be:

easily readable by humans,

accessible using Python’s inbuilt help method,

compatible with the existing Nektar++ doxygen-based documentation.

Traditionally, Python classes and functions are documented using a docstring – a string occurring as the very first statement after the function or class is defined. This string is then accessible as the __doc__attribute of the function or class. The conventions associated with Python docstrings are described in PEP 257 document [34].

Boost.Python provides an easy way to include docstrings in the wrapped methods and classes as shown in Listing 22.1. The included docstrings will appear when Python help method is used.

Listingx22.1: Example of class and method documentation in

            Boost.Python 1 voidxexport_Points()x 2 {x 3 xxxxxxxxpy::class_<PointsKey>("PointsKey",x 4 xxxxxxxx"CreatexaxPointsKeyxwhichxuniquelyxdefinesxquadraturexpoints.\n"x 5 xxxxxxxx"\n"x 6 xxxxxxxx"Args:\n"x 7 xxxxxxxx"\tnQuadPointsx(integer):xThexnumberxofxquadraturexpoints.\n"x 8 xxxxxxxx"\tpointsTypex(PointsTypexobject):xThextypexofxquadraturexpoints.",x 9 xxxxxxxxpy::init<constxint,xconstxPointsType&>())x 10 xxxxxxxxxxxxxxxxx 11 xxxxxxxxxxxx.def("GetNumPoints",x&PointsKey::GetNumPoints,x 12 xxxxxxxxxxxx"GetxthexnumberxofxquadraturexpointsxinxPointsKey.\n"x 13 xxxxxxxxxxxx"\n"x 14 xxxxxxxxxxxx"Args:\n"x 15 xxxxxxxxxxxx"\tNone\n"x 16 xxxxxxxxxxxx"Returns:\n"x 17 xxxxxxxxxxxx"\tIntegerxdefiningxthexnumberxofxquadraturexpointsxinxPointsKey.")x 18 }

In order to fully document the existing bindings a number of enumeration type classes such as PointsType had to have docstrings included which proved to be a challenge since Boost.Python does not provide a way to do this. Instead a direct call to Python C API has to be made and the method adapted from [52] was used, as shown in Listing 22.2. A downside of this solution is that it does requires the developer to manually update the Python documentation if the enumeration type is ever changed (e.g. adding a new type of point) as the code does not automatically gather information from the C++ class. In theory it could be possible to create a Python script which would generate Python docstrings based on the existing C++ documentation using regular expressions; however it would be difficult to integrate this solution into the existing framework.

Listingx22.2:

Code used to include dosctings in enumetation type classes - part of

          NekPyConfig.hpp 1 #definexNEKPY_WRAP_ENUM_STRING_DOCS(ENUMNAME,MAPNAME,DOCSTRING)xxxx\x 2 xxxx{xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\x 3 xxxxxxxxpy::enum_<ENUMNAME>xtmp(#ENUMNAME);xxxxxxxxxxxxxxxxxxxxxxxx\x 4 xxxxxxxxforx(intxax=x0;xax<x(int)SIZENAME(ENUMNAME);x++a)xxxxxxxxxx\x 5 xxxxxxxx{xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\x 6 xxxxxxxxxxxxtmp.value(MAPNAME[a].c_str(),x(ENUMNAME)a);xxxxxxxxxxxx\x 7 xxxxxxxx}xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\x 8 xxxxxxxxtmp.export_values();xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\x 9 xxxxxxxxPyTypeObjectx*xptox=xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\x 10 xxxxxxxxxxxxreinterpret_cast<PyTypeObject*>(tmp.ptr());xxxxxxxxxxxx\x 11 xxxxxxxxPyDict_SetItemString(pto->tp_dict,x"__doc__",xxxxxxxxxxxxxx\x 12 xxxxxxxxxxxxPyString_FromString(DOCSTRING));xxxxxxxxxxxxxxxxxxxxxxx\x 13 xxxx}

There are many docstrings conventions that are popular in Python such as Epytext, reST and Google therefore a choice had to be made as to which docstring style to use. After considering the criteria which the documentation had to fulfill it was decided to use Google Python Style [53] as it is highly readable by humans (and hence an excellent choice for documentation which will be primarily accessible by Python help method) and can be used to generate automated documentation pages with Sphinx (a tool for creating Python documentation).

Unfortunately, it proved to be difficult to include the documentation of NumPy package in the existing doxygen-based documentation due to the fact that the docstrings are generated by Boost.Python. It was decided that if the time constraints of the project permit this problem could be resolved at a later date and the possibility of accessing the documentation though inbuilt help method was deemed sufficient.

[next] [prev] [prev-tail] [front] [up] 23 Memory management in NekPy

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

Chapterx23Memory management in NekPy

A significant amount of effort has been invested into developing efficient memory management techniques, to enable the sharing of memory between Python and C++ for the ubiquitous Nektar++ Array structure which is used throughout the code. The computations carried out in Nektar++ are usually highly demanding on the machines they are run on, therefore any unnecessary data duplication or memory leaks would be detrimental to software performance.

This section first outlines some precusor information relating to the basics of memory management in C++ and Python, before explaining the strategy used in NekPy to effectively use the resources when passing data between the different layers of the bindings.

x18.1 Features and functionality

Memory management in C++ and Python

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

23.1 Memory management in C++ and Python

23.1.1 C++

In C++ the memory used by the application is divided into five containers [48]: text segment – containing the set of instructions for the program to execute;

data segment – containing static and global variables initialised with values, e.g. float pi = 3.14; – the size of the data segment is pre-determined at the time of the compilation of the program and depends on the size of the variables in the source code;

bss (block started by symbol) segment – containing static and global variables not explicitly initialised with any value, e.g. int n; – the size of the data segment is pre-determined at the time of the compilation of the program and depends on the size of the variables in the source code;

stack – a LIFO (last in, first out) structure containing the function calls and variables intialised in the functions – the size of the stack is pre-determined at the time of compilation of the program and depends on the operating system and development environment;

heap – containing all dynamically allocated memory (e.g. using instruction new in C++).

The access to data on the heap is maintained by pointers which store the memory address of said data. If the pointer ceases to exist (e.g. because the function which contained it finished running and hence was taken off the stack) the programmer has no way of referencing the data on the heap again. The unused, inaccessible data will stay in the memory, consuming valuable resources if not properly deallocated using delete instruction. Issues may arise if the same memory address is held by two different pointers – the deallocation of memory can render one of the pointers empty and possibly lead to errors.

In order to facilitate memory management, C++11 standard introduced shared pointers (shared_ptr)x[1]. Shared pointers maintain a reference counter which is increased when another shared pointer refers to the same memory address. The memory will only be deallocated when all shared pointers referencing the memory are destroyed, as shown in Figurex23.1.

PIC Figurex23.1: Memory management with C++11 shared_ptr. The two declared shared pointers reference an array of 4 integers. Reference counter shown in green.

Python

Python manages memory though a private heap containing all Python objectsx[30]. An internal memory manager is used to ensure the memory is properly allocated and deallocated while the script runs. In contrast to C++, Python’s memory manager tries to optimise the memory usage as much as possible. For example, the same object reference is allocated to a new variable if the object already exists in memory. Another important feature of Python’s memory manager is the use of garbage collector based on reference counting. When the object is no longer referenced by any variable the reference counter is set to zero, which triggers the garbage collector to free the memory (possibly at some later time). A disadvantage of this solution is slower execution time, since the garbage collector routines have to be called periodically. Features of the Python memory manager are schematically shown inxFigurex23.2 PIC

Figurex23.2: Memory management in Python. Memory is optimised as much as possible and garbage collector deallocates the memory if the reference counter (shown in green) reaches 0.

It is unusual for the Python programmer to manually modify the way Python uses memory resources – however it is sometimes necessary, as is the case with this project. In particular, the Python C API exposes several macros to handle reference counting, and developers can increase or decrease the reference counter of the object using Py_INCREF and Py_DECREF respectivelyx[31]. In Boost.Python, these calls are wrapped in the py::incref and

py::decref functions.

Passing C++ memory to Python

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

23.2 Passing C++ memory to Python

To highlight the technique for passing C++ memory natively into Python’s ndarray, we consider first the case of the native Nektar++ matrix structure. In many situations, matrices

created by Nektar++ (usually a shared_ptr of NekMatrix<D, StandardMatrixTag> type)

need to be passed to Python - for instance, when performing differentiation using e.g. Gauss quadrature rules a differentiation matrix must be obtained. In order to keep the program memory efficient, the data should not be copied into a NumPy array but rather be referenced by the Python interface. This, however, complicates the issue of memory management.

Consider a situation where C++ program no longer needs to work with the generated array and the memory dedicated to it is deallocated. If this memory has already been shares with Python, the Python interface may still require the data contained within the array. However since the memory has already been deallocated from the C++ side, this will typically cause an out-of-bounds memory exception. To prevent such situations a solution employing reference counting must be used.

Converter method

Boost.Python provides the methods to convert a C++ type element to one recognised by Python as well as to maintain appropriate reference counting. Listing 23.1 shows an abridged version of the converter method (for Python 2 only) with comments on individual parameters.

The object requiring conversion is a shared_ptr of NekMatrix<D, StandardMatrixTag> type

(named mat).

Listingx23.1: Converter method for converting the C++ arrays into Python NumPy

arrays. 1 #includex<LibUtilities/LinearAlgebra/NekMatrix.hpp>x 2 #includex<LibUtilities/LinearAlgebra/MatrixStorageType.h>x 3 #includex<NekPyConfig.hpp>x 4 x 5 usingxnamespacexNektar;x 6 usingxnamespacexNektar::LibUtilities;x 7 x 8 template<typenamexT,xtypenamexF>x 9 voidxNekMatrixCapsuleDestructor(voidx*ptr)x 10 {x 11 xxxx//xDestructorxforxshared_ptrxwhenxthexcapsulexisxdeallocated.x 12 xxxxstd::shared_ptr<NekMatrix<T,xF>>x*matx=x 13 xxxxxxxx(std::shared_ptr<NekMatrix<T,xF>>x*)ptr;x 14 xxxxdeletexmat;x 15 }x 16 x 17 template<typenamexT>x 18 structxNekMatrixToPythonx 19 {x 20 xxxxstaticxPyObjectx*convert(x 21 xxxxxxxxstd::shared_ptr<NekMatrix<T,xStandardMatrixTag>>xconstx&mat)x 22 xxxx{x 23 xxxxxxxx//xCreatexaxPyCObject,xwhichxwillxholdxaxshared_ptrxtoxthexmatrix.x 24 xxxxxxxx//xWhenxcapsulexisxdeallocatedxonxthexPythonxside,xitxwillxcallxthex 25 xxxxxxxx//xdestructorxfunctionxabovexandxreleasexthexshared_ptrxreference.x 26 xxxxxxxxpy::objectxcapsule(py::handle<>(PyCObject_FromVoidPtr(x 27 xxxxxxxxxxxxnewxstd::shared_ptr<NekMatrix<T,xStandardMatrixTag>>(mat),x 28 xxxxxxxxxxxxNekMatrixCapsuleDestructor<T,xStandardMatrixTag>)));x 29 x 30 xxxxxxxxintxnRowsx=xmat->GetRows(),xnColsx=xmat->GetColumns();x 31 x 32 xxxxxxxx//xincrementsxPythonxreferencexcounterxtoxavoidximmediatexdeallocationx 33 xxxxxxxxreturnxpy::incref(x 34 xxxxxxxxxxxxnp::from_data(x 35 xxxxxxxxxxxxxxxxmat->GetRawPtr(),x//xpointerxtoxdatax 36 xxxxxxxxxxxxxxxxnp::dtype::get_builtin<T>(),x//xdataxtypex 37 xxxxxxxxxxxxxxxxpy::make_tuple(nRows,xnCols),x//xshapexofxthexarrayx 38 xxxxxxxxxxxxxxxxpy::make_tuple(sizeof(T),xnRowsx*xsizeof(T)),x//xstridexofxthexarrayx 39 xxxxxxxxxxxxxxxxcapsule).ptr());x//xcapsulex-xcontainsxthexobjectxowningxthexdatax(preventingxitxfromxbeingxprematurelyxdeallocated)x 40 xxxx}x 41 }x 42 x 43 template<typenamexT>x 44 voidxexport_NekMatrix()x 45 {x 46 xxxxpy::to_python_converter<std::shared_ptr<NekMatrix<T,xStandardMatrixTag>>,x 47 xxxxxxxxxxxxxxxxxxxxxxxxxxxxNekMatrixToPython<T>>();x 48 }

Firstly we give a brief overview of the general process undertaken during type conversion. Boost.Python maintains a registry of known C++ to Python conversion types, which by default allows for fundamental data type conversions such as double and float. In this manner, many C++ functions can be automatically converted, for example when they are used in .def calls when registering a Python class. Clearly the convert function here contains much of the functionality. In order to perform automatic conversion between the NekMatrix and a 2D ndarray, we register the conversion function inside Boost.Python’s registry so that it is aware of the datatypes. We also note that throughout the conversion code (and elsewhere in NekPy), we make use of the Boost.NumPy bindings. These are a lightweight wrapper around the NumPy C API, which simplifies the syntax somewhat and avoids direct use of the API.

In terms of the conversion function itself, we first create a new Python capsule object. The capsule is designed to hold a C pointer and a callback function that is called when the Python object is deallocated. Since there is no Boost.Python wrapper around this, we use a handle<> to wrap it in a generic Boost.Python object. The strategy we therefore employ is to create a shared_ptr, which increases the reference counter of the NekMatrix. This will prevent it being destroyed if it is in use in Python, even if on the C++ side the memory is deallocated. The callback function simply deletes the shared_ptr when it is no longer required, cleaning up the memory appropriately. This process is shown in Figurex23.3. It is worth noting that the steps (c) and (d) can be reversed and the shared_ptr created by the Python binding can be removed first. In this case the memory will be deallocated only when the shared_ptr created by C++ is also removed.

PIC

(a) Nektar++ creates a data array referenced by a

shared_ptr

PIC

(b) Array is passed to Python binding which creates a new

shared_ptr to the data

PIC

(c) Nektar++ no longer needs the data - its shared_ptr is removed but the memory is not deallocated

PIC

(d) When the data is no longer needed in the Python interface the destructor is called and shared_ptr is removed.

Figurex23.3: Memory management of data created in C++ using shared_ptr and passed to Python. Reference counter shown in green.

Finally, the converter method returns a NumPy array using the np::from_data method. Note that the capsule object is passed as the ndarray base argument – i.e. the object that owns the data. In this case, when all ndarray objects that refer to the data become deallocated, NumPy will not directly deallocate the data but instead release its reference to the capsule. The capsule will then be deallocated, which will decrement the counter in the shared_ptr. We also note that data ordering is important for matrix storage; Nektar++ stores data in column-major order, whereas NumPy arrays are traditionally row-major. The stride of the array has to be passed into the np::from_data in a form of a tuple (a, b), where a is the number of bytes needed to skip to get to the same position in the next row and b is the number of bytes needed to skip to get to the same position in the next column. In order to stop Python from immediately destroying the resulting NumPy array, its reference counter is manually increased before the array is passed on to Boost.Python and eventually returned to the user’s code.

Testing

The process outlined above requires little manual intervention from the programmer. There are no almost no explicit calls to Python C API (aside from creating a PyObject – PyCObject_FromVoidPtr) as all operations are carried out by Boost.Python. Therefore, the testing focused mostly on correctness of returned data, in particular the order of the array. To this end, the Differentiation tutorials were used as tests. In order to correctly run the tutorials the Python wrapper needs to retrieve the differentiation matrix which, as mentioned before, has to be converted to a datatype Python recognises. The test runs the differentiation tutorials and compares the final result to the fixed expected value. The test is automatically run as a part of ctest command if both the Python wrapper and the tutorials have been built.

23.2.1 Passing Python data to C++

Conversely, a similar problem exists when data is created in Python and has to be passed to the C++ program. In this case, as the data is managed by Python, the main reference counter should be maintained by the Python object and incremented or decremented as appropriate using py::incref and py::decref methods respectively. Although we do not support this process for the NekMatrix as described above, we do use this process for the Array<OneD, > structure. When the array is no longer needed by the C++ program the reference counter on the Python side should be decremented in order for Python garbage collection to work appropriately – however this should only happen when the array was created by Python in the first place.

The files implementing the below procedure are:

     LibUtilities/Python/BasicUtils/SharedArray.cpp and

     LibUtilities/BasicUtils/SharedArray.hpp

Modifications to Array<OneD, const DataType> class template

In order to perform the operations described above, the C++ array structure should contain information on whether or not it was created from data managed by Python. To this end, two new attributes were added to C++ Array<OneD, const DataType> class template in the form of a struct:

      1 structxPythonInfox{x 2 xxxxvoidx*m_pyObject;x//xUnderlyingxPyObjectxpointerx 3 xxxxvoidx(*m_callback)(voidx*);x//xCallbackx 4 };

where:

m_pyObject is a pointer to the PyObject containing the data, which should be an

     ndarray;

m_callback is a function pointer to the callback function which will decrement the reference counter of the PyObject.

Inside Array<OneD, >, this struct is held as a double pointer, i.e.:

      1 PythonInfox**m_pythonInfo;

This is done because if the ndarray was created originally from a C++ to Python conversion (as outlined in the previous section), we need to convert the Array – and any other shared arrays that hold the same C++ memory – to reference the Python array. If we did not do this, then it is possible that the C++ array could be destroyed whilst it is still used on the Python side, leading to an out-of-bounds exception. By storing this as a double pointer, in a similar fashion to the underlying reference counter m_count, we can ensure that all C++ arrays are updated when necessary. We can keep track of Python arrays by checking *m_pythonInfo; if this is not set to nullptr then the array has been constructed throught the Python to C++ converter.

Adding new attributes to the arrays might cause a significantly increased memory usage or additional unforeseen overheads, although this was not seen in benchmarking. However to avoid all possibility of this, a preprocessor directive has been added to only include the additional arguments if NekPy had been built (using the option

NEKTAR_BUILD_PYTHON).

A new constructor has been added to the class template, as seen in Listing 23.2.

m_memory_pointer and m_python_decrement have been set to nullptr in the pre-existing

constructors. A similar constructor was added for const arrays to ensure that these can also be passed between the languages. Note that no calls to Nektar++ array initialisation policies are made in this constructor, unlike in the pre-existing ones, as there is no need for the new array to copy the data.

Listingx23.2: New constructor for initialising arrays created through the Python to C++ converter

method. 1 Array(unsignedxintxdim1Size,xDataType*xdata,xvoid*xmemory_pointer,xvoidx(*python_decrement)(voidx*))x:x 2 xxxxm_size(xdim1Sizex),x 3 xxxxm_capacity(xdim1Sizex),x 4 xxxxm_data(xdatax),x 5 xxxxm_count(xnullptrx),x 6 xxxxm_offset(x0x)x 7 {x 8 xxxxm_countx=xnewxunsignedxint();x 9 xxxx*m_countx=x1;x 10 x 11 xxxxm_pythonInfox=xnewxPythonInfox*();x 12 xxxx*m_pythonInfox=xnewxPythonInfo();x 13 xxxx(*m_pythonInfo)->m_callbackx=xpython_decrement;x 14 xxxx(*m_pythonInfo)->m_pyObjectx=xmemory_pointer;x 15 }

Changes have also been made to the destructor, as shown in Listingx23.3, in order to ensure that if the data was initially created in Python the callback function would decrement the reference counter of the NekPy array object. The detailed procedure for deleting arrays is described further in this section.

Listingx23.3: The modified destructor for C++

arrays. 1 ~Array()x 2 {x 3 xxxxif(xm_countx==xnullptrx)x 4 xxxx{x 5 xxxxxxxxreturn;x 6 xxxx}x 7 x 8 xxxx*m_countx-=x1;x 9 xxxxif(x*m_countx==x0x)x 10 xxxx{x 11 xxxxxxxxifx(*m_pythonInfox==xnullptr)x 12 xxxxxxxx{x 13 xxxxxxxxxxxxArrayDestructionPolicy<DataType>::Destroy(xm_data,xm_capacityx);x 14 xxxxxxxxxxxxMemoryManager<DataType>::RawDeallocate(xm_data,xm_capacityx);x 15 xxxxxxxx}x 16 xxxxxxxxelsex 17 xxxxxxxx{x 18 xxxxxxxxxxxx(*m_pythonInfo)->m_callback((*m_pythonInfo)->m_pyObject);x 19 xxxxxxxxxxxxdeletex*m_pythonInfo;x 20 xxxxxxxx}x 21 x 22 xxxxxxxxdeletexm_pythonInfo;x 23 xxxxxxxxdeletexm_count;x//xCleanxupxthexmemoryxusedxforxthexreferencexcount.x 24 xxxx}x 25 }

Creation of new arrays

The following algorithm has been proposed to create new arrays in Python and allow the C++ code to access their contents:

The NumPy array object to be converted is passed as an argument to a C++ method available in the Python wrapper.

The converter method is called to convert a Python NumPy array into C++

Array<OneD, const DataType> object.

If the NumPy array was created through the C++-to-Python process (which can be determined by checking the base of the NumPy array), then:

extract the Nektar++ Array from the capsule;

convert the Array (and all of its other references) to a Python array so that any C++ arrays that share this memory also know to call the appropriate decrement function;

set the NumPy array’s base to an empty object to destroy the original capsule.

Otherwise, the converter creates a new Array<OneD, const DataType> object with the following attribute values:

data points to the data contained by the NumPy array,

memory_pointer points to the NumPy array object,

python_decrement points to the function decrementing the reference counter

    of the PyObject,

subsequently, m_count is equal to 1.

The Python reference counter of the NumPy array object is increased.

If any new references to the array are created in C++ the m_count attribute is increased accordingly. Likewise, if new references to NumPy array object are made in Python the reference counter increases.

The process is schematically shown in Figure 23.4a and 23.4b.

Array deletion

The array deletion process relies on decrementing two reference counters: one on the Python side of the program (Python reference counter) and the other one on C++ side of the program. The former registers how many Python references to the data there are and if there is a C++ reference to the array. The latter (represented by m_count attribute) counts only the number of references on the C++ side and as soon as it reaches zero the callback function is triggered to decrement the Python reference counter so that registers that the data is no longer referred to in C++. Figure 23.5 presents the overview of the procedure used to delete the data.

In short, the fact that C++ uses the array is represented to Python as just an increment to the object reference counter. Even if the Python object goes out of scope or is explicitly deleted, the reference counter will always be non-zero until the callback function to decrement it is executed, as shown in Figure 23.4c. Similarly, if the C++ array is deleted first, the Python object will still exist as the reference counter will be non-zero (see Figure 23.4d).

Converter method

As with conversion from C++ to Python, a converter method was registered to make Python NumPy arrays available in C++ with Boost.Python, which can be found in the SharedArray.cpp bindings file. In essence, Boost.Python provides the used with a memory segment (all expressions containing rvalue_from_python are to do with doing that). The data has to be extracted from PyObject in order to be presented in a format C++ knows how to read – the get_data method allows the programmer to do it for NumPy arrays. Finally, care must be taken to manage memory correctly, thus the use of borrowed references when creating Boost.Python object and the incrementation of PyObject reference counter at the end of the method.

The callback decrement method is shown below in Listing 23.4. When provided with a pointer to PyObject it decrements it reference counter.

Listingx23.4: The decrement method called when the

m_count

of C++ array reaches

0. 1 staticxvoidxdecrement(voidx*objPtr)x 2 xxxx{x 3 xxxxxxxxPyObjectx*pyObjPtrx=x(PyObjectx*)objPtr;x 4 xxxxxxxxPy_XDECREF(pyObjPtr);x 5 xxxx}

PIC

(a) The NumPy array is created in Python. Note that the NumPy object and the data it contains are represented by two separate memory addresses.

PIC

(b) The C++ array is created through the converter method: its attributes point to the appropriate memory addresses and the reference counter of the memory address of NumPy array object is incremented.

PIC

(c) If the NumPy object is deleted first, the reference counter is decremented, but the data still exists in memory.

PIC

(d) If the C++ array is deleted first, the callback function decrements the reference counter of the NumPy object but the data still exists in memory.

Figurex23.4: Diagram showing the process of creation and deletion of arrays. Reference counters shown in green.

PIC Figurex23.5: Flowchart describing the procedure for array deletion from either side of the code.

Testing

As the process of converting arrays from Python to C++ required making direct calls to C API and relying on custom-written methods, more detailed testing was deemed necessary. In order to thoroughly check if the conversion works as expected, three tests were conducted to determine whether the array is:

referenced (not copied) between the C++ program and the Python wrapper,

accessible to C++ program after being deleted in Python code,

accessible in Python script after being deleted from C++ objects.

Python files containing test scripts are currently located in library\Demos\Python\tests. They should be converted into unit tests that should be run when Python components are built.

[next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up] 24 FieldConvert in NekPy

[next] [prev] [prev-tail] [tail] [up] Chapterx24FieldConvert in NekPy This chapter will describe the idea behind the FieldConvert utility in NekPy and discuss how the process is implemented through a Python FieldConverter class. As this part of the project has not been fully completed yet, this chapter also outlines the tasks that need to be done in order to show the proof-of-concept, as well as some ideas for further improving the tool.

x21.2 Basic class wrapping

xx21.2.1 py::class_<>

xx21.2.2 Wrapping member functions xx21.2.3 Inheritance xx21.2.4 Wrapping enums

Idea and motivation

[next] [prev] [prev-tail] [tail] [up] 24.1 Idea and motivation The idea behind porting FieldConvert utility to Python is to allow the user to execute a seamless workflow, from converting the mesh and running calculations to preparing data visualisations. This solution is a potential improvement over the original FieldConvert tool, as the Python-based workflow can potentially be executed from a single Python script. Another advantage of the Python version of the tool is the possibility of interacting with other software featuring Python interface, e.g. Paraview. This could make the process of visualising the results of computations done with Nektar++ even easier.

Design and implementation

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

24.2 Design and implementation

FieldConvert utility in Python was designed to provide the user with an easy workflow. Hence, a minimum effort is required to set up the conversion and the majority of work is done behind the scenes.

This section will discuss the design of the utility, including the description of work yet to be done, as well as the user workflow required to execute the basic conversions.

24.2.1 FieldConverter class

The entire procedure outlined in the original FieldConvert utility is managed by the

FieldConverter class located in utilities/FieldConvert/Python/FieldConvert.py. The

class has the following attributes:

fieldSharedPtr: a shared pointer to the field, necessary to initialise a module;

TO-DO: FieldSharedPtr type remains to be wrapped.

TO-DO: Whether a separate FieldSharedPtr is needed for each module initialisation needs to be discussed.

moduleFactory: would be equivalent to GetModuleFactory() in the original utility; TO-DO: Whether this is possible needs to be discussed. TO-DO: It would be best to wrap the general template NekFactory and initialise a module factory this way.

availableModuleList: holds a list of available modules in order to assert that the modules requested by the user are available;

TO-DO: PrintAvailableClasses method needs to be wrapped. TO-DO: It would definitely be neater if the available modules existed as an enum rather than just strings containing names.

sessionFile, inputFile, outputFile: hold the necessary filenames;

moduleList: holds a list of Module objects, created as requested by the user; variableMap: equivalent of vm variable from the original utility.

TO-DO: This variable map needs to be constructed from user input as it needs to be passed into Process method of Module class.

24.2.2 User workflow The currently suggested user workflow is as follows:

Initialise an instance of FieldConverter class. Add a session file as well as an input and an output file.

The converter will assert that the session file and the input file exist, assert that the file extensions are supported and create appropriate modules.

Add any modules, e.g. vorticity.

Currently the argument passed into addProcessModule is a string containing the module name.

TO-DO: As mentioned before, it would be neater if said argument was an enum.

TO-DO: Some more thought is required about how to support modules requiring or permitting parameters. One solution would be to pass a tuple containing module name and parameters. Care needs to be taking in asserting the validity of parameters.

As before, the converter will assert that the module exists and create an appropriate Module class object.

Run the conversion.

The code showcasing the suggested workflow can be found in the main function of the

FieldConvert.py file.

24.2.3 Conversion process Further development and improvement

[next] [prev] [prev-tail] [tail] [up] 24.3 Further development and improvement [next] [prev] [prev-tail] [front] [up] [next] [prev] [prev-tail] [front] [up]

xx

Bibliography

[prev] [prev-tail] [tail] [up] Bibliography

[1] xxxstd::shared_ptr - cppreference.com. [Accessed 21 March 2018].

[2] xxxboost.python/howto - python wiki, 2015. [Accessed 1st May 2018].

[3] xxxosgboostpython/manualwrappingrationale.md at wiki · skylark13/osgboostpython · github, 2015. [Accessed 7th May 2018].

[4] xxxosgboostpython/wrappingcookbook.md at wiki · skylark13/osgboostpython · github, 2015. [Accessed 7th May 2018].

[5] xxxGithub - tng/boost-python-examples: Some examples for the use of boost::python, 2016. [Accessed 7th May 2018].

[6] xxxMark Ainsworth. Pyramid algorithms for bernstein-bxzier finite elements of high, nonuniform order in any dimension. SIAM Journal of Scientific Computing, 36:A543–A569, 2014.

[7] xxxMark Ainsworth, Gaelle Andriamaro, and Oleg Davydov. Bernstein-bxzier finite elements of arbitrary order and optimal assembly procedures. SIAM Journal of Scientific Computing, 33:3087–3109, 2011.

[8] xxxR.xAris. Vectors, tensors, and the basic equations of fluid mechanics. Dover Pubns, 1989.

[9] xxxIvo Babuska, BarnaxA Szabo, and IxNorman Katz. The p-version of the finite element method. SIAM journal on numerical analysis, 18(3):515–545, 1981.

[10]xxxWolfgang Bangerth, Ralf Hartmann, and Guido Kanschat. deal.II–a general-purpose object-oriented finite element library. ACM Transactions on Mathematical Software (TOMS), 33(4):24, 2007.

[11]xxxHughxM Blackburn and SJxSherwin. Formulation of a galerkin spectral element–fourier method for three-dimensional incompressible flows in cylindrical geometries. Journal of Computational Physics, 197(2):759–778, 2004.

[12]xxxA.xBolis, C.D. Cantwell, R.M. Kirby, and S.J. Sherwin. h-to-p efficiently: Optimal implementation strategies for explicit time-dependent problems using the spectral/hp element method. International Journal for Numerical Methods in Fluids, 75:591–607, 2014.

[13]xxxA.I. Borisenko, I.E. Tarapov, and R.A.xSilverman (Translator). Vector and Tensor Analysis with Applications. Dover Books on Mathematics, 2012.

[14]xxxJ.xC. Butcher. General linear methods. Acta Numerica, 15:157–256, 2006.

[15]xxxC.D. Cantwell, D.xMoxey, A.xComerford, A.xBolis, G.xRocco, G.xMengaldo, D.xdexGrazia, S.xYakovlev, J-E Lombard, D.xEkelschot, B.xJordi, H.xXu, Y.xMohamied, C.xEskilsson, B.xNelson, P.xVos, C.xBiotto, R.M. Kirby, and S.J. Sherwin. Nektar++: An open-source spectral/hp element framework. Computer Physics Communications, 192:205–219, 2015.

[16]xxxC.D. Cantwell, S.J. Sherwin, R.M. Kirby, and P.H. Kelly. From h-to-p efficiently: Selecting the optimal spectral/hp discretisation in three dimensions. Math. Model. Nat. Phenom., 6(3):84–96, 2011.

[17]xxxC.D. Cantwell, S.J. Sherwin, R.M. Kirby, and P.H.J. Kelly. From h-to-p efficiently: Strategy selection for operator evaluation on hexahedral and tetrahedral elements. Computers and Fluids, 43:23–28, 2011.

[18]xxxC.D. Cantwell, S.xYakovlev, R.M. Kirby, N.S. Peters, and S.J. Sherwin. High-order continuous spectral/hp element discretisation for reaction-diffusion problems on a surface. Journal of Computational Physics, 257:813–829, 2014.

[19]xxxC.xCanuto, M.Y. Hussaini, A.xQuarteroni, and T.A. Zang. Spectral Methods in Fluid Mechanics. Springer-Verlag, New York, 1987.

[20]xxxBernardo Cockburn, George Karniadakis, and Chi-Wang Shu. Discontinuous Galerkin Methods: Theory, Computation and Applications. Springer-Verlag, 2000.

[21]xxxJ.xAustin Cottrell, Thomas J.xR. Hughes, and Yuri Bazilevs. Isogeometric Analysis: Toward Integration of CAD and FEA. John Wiley and Sons, 2009.

[22]xxxAbrahamsxD. dexGuzmanxJ. Boost.python tutorial - 1.67.0, 2018. [Accessed 1st May 2018].

[23]xxxAndreas Dedner, Robert Klxfkorn, Martin Nolte, and Mario Ohlberger. A generic interface for parallel and adaptive discretization schemes: abstraction principles and the DUNE-FEM module. Computing, 90(3-4):165–196, 2010.

[24]xxxJamesxW. Demmel. Applied Numerical Linear Algebra. SIAM, Philadelphia, PA, USA, 1997.

[25]xxxM.O. Deville, P.F. Fisher, and E.H. Mund. High-Order Methods for Incompressible Fluid Flow. Cambridge University Press, 2002.

[26]xxxJulia Docampo-Sxnchez, JenniferxK Ryan, Mahsa Mirzargar, and RobertxM Kirby. Multi-dimensional filtering: Reducing the dimension through rotation. SIAM Journal on Scientific Computing, 39(5):A2179–A2200, 2017.

[27]xxxM.xDubiner. Spectral methods on triangles and other domains. J. Sci. Comp., 6:345, 1991.

[28]xxxM.G. Duffy. Quadrature over a pyramid or cube of integrands with a singularity at a vertex. SIAM J. Numer. Anal., 19:1260, 1982.

[29]xxxPaul Fischer, James Lottes, Stefan Kerkemeier, Oana Marin, Katherine Heisey, Aleks Obabko, Elia Merzari, and Yulia Peet. Nek5000 User Manual. ANL/MCS-TM-351, 2014.

[30]xxxPythonxSoftware Foundation. Memory management - python 2.7.14 documentation. [Accessed 21 March 2018].

[31]xxxPythonxSoftware Foundation. Reference counting - python 2.7.14 documentation. [Accessed 21 March 2018].

[32]xxxD.xFunaro. Polynomial Approximations of Differential Equations: Lecture Notes in Physics, Volume 8. Springer-Verlag, New York, 1992.

[33]xxxF.X. Giraldo, J.F. Kelly, and E.M. Constantinescu. Implicit explicit formulations of a three dimensional non-hydrostatic unified model of the atmosphere (NUMA). SIAM Journal of Scientific Computing, 35:1162–1194, 2013.

[34]xxxvan RossumxG. GoodgerxD. Pep 257 – docstring conventions, 2001. [Accessed 16th May 2018].

[35]xxxMichaelxT. Heath. Scientific Computing: An Introductory Survey. McGraw-Hill Companies, 2002.

[36]xxxJan Hesthaven, Sigal Gottlieb, and David Gottlieb. Spectral Methods for Time-Dependent Problems. Cambridge University Press, 2007.

[37]xxxJanxS Hesthaven and Tim Warburton. Nodal discontinuous Galerkin methods: algorithms, analysis, and applications, volumex54. Springer, 2007.

[38]xxxJ.S. Hesthaven. From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex. SIAM J. Numer. Anal., 35(2):655–676, 1998.

[39]xxxJ.S. Hesthaven and T.C. Warburton. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Springer Texts in Applied Mathematics 54. Springer Verlag: New York, 2008.

[40]xxxT.xJ.xR. Hughes. The Finite Element Method. Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1987.

[41]xxxAshok Jallepalli, Julia Docampo-Sxnchez, JenniferxK Ryan, Robert Haimes, and RobertxM Kirby. On the treatment of field quantities and elemental continuity in FEM solutions. IEEE Transactions on Visualization and Computer Graphics, 24(1):903–912, 2017.

[42]xxxCem Kaner, Jack Falk, and HungxQuoc Nguyen. Testing Computer Software. John Wiley & Sons, 2010.

[43]xxxGeorgexEm Karniadakis and RobertxM. Kirby. Parallel Scientific Computing in C++ and MPI. Cambridge University Press, New-York, NY, USA, 2003.

[44]xxxGeorgexEm Karniadakis and SpencerxJ. Sherwin. Spectral/hp element methods for Computational Fluid Dynamics (Second Edition). Oxford University Press, 2005.

[45]xxxRobertxM. Kirby and SpencerxJ. Sherwin. Aliasing errors due to quadratic non-linearities on triangular spectral/hp element discretisations. Journal of Engineering Mathematics, 56:273–288, 2006.

[46]xxxDavidxA. Kopriva. Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engineers. Springer, 2009.

[47]xxxAnders Logg, Kent-Andre Mardal, and GarthxWells (editors). Automated Solution of Differential Equations by the Finite Element Method. Springer Lecture Notes in Computational Science and Engineering, Volume 84, 2012.

[48]xxxDaconta M. C++ pointers and dynamic memory management. New York: Wiley, 1995.

[49]xxxReddy M. API Design for C++. Burlington: Elsevier, 2011.

[50]xxxA.T.T. McRae, G.-T. Bercea, L.xMitchell, D.A. Ham, and C.J. Cotter. Automated generation and symbolic manipulation of tensor product finite elements. SIAM Journal on Scientific Computing, 38(5):S25–S47, 2016.

[51]xxxScott Meyers. Effective C++: 55 Specific Ways to Improve Your Programs and Designs (Third Edition). Addison-Wesley Professional, 2005.

[52]xxxJaroszyxski P. Piotr jaroszyxski’s blog: Boost.python: docstrings in enums, 2007. [Accessed 16th May 2018].

[53]xxxJhong E. etxal. PatelxA., PicardxA. Google python style guide. [Accessed 16th May 2018].

[54]xxxCh. Schwab. p– and hp– Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics. Oxford University Press, 1999.

[55]xxxBjarne Stroustrup. The C++ Programming Language (Fourth Edition). Addison-Wesley Professional, 2013.

[56]xxxM.xTaylor and B.A. Wingate. The fekete collocation points for triangular spectral elements. Journal on Numerical Analysis, 1998.

[57]xxxM.xTaylor, B.A. Wingate, and R.E. Vincent. An algorithm for computing Fekete points in the triangle. SIAM J. Num. Anal., 38:1707–1720, 2000.

[58]xxxLloydxN. Trefethen. Is gauss quadrature better than clenshaw-curtis? SIAM Review, 50:67–87, 2008.

[59]xxxLloydxN. Trefethen and III DavidxBau. Numerical Linear Algebra. SIAM, Philadelphia, PA, USA, 1997.

[60]xxxTomxx Vejchodskx, Pavel xolxn, and Martin Zxtka. Modular hp-FEM system HERMES and its application to Maxwell’s equations. Mathematics and Computers in Simulation, 76(1):223–228, 2007.

[61]xxxPeter E.xJ. Vos, SpencerxJ. Sherwin, and RobertxM. Kirby. h-to-p efficiently: Implementing finite and spectral/hp element methods to achieve optimal performance for low- and high-order discretisations. Journal of Computational Physics, 229:5161–5181, 2010.

[62]xxxPeterxE.J. Vos, Sehun Chun, Alessandro Bolis, Claes Eskilsson, RobertxM. Kirby, and SpencerxJ. Sherwin. A generic framework for time-stepping pdes: General linear methods, object-oriented implementations and applications to fluid problems. International Journal of Computational Fluid Dynamics, 25:107–125, 2011.

[63]xxxF.D. Witherden, P.E. Vincent, and A.xJameson. Chapter 10 – high-order flux reconstruction schemes. Handbook of Numerical Analysis, 17:227–263, 2016.

[64]xxxFRxWitherden, AMxFarrington, and PExVincent. PyFR: An open source framework for solving advection-diffusion type problems on streaming architectures using the flux reconstruction approach. Computer Physics Communications, 185:3028–3040, 2014.

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