Welcome to Astrix!

Astrix (AStrophysical fluid dynamics on TRIangular eXtreme grids) is a CUDA/C++ implementation of a two-dimensional residual distribution scheme aimed at tackling problems in astrophysical fluid dynamics.

Contents:

Changelog

This changelog includes the most important changes in recent updates.

Version 1.1

  • Choose conservation law from command line
  • Source term integration
  • Contour integration for general conservation laws
  • Support for isothermal hydrodynamics
  • Symmetry boundaries
  • Limit LDA contribution if change too much
  • Various bug fixes in Mesh generation

Version 1.0.1

  • Bugfix in documentation generation

Installation

Astrix is completely open-source and freely available from GitHub. It runs on both Mac and Linux machines (Windows is not supported) as long as CUDA is installed.

Installing CUDA

Note that while Astrix does not require a GPU to run, it needs the CUDA compiler in order to be built. CUDA can be obtained for free from https://developer.nvidia.com/cuda-downloads. There are three parts that can be installed:

  • CUDA toolkit
  • CUDA samples
  • CUDA driver

Only the CUDA toolkit is required for building Astrix. Installing the driver is only necessary on machines that you plan to use the GPU of (requires root access); installing the samples is optional but they include some useful tests to see if everything is working correctly. Unfortunately for Mac users, CUDA only works with specific versions of Xcode, see http://docs.nvidia.com/cuda/cuda-installation-guide-mac-os-x/index.html#prerequisites. For supported versions of Linux see http://docs.nvidia.com/cuda/cuda-installation-guide-linux/index.html#system-requirements.

Installing Astrix

Once CUDA has been installed correctly, it is time to download Astrix. At a terminal, type:

git clone http://github.com/SijmeJan/Astrix

This will create a directory Astrix. Now do:

cd Astrix

Astrix needs to know where CUDA was installed. The easiest way is to make sure that the CUDA bin directory (for example /usr/local/cuda/bin) is in your PATH environmental variable. You can check if this is the case by typing which nvcc. If this gives a result, you are good to go. If not, add the CUDA bin directory to your path, and then type:

make astrix

Sit back and relax because this may take a while. When finished, Astrix is ready to use. The executable will be located in Astrix/bin.

Advanced options

The compilation process can be sped up is you know what GPU you are going to use. If you are planning to run on a Tesla K20m, which has CUDA compute capability 3.5, issue:

make astrix CUDA_COMPUTE=35

Not only will this speed up compilation, it will also speed up the code itself since the compiler can use a few tricks not available on devices with compute capability < 3.5. The minimum compute capability currently supported by nVidia is 3.0; while Astrix should work on devices with compute capability 2.x this is not supported by the standard build.

It is possible to override the standard CUDA installation and build with a specific version located in path/to/cuda:

make astrix CUDA_INSTALL_PATH=path/to/cuda

This can be useful when testing a new CUDA version.

By default, Astrix will perform floating point computations in single precision. It is possible to force double presicion through:

make astrix ASTRIX_DOUBLE=1

Note that switching between double and single precision requires a complete rebuild: therefore first remove a previous build by entering:

make clean

A simple visualisation program is included and can be built by:

make visAstrix

This requires the OpenGL and glut libraries to be installed.

If you are interested in building this documentation locally, it can be built through:

make doc

In order to clean up everything, enter:

make allclean

which removes the Astrix build, any visAstrix build and all documentation.

Running Astrix

Command line options

Issueing astrix gives:

Usage: astrix [-d] [-v verboseLevel] [-D debugLevel] [-r restartNumber] [-cl conservationLaw] filename
-d                  : run on GPU device
-v verboseLevel     : amount of output to stdout (0 - 2)
-D debugLevel       : amount of extra checks for debugging
-r restartNumber    : try to restart from previous dump
-cl conservationLaw : use different conservation law. Can be
                      either "advect", "burgers"
                      "cart_iso" or "cart_euler"
filename            : input file name

Test problems

A few test problems for the Euler equations are supplied in the Astrix/run/euler directory. Each directory contains an input file astrix.in specifying parameters of the simulation. See the classes MeshParameter and SimulationParameter for details. Note that all parameters have to be present in the input file. For each test problem, run Astrix in the respective directory:

  • blast/ : A one-dimensional problem of two interacting blast waves.
  • cyl/ : supersonic flow around a cylinder.
  • kh/ : Kelvin-Helmholtz instability.
  • linear/ : A one-dimensional problem of a linear sound wave.
  • noh/ : The Noh test problem.
  • riemann/ : Two-dimensional Riemann problem.
  • sod/ : A one-dimensional shock tube.
  • source/ : Rayleigh-Taylor instability.
  • vortex/ : Isentropic vortex test.

In addition, some scalar equation tests can be found in Astrix/run/scalar. A suite of test problems can be run by entering, in the Astrix directory:

python python/astrix/testsuite.py ./

which will generate a pdf document with outputs from most test problems.

Creating a new setup

Going beyond the test problems supplied requires a bit of coding. It is instructive to look at places in the code that depend on the ProblemDefinition specifying which test is being run. For example, in the Astrix directory enter:

grep -r PROBLEM_CYL src/astrix/*

to see a list of places where code specific to the test problem of the flow around a cylinder is executed. These will appear as part of the Mesh class, setting up the initial mesh with the inner hole, as well as in the Simulation class, most notably to set initial and boundary conditions. Additional setups can be added by adding an entry to the ProblemDefinition enum in definitions.h, and create specific mesh generation instructions, initial and boundary conditions for this new setup.

Output and Visualisation

Simulation output comes in three kinds:

  • Raw data of both Mesh and Simulation. Every save interval, both the Mesh and the state are written to disc. Mesh information is written in three files: vert####.dat, containing vertex coordinates, tria####.dat, containing triangle information (vertices and edges), and edge####.dat, containing edge information (triangles). Here and in the following, #### stands for a four-digit number, e.g. 0001, 0199. The state vector is written in four files dens####.dat, containing the density, momx####.dat containing the x-momentum, momy####.dat containing the y-momentum and ener####.dat containing the total energy.
  • Fine grain global data in simulation.dat. Every fine grain save interval, a new line is added to this ASCII file, containing global simulation quantities (simulation time plus other quantities that might be interesting to monitor).
  • When desired, Astrix can output legacy VTK files for easy visualisation for example with the open source package VisIt (available from https://wci.llnl.gov/simulation/computer-codes/visit)

All raw data files are binary files. The format for each is as follows:

  • vert####.dat: three 4 byte integers specifying the number of spatial dimensions (fixed to two for now), the size of each element (can be 4 byte float or 8 byte double), and the total number of vertices Nv. These are followed by the actual data: Nv x coordinates followed by Nv y coordinates.
  • tria####.dat: a 4 byte integer specifying the total number of triangles Nt. This is followed by Nt integers specifying the first vertex for each triangle, followed by Nt integers specifying the second vertex for each triangle, followed by Nt integers specifying the third vertex for each triangle. After this, the edges for every triangle: three sets of Nt integers specifying the first, second and third edge for every triangle.
  • edge####.dat: a 4 byte integer specifying the total number of edges Ne. This is followed by Ne integers specifying the first neighbouring triangle, followed by Ne integers specifying the second neighbouring triangle. These entries can be -1 if the edge happens to have only one neighbouring triangle.
  • dens####.dat (and similar for the other state vector files): a 4 byte integer specifying the size of the output (can be float or double), a float or double specifying the current simulation time, followed by an 4 byte integer specifying the number of time steps taken so far. Then the actual data follows for each vertex a double or a float.

Constructing periodic meshes from the raw outputs can be tricky. A python script to read raw data files into a format that can be used for example with matplotlib is provided in the Astrix/python/astrix/ directory. The module readfiles contains a function readall that reads in both Mesh and Simulation data. Provided Astrix/python/ is in your PYTHONPATH environmental variable, this function can be used for plotting as follows:

import numpy as np
import matplotlib.pyplot as plt
import astrix.readfiles as a

# Read in data
coords, triang, state = a.readall("path/to/data/", 0)

# Vertex coordinates
x = coords[:, 0]
y = coords[:, 1]

# Density
d = state[:, 0]

# Contour levels
levels = np.linspace(np.min(d), np.max(d), 100)

fig = plt.figure(figsize=(5.5, 5.5))
ax = fig.add_subplot(1, 1, 1)

# Plot triangulation
ax.triplot(x, y, triang, 'k-')
# Contour plot density
ax.tricontourf(x, y, triang, d, levels=levels, cmap=cm.bwr)

plt.show()

The complete content of the astrix.readfiles module is given below. In most cases, the readall function is all that is required.

astrix.readfiles.CanVertexBeTranslatedX(a, N)[source]

Check whether a vertex is a periodic x vertex.

For a mesh that is periodic in x, some triangles ‘wrap around’ and therefore need to look at the other side of the mesh for the coordinates of one of their vertices. This function checks if this is the case for vertex a.

Parameters:
  • a (int) – Vertex input to consider.
  • N (int) – Total number of vertices in Mesh
Returns:

1 if a can be found one period away towards positive x, -1 if a can be found one period away towards negative x, and zero otherwise.
Return type:

int
astrix.readfiles.CanVertexBeTranslatedY(a, N)[source]

Check whether a vertex is a periodic y vertex.

For a mesh that is periodic in y, some triangles ‘wrap around’ and therefore need to look at the other side of the mesh for the coordinates of one of their vertices. This function checks if this is the case for vertex a.

Parameters:
  • a (int) – Vertex input to consider.
  • N (int) – Total number of vertices in Mesh
Returns:

1 if a can be found one period away towards positive y, -1 if a can be found one period away towards negative y, and zero otherwise.
Return type:

int
astrix.readfiles.GetCoordinates(a, vertX, vertY, Px, Py)[source]

Get coordinates of vertex number a.

For a mesh that is periodic in x, some triangles ‘wrap around’ and therefore need to look at the other side of the mesh for the coordinates of one of their vertices. This function gets the ‘proper’ coordinates of vertex a.

Parameters:
  • a (int) – Vertex input to consider.
  • vertX (ndarray) – array of vertex x coordinates
  • vertY (ndarray) – array of vertex y coordinates
  • Px (float) – Period in x
  • Py (float) – Period in y
Returns:

x and y coordinate of vertex a.
Return type:

float, float
astrix.readfiles.readState(read_direc, read_indx, nVertex)[source]

Reads in simulation data from Astrix simulation at specified time.

Read the simulation output (i.e. the state at every vertex) contained in directory read_direc at snapshot read_indx.

Parameters:
  • read_direc (string) – Directory containing Astrix output.
  • read_indx (int) – Output number to read.
  • nVertex (int) – total number of vertices in Mesh
Returns:

four arrays containing density, x-velocity, y-velocity and total energy
Return type:

ndarray, ndarray, ndarray, ndarray
astrix.readfiles.readTriangle(read_direc, read_indx)[source]

Reads in triangle data from Astrix simulation at specified time.

Read the triangle output (i.e. the vertices belonging to every triangle) contained in directory read_direc at snapshot read_indx. Note that in case of periodic meshes, the entries may be smaller than zero or larger than the number of vertices.

Parameters:
  • read_direc (string) – Directory containing Astrix output.
  • read_indx (int) – Output number to read.
Returns:

a connectivity array of length 3 times the number of triangles. The first three entries represent the first triangle, etc.
Return type:

ndarray
astrix.readfiles.readVertex(read_direc, read_indx)[source]

Reads in vertex data from Astrix simulation at specified time.

Read the vertex output contained in directory read_direc at snapshot read_indx.

Parameters:
  • read_direc (string) – Directory containing Astrix output.
  • read_indx (int) – Output number to read.
Returns:

x-coordinates of the vertices; y-coordinates of the vertices
Return type:

ndarray, ndarray
astrix.readfiles.readall(read_direc, read_indx, Px=1.0, Py=1.0)[source]

Reads in data from Astrix simulation at specified time.

Read the simulation output contained in directory read_direc at snapshot read_indx. If the mesh is periodic in x or y or both, the periods must be supplied as Px and Py so that we can create the full mesh.

Parameters:
  • read_direc (string) – Directory containing Astrix output.
  • read_indx (int) – Output number to read.
  • Px (float) – Optional distance in x over which the mesh is periodic.
  • Py (float) – Optional distance in y over which the mesh is periodic.
Returns:

Coordinates (x,y) of the vertices as a (Nv, 2) array, where Nv is the number of vertices; triangulation as a (Nt, 3) array, where Nt is the number of triangles; state as a (Nv, 4) array, containing density, two velocities and the total energy.
Return type:

ndarray(Nv,2), ndarray(Nt,3), ndarray(Nv,4)

Class overview

Simulation

template <class realNeq, ConservationLaw CL>
class

Simulation: class containing simulation.

This is the basic class needed to run an Astrix simulation.

Public Functions

template astrix::Simulation< realNeq, CL >::Simulation(int _verboseLevel, int _debugLevel, char * fileName, Device * _device, int restartNumber)

Constructor for Simulation object.

Define Arrays, create Mesh and setup simulation.

Parameters
  • _verboseLevel -

    How much information to output to stdout in Astrix.

  • _debugLevel -

    Level of extra checks for correct mesh.

  • *fileName -

    Input file name

  • *device -

    Device to be used for computation.

  • restartNumber -

    Number of saved file to restore from

template astrix::Simulation< realNeq, CL >::~Simulation()

Destructor, releases all dynamically allocated memory.

template void astrix::Simulation< realNeq, CL >::Run(real maxWallClockHours)

Run simulation.

Run simulation, possibly restarting from saved state, for a maximum amount of wall clock hours.

Parameters
  • maxWallClockHours -

    Maximum amount of wall clock hours to run for.

SimulationParameter

class

Class containing parameters for the Simulation.

Various parameters governing the Simulation are read from an input file and stored in this class

Public Functions

astrix::SimulationParameter::SimulationParameter()

Constructor, set all data to invalid values.

astrix::SimulationParameter::~SimulationParameter()

Destructor.

void astrix::SimulationParameter::ReadFromFile(const char *fileName, ConservationLaw CL)

Read in data from file.

Read in data from input file. File is read line by line; input parameters can appear in any order but all must be present. An exception is thrown if any parameter is missing or has an invalid value.

Parameters
  • *fileName -

    Pointer to input file name

Public Members

ProblemDefinition astrix::SimulationParameter::problemDef

Problem specification (see Common/definitions.h)

Converted to int using definitions in definitions.h.

int astrix::SimulationParameter::nSpaceDim

Number of space dimensions (fixed to 2)

IntegrationScheme astrix::SimulationParameter::intScheme

Integration scheme (see Common/definitions.h)

Converted to int using definitions in definitions.h.

int astrix::SimulationParameter::integrationOrder

Order of accuracy in time (1 or 2)

int astrix::SimulationParameter::massMatrix

Mass matrix formulation to use (1, 2, 3 or 4)

int astrix::SimulationParameter::selectiveLumpFlag

Flag whether to use selective lumping.

real astrix::SimulationParameter::CFLnumber

Courant number.

int astrix::SimulationParameter::preferMinMaxBlend

Preference for using minimum/maximum value of blend parameter.

real astrix::SimulationParameter::specificHeatRatio

Ratio of specific heats.

real astrix::SimulationParameter::maxSimulationTime

Maximum simulation time.

real astrix::SimulationParameter::saveIntervalTime

Time between 2D saves.

real astrix::SimulationParameter::saveIntervalTimeFine

Time between 0D saves.

int astrix::SimulationParameter::writeVTK

Flag whether do output VTK files.

Device

class

Simple class containing information about device.

This class is used to hold some very basic information about the machine the simulation is run on: whether we want to use any CUDA capable device, and how many CUDA-capable devices there are in total.

Public Functions

astrix::Device::Device(int _cudaFlag)

Constructor.

Construct Device object. Count number of CUDA-capable devices and display capabilities on screen. By default, device 0 is used.

Parameters
  • _cudaFlag -

    Flag whether to run on CUDA device. If set to zero, still count CUDA devices but do not use them to for computation.

astrix::Device::~Device()

Destructor.

Free Device object. If using CUDA, reset device for clean exit.

int astrix::Device::GetDeviceCount()

Return number of CUDA devices.

int astrix::Device::GetCudaFlag()

Return flag whether using CUDA.

Mesh

class

Class containing functions acting on the whole mesh.

A Mesh object can be thought of representing the Delaunay mesh on which the computations are performed.

Public Functions

astrix::Mesh::Mesh(int meshVerboseLevel, int meshDebugLevel, int meshCudaFlag, const char *fileName, Device *device, int restartNumber)

Constructor.

Create Mesh with parameters specified in file fileName.

Parameters
  • meshVerboseLevel -

    Level of output to stdout

  • meshDebugLevel -

    Level of debugging in Mesh construction and maintenance

  • meshCudaFlag -

    Flag whether to use CUDA device

  • *fileName -

    Input file name

  • *device -

    Pointer to Device class containing information about any CUDA device present

  • restartNumber -

    Number of save file to restore from

astrix::Mesh::~Mesh()

Destructor; releases memory.

template <class realNeq, ConservationLaw CL>
template int astrix::Mesh::ImproveQuality< real4, CL_CART_EULER >(Array < realNeq > * vertexState, real specificHeatRatio, int nTimeStep)

Refine mesh.

Refine Delaunay mesh, depending on the state vector (i.e. density, momentum, etc). If *vertexState != 0, we first calculate an estimate of the discretization error and flag triangles to be refined; otherwise, all triangles can be refined, necessary for example when building the mesh for the first time. It returns the number of vertices added.

Parameters
  • *vertexState -

    Pointer to state vector at vertices

  • specificHeatRatio -

    Ratio of specific heats

  • nTimeStep -

    Number of time steps taken so far. Used in combination with nStepSkipRefine to possibly avoid refining every timestep

template <class realNeq, ConservationLaw CL>
template int astrix::Mesh::RemoveVertices< real4, CL_CART_EULER >(Array < realNeq > * vertexState, real specificHeatRatio, int nTimeStep)

Coarsen mesh.

Coarsen mesh. First calculate an estimate of the discretization error and flag triangles that can be coarsened. Then remove as many vertices as possible from mesh.

Parameters
  • *vertexState -

    Pointer to state vector at vertices

  • specificHeatRatio -

    Ratio of specific heats

  • nTimeStep -

    Number of time steps taken so far. Used in combination with nStepSkipCoarsen to possibly avoid coarsening every timestep

void astrix::Mesh::Save(int nSave)

Save mesh to disk.

Save vertexCoordinates in a vertex file, triangleVertices and triangleEdges in a triangle file, and edgeTriangles in an edge file.

Parameters
  • nSave -

    Number of save, used to generate file names

void astrix::Mesh::ReadFromDisk(int nSave)

Read previously created mesh from disk.

Read mesh from disk as it was saved under number nSave. In addition, calculate triangle normals, vertex areas and find boundary vertices.

Parameters
  • nSave -

    Number of save to restore

int astrix::Mesh::GetNVertex()

Return number of vertices.

int astrix::Mesh::GetNTriangle()

Return number of triangles.

int astrix::Mesh::GetNEdge()

Return number of edges.

int astrix::Mesh::IsAdaptive()

Return if mesh is adaptive.

real astrix::Mesh::GetPx()

Return size of domain in x.

real astrix::Mesh::GetPy()

Return size of domain in y.

real astrix::Mesh::GetMinX()

Return minimum x.

real astrix::Mesh::GetMaxX()

Return maximum x.

real astrix::Mesh::GetMinY()

Return minimum y.

real astrix::Mesh::GetMaxY()

Return maximum y.

real astrix::Mesh::GetTotalArea()

Return total vertex area.

const real2 *astrix::Mesh::VertexCoordinatesData()

Return pointer to vertex coordinates data.

const int *astrix::Mesh::VertexBoundaryFlagData()

Return pointer to vertex boundary flag data.

const real *astrix::Mesh::VertexAreaData()

Return pointer to vertex area data.

const int3 *astrix::Mesh::TriangleVerticesData()

Return pointer to triangle vertices data.

const int3 *astrix::Mesh::TriangleEdgesData()

Return pointer to triangle edges data.

const real2 *astrix::Mesh::TriangleEdgeNormalsData(int dim)

Return pointer to triangle edge normals data.

const real3 *astrix::Mesh::TriangleEdgeLengthData()

Return triangle edge length data.

const int2 *astrix::Mesh::EdgeTrianglesData()

Return edge triangles data.

void astrix::Mesh::Transform()

Transform all Arrays.

MeshParameter

class

Class containing parameters for the Mesh.

Various parameters governing the resolution and quality of the Mesh are read from an input file and stored in this class

Public Functions

astrix::MeshParameter::MeshParameter()

Constructor, set all data to invalid values.

astrix::MeshParameter::~MeshParameter()

Destructor.

void astrix::MeshParameter::ReadFromFile(const char *fileName)

Read in data from file.

Read in data from input file. File is read line by line; input parameters can appear in any order but all must be present. An exception is thrown if any parameter is missing or has an invalid value.

Parameters
  • *fileName -

    Pointer to input file name

Public Members

ProblemDefinition astrix::MeshParameter::problemDef

Definition of test problem.

int astrix::MeshParameter::equivalentPointsX

Approximate number of vertices in x direction.

real astrix::MeshParameter::qualityBound

Quality bound on triangles; must be >= 1.

int astrix::MeshParameter::periodicFlagX

Flag to create periodic domain in x.

int astrix::MeshParameter::periodicFlagY

Flag to create periodic domain in y.

real astrix::MeshParameter::minx

Position of left x boundary.

real astrix::MeshParameter::maxx

Position of right x boundary.

real astrix::MeshParameter::miny

Position of left y boundary.

real astrix::MeshParameter::maxy

Position of right y boundary.

int astrix::MeshParameter::structuredFlag

Flag whether using structured mesh.

int astrix::MeshParameter::adaptiveMeshFlag

Flag whether mesh is adaptive.

int astrix::MeshParameter::maxRefineFactor

Maximum factor to increase resolution over base resolution if adaptive mesh is used.

int astrix::MeshParameter::nStepSkipRefine

Number of time steps without checking if refinement is needed.

int astrix::MeshParameter::nStepSkipCoarsen

Number of time steps without checking if coarsening is needed.

real astrix::MeshParameter::minError

If discretization error smaller than minError, coarsen Mesh.

real astrix::MeshParameter::maxError

If discretization error larger than maxError, refine Mesh.

real astrix::MeshParameter::baseResolution

Triangle size for initial Mesh (derived from equivalentPointsX)

real astrix::MeshParameter::maxResolution

Triangle size for adaptive mesh (derived from baseResolution and maxRefineFactor)

Connectivity

class

Class containing Mesh data structure.

Class containing coordinates and connectivity of Mesh; data needed by all Mesh-related classes. The class is essentially data-only, plus a few functions to move data between host and device. All data members are public, which can be unsafe. It is assumed that at the end of any function modifying the Mesh, the Connectivity represents a valid triangulation (not necessarily Delaunay), and that the sizes of the Arrays are properly set, i.e. the size of vertexCoordinates is the number of vertices, the size of both triangleVertices and triangleEdges equals the number of triangles, and the size of edgeTriangles equals the number of edges.

Public Functions

astrix::Connectivity::Connectivity(int _cudaFlag)

Constructor.

Constructor for Connectivity class. Memory is allocated in large chunks to minimise any further calls to cudaMalloc when improving the Mesh.

Parameters
  • _cudaFlag -

    Flag whether Arrays reside on host (0) or device (1)

astrix::Connectivity::~Connectivity()

Destructor; releases memory.

void astrix::Connectivity::Transform()

Transform from device to host or vice versa.

Move whole class to device or host, depending on current value of cudaFlag

void astrix::Connectivity::CopyToHost()

Copy data to host.

Copy all data currently residing on device to host. Unlike during a transform, cudaFlag is not changed

void astrix::Connectivity::CopyToDevice()

Copy data to device.

Copy all data currently residing on the host to the device. Unlike during a transform, cudaFlag is not changed

void astrix::Connectivity::CheckEdgeTriangles()

Check if the neighbouring triangles of every edge e do have e as an edge. This is done on the host

void astrix::Connectivity::CalcVertexArea(real Px, real Py)

Calculate area associated with vertices (Voronoi cells)

Every triangle contributes one third of its area to the area of the Voronoi cell associated with its vertices. Atomically add this contribution to vertexArea

Public Members

Array<real2> *astrix::Connectivity::vertexCoordinates

Coordinates of vertices.

Array<int3> *astrix::Connectivity::triangleVertices

Vertices belonging to triangles.

Array<int3> *astrix::Connectivity::triangleEdges

Edges belonging to triangles.

Array<int2> *astrix::Connectivity::edgeTriangles

Triangles associated with edges.

Array<real> *astrix::Connectivity::vertexArea

Vertex area (area of Voronoi cell)

Delaunay

class

Class for turning Mesh into Delaunay mesh.

The public function of this class turns the Mesh into a Delaunay mesh by edge flipping.

Public Functions

astrix::Delaunay::Delaunay(int _cudaFlag, int _debugLevel)

Constructor.

Constructor for Delaunay class.

Parameters
  • _cudaFlag -

    Flag indicating whether to use device (1) or host (0)

  • _debugLevel -

    Level of extra checks to do

astrix::Delaunay::~Delaunay()

Destructor; releases memory.

template <class realNeq, ConservationLaw CL>
template void astrix::Delaunay::MakeDelaunay< real4, CL_CART_EULER >(Connectivity *const connectivity, Array < realNeq > *const vertexState, const Predicates * predicates, const MeshParameter * meshParameter, const int maxCycle, Array < int > *const edgeNeedsChecking, const int nEdgeCheck, const int flopFlag)

Turn Mesh into Delaunay nesh.

Transform triangulated Mesh into Delaunay Mesh. This is achieved by flipping edges that do not have the Delaunay property. First, we make a list of edges that are not Delaunay, then we select those that can be flipped in parallel, we adjust the state vector in order to conserve mass, momentum and energy, and finally we flip the edges. A repair step ensures all edges have the correct neighbouring triangles. This is repeated until all edges are Delaunay.

Parameters
  • *connectivity -

    Pointer to basic Mesh data

  • *vertexState -

    Pointer to state vector

  • *predicates -

    Pointer to Predicates object, used to check Delaunay property without roundoff error

  • *meshParameter -

    Pointer to Mesh parameters

  • maxCycle -

    Maximum number of cycles. If <= 0, cycle until all edges are Delaunay

Refine

class

Public Functions

astrix::Refine::Refine(int _cudaFlag, int _debugLevel, int _verboseLevel)

Constructor.

Constructor for Refine class.

Parameters
  • _cudaFlag -

    Flag indicating whether to use device (1) or host (0)

  • _debugLevel -

    Level of extra checks to do

  • _verboseLevel -

    Level of output to stdout

astrix::Refine::~Refine()

Destructor; releases memory.

template <class realNeq, ConservationLaw CL>
template int astrix::Refine::ImproveQuality< real4, CL_CART_EULER >(Connectivity *const connectivity, const MeshParameter * meshParameter, const Predicates * predicates, Morton *const morton, Delaunay *const delaunay, Array < realNeq > *const vertexState, const real specificHeatRatio, Array < int > *const triangleWantRefine)

Add vertices to Mesh until quality constraints met.

Function to improve quality of Mesh by adding new vertices, until the requirements as specified in MeshParameter are met. Returns the number of vertices that were added.

Parameters
  • *connectivity -

    Pointer to basic Mesh data: vertices, triangles, edges

  • *meshParameter -

    Pointer to Mesh parameters, read from input file

  • *predicates -

    Exact geometric predicates

  • *morton -

    Pointer to Morton object, used for sorting to improve data locality

  • *delaunay -

    Pointer to Delaunay object, used to maintain Delaunay triangulation

  • *vertexState -

    Pointer to state at vertices. If refining at t > 0, we need to interpolate the state at new vertices. Otherwise, this pointer needs to be 0 and will be ignored.

  • specificHeatRatio -

    Ratio of specific heats. Needed when interpolating the state.

  • *triangleWantRefine -

    Pointer to flags if triangle needs to be refined based on current state. Only used when t > 0; otherwise it needs to be 0.

int astrix::Refine::AddVertices(Connectivity *const connectivity, const MeshParameter *meshParameter, const Predicates *predicates, Delaunay *const delaunay, Array<real2> *const vertexBoundaryCoordinates, Array<int> *const vertexOrder)

Add list of vertices to Mesh.

Add a list of vertices to the Mesh. Used in cases where vertex locations are known in advance, i.e. when inserting boundaries. Can only be used at t = 0, so no interpolation of the state is necessary.

Parameters
  • *connectivity -

    Pointer to basic Mesh data: vertices, triangles, edges

  • *meshParameter -

    Pointer to Mesh parameters, read from input file

  • *predicates -

    Exact geometric predicates

  • *delaunay -

    Pointer to Delaunay object, used to maintain Delaunay triangulation

  • *vertexBoundaryCoordinates -

    Pointer to coordinates of vertices to be added

  • *vertexOrder -

    Output: the order in which vertices were inserted.

Morton

class

Class containing functions for Morton ordering.

A Morton object can be used to reorder vertices, triangles and edges to improve data locality.

Public Functions

astrix::Morton::Morton(int _cudaFlag)

Constructor.

Constructor for Morton object. Allocates Arrays of standard size.

Parameters
  • _cudaFlag -

    Flag whether to use device (1) or host (0) to compute

astrix::Morton::~Morton()

Destructor; releases memory.

template <class realNeq, ConservationLaw CL>
template void astrix::Morton::Order< real4, CL_CART_EULER >(Connectivity *const connectivity, Array < int > *const triangleWantRefine, Array < realNeq > *const vertexState)

Reorder Arrays for maximum locality.

Sort vertices, edges and triangles according to their Morton value to improve data locality. We first compute the Morton values for the vertex coordinates and use those to sort vertices, triangles and edges.

Parameters
  • *connectivity -

    Pointer to basic Mesh data

  • *triangleWantRefine -

    Pointer to flags whether triangles need to be refined

  • *vertexState -

    Pointer to state vector

Predicates

class

Class for exact geometric predicates.

Creating and updating Delaunay triangulations requires exact evaluation of certain geometric tests, in particular whether a point d lies inside or outside the circle through three other points a, b and c, and whether three points a, b and c are orientated in anti-clockwise direction or not. This class uses the algorithms by Shewchuck (1997).

Public Functions

astrix::Predicates::Predicates(Device *device)

Constructor.

Construct Predicates object and compute parameter vector on host and device.

Parameters
  • *device -

    Pointer to Device object.

astrix::Predicates::~Predicates()

Destructor.

real astrix::Predicates::incircle(real ax, real ay, real bx, real by, real cx, real cy, real dx, real dy, const real *const pParam)
const

Test whether point d lies in circle formed by a, b and c.

Test whether point (dx, dy) lies inside a circle defined by points (ax, ay), (bx, by) and (cx, cy). Return value > 0 if d lies inside circle, < 0 if outside circle, and = 0 if exactly on circle.

Parameters
  • ax -

    x coordinate of point a

  • ay -

    y coordinate of point a

  • bx -

    x coordinate of point b

  • by -

    y coordinate of point b

  • cx -

    x coordinate of point c

  • cy -

    y coordinate of point c

  • dx -

    x coordinate of point d

  • dy -

    y coordinate of point d

  • *pParam -

    pointer to parameter vector

real astrix::Predicates::orient2d(real ax, real ay, real bx, real by, real cx, real cy, const real *const pParam)
const

Test whether points a, b and c lie in counterclockwise orientation.

Test whether points (ax, ay), (bx, by) and (cx, cy) are orientated in counterclockwise direction. Return value > 0 if a, b, c occur in counterclockwise order, < 0 if in clockwise order, and = 0 if points are collinear.

Parameters
  • ax -

    x coordinate of point a

  • ay -

    y coordinate of point a

  • bx -

    x coordinate of point b

  • by -

    y coordinate of point b

  • cx -

    x coordinate of point c

  • cy -

    y coordinate of point c

  • *pParam -

    pointer to parameter vector

real *astrix::Predicates::GetParamPointer(int cudaFlag)
const

Get pointer to parameter vector.

Return pointer to parameter vector, either device pointer (if cudaFlag = 1) or host pointer (if cudaFlag = 0). Note that both exist unless there is no CUDA capable device.

Parameters
  • cudaFlag -

    Flag whether to return device pointer (= 1) or host pointer (= 0).

Array

template <class T>
class

Basic vector-like class used in Astrix.

An Array object can be thought of as a vector that can either live on the host or on the device.

Public Functions

template astrix::Array< T >::Array()

Basic constructor for 1D Array on host.

Construct 1D array on host of size dynArrayStep.

template astrix::Array< T >::Array(unsigned int _nDims, int _cudaFlag)

Constructor for multidimensional Array, possibly on device.

Construct Array of dimensions _nDims, on device if _cudaFlag = 1, of size dynArrayStep.

Parameters
  • _nDims -

    Number of dimensions

  • _cudaFlag -

    Create Array on host (=0) or on device (=1)

template astrix::Array< T >::Array(unsigned int _nDims, int _cudaFlag, unsigned int _size)

Constructor for multidimensional Array of specified size, possibly on device.

Construct Array of dimensions _nDims, on device if _cudaFlag = 1, of size _size.

Parameters
  • _nDims -

    Number of dimensions

  • _cudaFlag -

    Create Array on host (=0) or on device (=1)

  • _size -

    Size of every dimension

template astrix::Array< T >::Array(unsigned int _nDims, int _cudaFlag, unsigned int _size, int _dynArrayStep)

Constructor for multidimensional Array of specified size and dynArrayStep, possibly on device.

Construct Array of dimensions _nDims, on device if _cudaFlag = 1, of size _size, with dynamical size increase step _dynArrayStep.

Parameters
  • _nDims -

    Number of dimensions

  • _cudaFlag -

    Create Array on host (=0) or on device (=1)

  • _size -

    Size of every dimension

  • _dynArrayStep -

    Increase physical size of array in these steps

template astrix::Array< T >::~Array()

Destructor, releases allocated memory.

Destroy Array object, releasing both host and device memory

template void astrix::Array< T >::TransformToDevice()

Transform from host vector to device vector.

template void astrix::Array< T >::TransformToHost()

Transform from device vector to host vector.

template<>
int astrix::Array<T>::GetCudaFlag()
const

Return whether data currently resides on host or device.

template unsigned int astrix::Array< T >::GetSize() const

Return size of array.

template unsigned int astrix::Array< T >::GetRealSize() const

Return realSize of array.

template unsigned int astrix::Array< T >::GetDimension() const

Return number of dimensions.

template<>
void astrix::Array<T>::SetSizeHost(unsigned int _size)

Set size of Array on host.

template<>
void astrix::Array<T>::SetSizeDevice(unsigned int _size)

Set size of Array on device.

template void astrix::Array< T >::SetSize(unsigned int _size)

Set size of Array, either on host or device depending on cudaFlag.

template<>
void astrix::Array<T>::SetToValue(T value)

Set all Array entries to value.

template<>
void astrix::Array<T>::SetToValue(T value, unsigned int startIndex, unsigned int endIndex)

Set Array entries from startIndex to endIndex to value.

template void astrix::Array< T >::SetToSeries()

Set a[0] = 0, a[1] = 1, etc.

template void astrix::Array< T >::SetToSeries(unsigned int startIndex, unsigned int endIndex)

Set a[i] = i for all i from startIndex to endIndex.

template void astrix::Array< T >::SetEqual(const Array * B)

Set all entries equal to those in Array *B.

template void astrix::Array< T >::SetEqual(const Array * B, unsigned int N, unsigned int M)

Set all entries of dimension N of Array equal to dimension M of Array B.

template void astrix::Array< T >::SetEqual(const Array * B, int startPosition)

Set all entries starting from startPosition equal to those of Array *B, i.e. a[startPosition] = b[0] etc.

template void astrix::Array< T >::Reindex(unsigned int * reindex)

Reindex array: a[i] = a[reindex[i]].

template void astrix::Array< T >::Reindex(unsigned int * reindex, unsigned int N)

Reindex array: a[i] = a[reindex[i]] for the first N elements.

template<>
void astrix::Array<T>::InverseReindex(unsigned int *reindex, int maxValue, bool ignoreValue)

Inverse reindex array: a[i] = reindex[a[i]].

Set a[i] = reindex[a[i]]. If a[i] = -1, leave it at -1. If a[i] >= maxValue, subtract maxValue n times until a[i] < maxValue, and set a[i] = a[reindex[a[i]-n*maxValue]] + n*maxValue

template void astrix::Array< T >::Compact(int nKeep, Array < int > * keepFlag, Array < int > * keepFlagScan)

Compact; keep only entries where keepFlag == 1.

template void astrix::Array< T >::CopyToDevice()

Copy data from host to device.

template void astrix::Array< T >::CopyToHost()

Copy data from device to host.

template<>
void astrix::Array<T>::SetSingleValue(T value, int position)

Set a[position] = value.

template<>
void astrix::Array<T>::SetSingleValue(T value, int position, unsigned int N)

Set a[position] = value for dimension N.

template<>
void astrix::Array<T>::GetSingleValue(T *value, int position)

Real a[position] into *value.

template<>
void astrix::Array<T>::GetSingleValue(T *value, int position, unsigned int N)

Read a[position] for dimension N into *value.

template void astrix::Array< T >::AddValue(T value, unsigned int startIndex, unsigned int endIndex)

Add value to all entries from startIndex to endIndex.

template void astrix::Array< T >::Sort(Array < T > * arrayB)

Sort array, together with arrayB.

template <class S>
template<>
void astrix::Array<T>::SortByKey(Array<S> *indexArray)

Create index array for sorting.

template <class S>
template<>
void astrix::Array<T>::SortByKey(Array<S> *indexArray, unsigned int N)

Create index array for sorting dimension N.

template<>
T astrix::Array<T>::ExclusiveScan(Array<T> *result)

Perform exclusive scan.

template<>
T astrix::Array<T>::ExclusiveScan(Array<T> *result, unsigned int N)

Perform exclusive scan on dimension N.

template void astrix::Array< T >::Gather(Array < T > * in, Array < int > * map, int maxIndex)

Set out[i] = in[map[i]].

template void astrix::Array< T >::GatherIf(Array < T > * in, Array < int > * map, int value, int maxIndex)

If map[i] != value then out[i] = in[map[i]].

template void astrix::Array< T >::Scatter(Array < T > * in, Array < int > * map, int maxIndex)

Set out[map[i]] = in[i].

template <class S>
template void astrix::Array< T >::ScatterSeries(Array < S > * map, unsigned int maxIndex)

Set out[map[i]] = i.

template void astrix::Array< T >::SetToRandom()

Set to random values using rand()

template void astrix::Array< T >::SetToDiff(Array < T > * A, Array < T > * B)

Set a[i] = a[i] - b[i].

template int astrix::Array< T >::Minimum()

Return minimum of array.

template double astrix::Array< T >::Minimum(int N)

Return minimum of dimension N of array.

template int astrix::Array< T >::Maximum()

Return maximum of array.

template double astrix::Array< T >::Maximum(int N)

Return maximum of dimension N of array.

template unsigned int astrix::Array< T >::Sum()

Return sum of elements.

template void astrix::Array< T >::Invert()

Swap values 0 and 1.

template void astrix::Array< T >::Concat(Array < T > * A)

Join with Array A.

template int astrix::Array< T >::RemoveEvery(int start, int step)

Remove every entry start+i*step, compact array.

template <class S>
template<>
int astrix::Array<T>::RemoveEvery(int start, int step, Array<S> *A)

Remove every entry start+i*step, compact array and inverse reindex A.

template<>
int astrix::Array<T>::RemoveValue(T value)

Remove entries equal to value from Array.

template void astrix::Array< T >::ScatterUnique(Array < int > * A, Array < int > * B, int maxIndex, int ignoreValue, T value)

At a unique entry of A, i, (ignoring ignoreValue) set hostVec[B[i]] = value.

template void astrix::Array< T >::Shuffle()

Shuffle array (non-random!)

template<>
T *astrix::Array<T>::GetHostPointer()
const

Return pointer to host memory.

template<>
T *astrix::Array<T>::GetHostPointer(unsigned int _dim)
const

Return pointer to host memory for dimension _dim.

template<>
T *astrix::Array<T>::GetDevicePointer()
const

Return pointer to device memory.

template<>
T *astrix::Array<T>::GetDevicePointer(unsigned int _dim)
const

Return pointer to device memory for dimension _dim.

template<>
T *astrix::Array<T>::GetPointer()
const

Return pointer to either host or device memory, depending on cudaFlag.

template<>
T *astrix::Array<T>::GetPointer(unsigned int _dim)
const

Return pointer to either host or device memory for dimension _dim.

Public Static Attributes

template<>
int64_t astrix::Array<T>::memAllocatedHost

Total amount of memory (bytes) allocated on host in all Array‘s.

template<>
int64_t astrix::Array<T>::memAllocatedDevice

Total amount of memory (bytes) allocated on device in all Array‘s.

Indices and tables