SOLVCON¶

SOLVCON is a collection of conservation-law solvers that use the space-time Conservation Element and Solution Element (CESE) method [Chang95]. The equations to be solved are formulated as:

$\dpd{\mathbf{u}}{t} + \sum_{k=1}^3 \mathrm{A}^{(k)}(\mathbf{u})\dpd{\mathbf{u}}{x_k} = \mathbf{s}(\mathbf{u})$

where $\mathbf{u}$ is the unknown vector, $\mathrm{A}^{(1)}$ , $\mathrm{A}^{(2)}$ , and $\mathrm{A}^{(3)}$ the Jacobian matrices, and $\mathbf{s}$ the source term.

Install¶

Clone the hg repository from https://bitbucket.org/solvcon/solvcon:

$ hg clone https://bitbucket.org/solvcon/solvcon

SOLVCON needs the following packages: gcc 4.3+ (clang on OSX works as well), SCons 2+, Python 2.7, Cython 0.16+, Numpy 1.5+, LAPACK, NetCDF 4+, SCOTCH 6.0+, Nose 1.0+, Paramiko 1.14+, boto 2.29.1+, gmsh 2.5+, and VTK 5.6+.

In the contrib/ directory, you can find the scripts for installing these dependencies:

aptget.*.sh for Debian/Ubuntu
conda.sh for Miniconda/Anaconda

The binary part of SOLVCON should be built with SCons:

$ scons scmods

Then add the installation path to the environment variables $PATH and $PYTHONPATH.

Additional build and tests:

Building document requires Sphinx 1.1.2+, Sphinxcontrib issue tracker 0.11, and graphviz 2.28+. Once the binary of SOLVCON is built, the following commands can build the document:
```
$ cd doc
$ make html
```
The built document will be available at doc/build/html/.
Unit tests should be run with Nose:
```
$ nosetests
```
Another set of tests are collected in ftests/ directory, and can be run with:
```
$ nosetests ftests/*
```
Some tests in ftests/ involve remote procedure call (RPC) that uses ssh. You need to set up the public key authentication to properly run them.

Documents¶

Solver Architecture¶

SOLVCON is built upon two keystones: (i) unstructured meshes for spatial discretization and (ii) two-level loop structure of partial differential equation (PDE) solvers.

Unstructured Meshes¶

We usually discretize the spatial domain of interest before solving PDEs with digital computers. The discretized space is called a mesh [Mavriplis97]. When discretization is done by exploiting regularity in space, like cutting along each of the Cartesian coordinate axes, the discretized space is called a structured mesh. If the discretization does not follow any spatial order, we call the spatial domain an unstructured mesh. Both meshing strategies have their strength and weakness. Sometimes a structured mesh is also call a grid. Numerical methods that rely on spatial discretization are called mesh-based or grid-based. Most PDE-solving methods in production uses are mesh-based, but meshless methods have their advantages.

To accommodate complex geometry, SOLVCON chose to use unstructured meshes of mixed elements. Because no structure is assumed for the geometry to be modeled, the mesh can be automatically generated by using computer programs. For example, the following image shows a triangular mesh of a two-dimensional irregular domain:

Figure 1: Two-dimensional sample mesh

which is generated by using the Gmsh commands listed in ustmesh_2d_sample.geo. On the other hand, creation of structured meshes often needs a large amount of manual operations and will not be discussed in this document.

In SOLVCON, we assume a mesh is fully covered by a finite number of non-overlapping sub-regions, and only composed by these sub-regions. The sub-regions are called mesh elements. In one-dimensional space, SOLVCON also defines one type of mesh elements, line, as shown in Figure One-dimensional mesh element.

Figure 2: One-dimensional mesh element

SOLVCON allows two types of two-dimensional mesh elements, quadrilaterals and triangles, as shown in Figure Two-dimensional mesh elements, and four types of three-dimensional mesh elements, hexahedra, tetrahedra, prisms, and pyramids, as shown in Figure Three-dimensional mesh elements.

Figure 3: Two-dimensional mesh elements

Figure 4: Three-dimensional mesh elements

The numbers annotated in the figures are the order of the vertices of the elements. A SOLVCON mesh can be a mixture of elements of the same dimension, although it is often composed of one type of element. Two modules provide the support of the meshes: (i) solvcon.block defines and manages various look-up tables that form the data structure of the mesh in Python, and (ii) solvcon.mesh serves as the interface of the mesh data in C.

Entities¶

Before explaining the data structure of the meshes, we need to introduce some basic terminologies and definitions. In SOLVCON, a cell means a mesh element. The dimensionality of a cell equals to that of the mesh it belongs to, e.g., a two-dimensional mesh is composed by two-dimensional cells. A cell is assumed to be concave, and enclosed by a set of faces. The dimensionality of a face is one less than that of a cell. A face is also assumed to be concave, and formed by connecting a sequence of nodes. The dimensionality of a node is at least one less than that of a face. Cells, faces, and nodes are the basic constructs, which we call entities, of a SOLVCON mesh.

Defining the term “entity” for SOLVCON facilitates a unified treatment of two- and three-dimensional meshes and the corresponding solvers [1]. A cell can be either two- or three-dimensional, and the associated faces become one- or two-dimensional, respectively. Because a face is either one- or two-dimensional, it can always be formed by a sequence of points, which is zero-dimensional. In this treatment, a point is equivalent to a node defined in the previous passage.

Take the two-dimensional mesh shown above as an example, triangular elements are used as the cells. The triangles are formed by three lines (one-dimensional shapes), which are the faces. Each line has two points (zero-dimensional). If we have a three-dimensional mesh composed by hexahedral cells, then the faces should be quadrilaterals (two-dimensional shapes).

All the mesh elements supported by SOLVCON are listed in the following table. The first column is the name of the element, and the second column is the type ID used in SOLVCON. The third column lists the dimensionality. The fourth, fifth, and sixth columns show the number of zero-, one-, and two-dimensional sub-entities belong to the element type, respectively. Note that the terms “point” and “line” appear in both the first row and first column, for they are the only element type in the space of the corresponding dimensionality.

Name	Type	Dim	Point#	Line#	Surface#
Point	0	0	0	0	0
Line	1	1	2	0	0
Quadrilateral	2	2	4	4	0
Triangle	3	2	3	3	0
Hexahedron	4	3	8	12	6
Tetrahedron	5	3	4	4	4
Prism	6	3	6	9	5
Pyramid	7	3	5	8	5

Although SOLVCON doesn’t support one-dimensional solvers, for completeness, we define the relation between one-dimensional cells (lines) and their sub-entities as:

Shape (type)	Face	= Point
Line (0)	0	$\cdot$ 0
Line (0)	1	$\cdot$ 1

That is, as shown in Figure One-dimensional mesh element, a one-dimensional “cell” (line) has two “faces”, which are essentially point 0 and point 1. Symbol $\cdot$ indicates a point.

It will be more practical to illustrate the relation between two-dimensional cells and their sub-entities in a table (see Figure Two-dimensional mesh elements for point locations):

Shape (type)	Face	= Line formed by points
Quadrilateral (2)	0	$\diagup$ 0 1
	1	$\diagup$ 1 2
	2	$\diagup$ 2 3
	3	$\diagup$ 3 0
Triangle (3)	0	$\diagup$ 0 1
	1	$\diagup$ 1 2
	2	$\diagup$ 2 0

Symbol $\diagup$ indicates a line. The orientation of lines of each two-dimensional shape is defined to follow the right-hand rule. The shape enclosed by the lines has an area normal vector points to the direction of $+z$ (outward paper/screen).

The relation between three-dimensional cells and their sub-entities is defined in the table (see Figure Three-dimensional mesh elements for point locations):

Shape (type)	Face	= Surface formed by points
Hexahedron (4)	0	$\square$ 0 3 2 1
	1	$\square$ 1 2 6 5
	2	$\square$ 4 5 6 7
	3	$\square$ 0 4 7 3
	4	$\square$ 0 1 5 4
	5	$\square$ 2 3 7 6
Tetrahedron (5)	0	$\triangle$ 0 2 1
	1	$\triangle$ 0 1 3
	2	$\triangle$ 0 3 2
	3	$\triangle$ 1 2 3
Prism (6)	0	$\triangle$ 0 1 2
	1	$\triangle$ 3 5 4
	2	$\square$ 0 3 4 1
	3	$\square$ 0 2 5 3
	4	$\square$ 1 4 5 2
Pyramid (7)	0	$\triangle$ 0 4 3
	1	$\triangle$ 1 4 0
	2	$\triangle$ 2 4 1
	3	$\triangle$ 3 4 2
	4	$\square$ 0 3 2 1

Symbol $\square$ indicates a quadrilateral, while symbol $\triangle$ indicates a triangle.

Because a face is associated to two adjacent cells unless it’s a boundary face, it needs to identify to which cell it belongs, and to which cell it is neighbor. The area normal vector of a face is always point from the belonging cell to neighboring cell. The same rule applies to faces of two-dimensional meshes (lines) too.

Data Structure Defined in `solvcon.block`¶

Real data of unstructured meshes are stored in module solvcon.block. A simple table for all element types is defined as elemtype:

solvcon.block.elemtype¶: A numpy.ndarray object of shape (8, 5) and type int32. This array is a reference table for element types in SOLVCON. The content is shown in the first table in Section Entities. Each row represents an element type. The first column is the index of the element type, the second the dimensionality, the third column the number of points, the fourth the number lines, and the fifth the number of surfaces.

Class Block contains descriptive information, look-up tables, and other miscellaneous information for a SOLVCON mesh. There are three steps required to fully construct a Block object: (i) instantiation, (ii) definition, and (iii) build-up. In the first step, when instantiating an object, shape information must be provided to the constructor to allocate arrays for look-up tables:

from solvcon.block import Block
blk = Block(ndim=2, nnode=4, ncell=3)

Second, we fill the definition of the look-up tables into the object. We at least need to provide the node coordinates and the node lists of cells:

blk.ndcrd[:,:] = (0,0), (-1,-1), (1,-1), (0,1)
blk.cltpn[:] = 3
blk.clnds[:,:4] = (3, 0,1,2), (3, 0,2,3), (3, 0,3,1)

Third and finally, we build up the rest of the object by calling:

blk.build_interior()
blk.build_boundary()
blk.build_ghost()

By running the additional code, the block can be saved as a VTK file for viewing:

from solvcon.io.vtkxml import VtkXmlUstGridWriter
iodev = VtkXmlUstGridWriter(blk)
iodev.write('block_2d_sample.vtu')

Figure 5: A simple Block object

class solvcon.block.Block(ndim=0, nnode=0, nface=0, ncell=0, nbound=0, use_incenter=False)¶

This class represents the unstructured meshes used in SOLVCON. As such, in SOLVCON, an unstructured mesh is also called a block. The following six attributes can be passed into the constructor. ndim, nnode, and ncell need to be non-zero to instantiate a valid block. nface and nbound might be different to the given value after building up the object. use_incenter is an optional flag.

ndim¶

Type:	`int`

Number of dimensionalities of this mesh. Read only after instantiation.

nnode¶

Type:	`int`

Total number of (non-ghost) nodes of this mesh. Read only after instantiation.

nface¶

Type:	`int`

Total number of (non-ghost) faces of this mesh. Read only after instantiation.

ncell¶

Type:	`int`

Total number of (non-ghost) cells of this mesh. Read only after instantiation.

nbound¶

Type:	`int`

Total number of boundary faces or ghost cells of this mesh. Read only after instantiation.

use_incenter¶

Type:	`bool`

Indicates calculating incenters instead of centroids for cells. Default is False (using centroids of cells).

To construct a block object, SOLVCON needs to know the dimensionalities (ndim), the number of nodes (nnode), faces (nface), and cells (ncell), and the number of boundary faces (nbound) of the mesh. These keyword parameters are taken to initialize the following properties:

The meshes are mainly defined by three sets of look-up tables (arrays). The first set is the geometry arrays, which store the coordinate values of mesh elements:

ndcrd¶: Coordinates of nodes. It’s a two-dimensional numpy.ndarray array of shape (nnode, ndim) of type float64.

fccnd¶: Centroids of faces. It’s a two-dimension numpy.ndarray of shape (nface, ndim) of type float64.

fcnml¶: Unit normal vectors of faces. It’s a two-dimension numpy.ndarray of shape (nface, ndim) of type float64.

fcara¶: Areas of faces. The value should always be non-negative. It’s a one-dimension numpy.ndarray of shape (nface,) of type float64.

clcnd¶: Centroids of cells. It’s a two-dimension numpy.ndarray of shape (ncell, ndim) of type float64.

clvol¶: Volumes of cells. It’s a one-dimension numpy.ndarray of shape (ncell,) of type float64.

The second set is the meta-data or type data arrays:

fctpn¶: Type ID of faces. It’s a one-dimensional numpy.ndarray of shape (nface,) of type int32.

cltpn¶: Type ID of cells. It’s a one-dimensional numpy.ndarray of shape (ncell,) of type int32.

clgrp¶: Group ID of cells. It’s a one-dimensional numpy.ndarray of shape (ncell,) of type int32. For a new Block object, it should be initialized with -1.

The third and last set is the connectivity arrays:

fcnds¶: Lists of the nodes of each face. It’s a two-dimensional numpy.ndarray of shape (nface, FCMND+1) and type int32.

fccls¶: Lists of the cells connected by each face. It’s a two-dimensional numpy.ndarray of shape (nface, 4) and type int32.

clnds¶: Lists of the nodes of each cell. It’s a two-dimensional numpy.ndarray of shape (ncell, CLMND+1) and type int32.

clfcs¶: Lists of the faces of each cell. It’s a two-dimensional numpy.ndarray of shape (ncell, CLMFC+1) and type int32.

Every look-up array has two associated arrays distinguished by different prefixes: (i) gst (denoting for “ghost”) and (ii) sh (denoting for “shared”). SOLVCON uses the technique of ghost cells to treat boundary conditions [Mavriplis97], and the gst arrays store the information for ghost cells. However, to facilitate efficient indexing in solvers, each of the ghost arrays should be put in a continuous block of memory adjacent to its interior counterpart. In SOLVCON, the sh arrays are the continuous memory blocks for both ghost and interior look-up tables, and a pair of gst and normal arrays is simply the views of two consecutive, non-overlapping sub-regions of a memory block. More details of the technique of ghost cells will be given in module solvcon.mesh.

There are some attributes associated with ghost cells:

ngstnode¶

Type:	`int`

Number of nodes only associated with ghost cells. Only valid after build-up. Read only.

ngstface¶

Type:	`int`

Number of faces only associated with ghost cells. Only valid after build-up. Read only.

ngstcell¶

Type:	`int`

Number of ghost cells. Only valid after build-up. Read only.

Three arrays need to be defined before we can build up a Block object: (i) ndcrd, (ii) cltpn, and (iii) clnds. With these information, build_interior() builds up the interior arrays for a Block object. build_boundary() then organizes the information for boundary conditions. Finally, build_ghost() builds up the shared and ghost arrays for the Block object. Only after the build-up process, the Block object can be used by solvers.

build_interior()¶

Returns:	Nothing.

Building up a Block object includes two steps. First, the method extracts arrays clfcs, fctpn, fcnds, and fccls from the defined arrays cltpn and clnds. If the number of extracted faces is not the same as that passed into the constructor, arrays related to faces are recreated.

Second, the method calculates the geometry information and fills the corresponding arrays.

build_boundary(unspec_type=None, unspec_name='unspecified')¶

Parameters:	unspec_type (`type`) – BC type for the unspecified boundary faces. Set to `None` indicates the default to `solvcon.boundcond.unspecified`. unspec_name (`str`) – Name for the unspecified BC.
Returns:	Nothing.

This method iterates over each of the solvcon.boundcond.BC objects listed in bclist to collect boundary-condition information and build boundary faces. If a face belongs to only one cell (i.e., has no neighboring cell), it is regarded as a boundary face.

Unspecified boundary faces will be collected to form an additional solvcon.boundcond.BC object. It sets bndfcs for later use by build_ghost().

build_ghost()¶

Returns:	Nothing.

This method creates the shared arrays, calculates the information for ghost cells, and reassigns interior arrays as the right portions of the shared arrays.

A Block object also contains three instance variables for boundary-condition treatments:

bclist¶

Type:	`list`

The list of associated solvcon.boundcond.BC objects.

nbound

Type:	`int`

Number of boundary faces. Only valid after build-up. It should equals to ngstcell.

bndfcs¶

Type:	`numpy.ndarray`

The array is of shape (nbound, 2) and type int32. Each row contains the data for a boundary face. The first column is the 0-based index of the face, while the second column is the serial number of the associated solvcon.boundcond.BC object.

create_msh()¶

Returns:	An object contains the `sc_mesh_t` variable for C code to use data in the `Block` object.
Return type:	`solvcon.mesh.Mesh`

The following code shows how and when to use this method:

>>> blk = Block(ndim=2, nnode=4, nface=6, ncell=3, nbound=3)
>>> blk.ndcrd[:,:] = (0,0), (-1,-1), (1,-1), (0,1)
>>> blk.cltpn[:] = 3
>>> blk.clnds[:,:4] = (3, 0,1,2), (3, 0,2,3), (3, 0,3,1)
>>> blk.build_interior()
>>> # it's OK to get a msh when its content is still invalid.
>>> msh = blk.create_msh()
>>> blk.build_boundary()
>>> blk.build_ghost()
>>> # now the msh is valid for the blk is fully built-up.
>>> msh = blk.create_msh()

In class Block there are also useful constants defined:

Block.FCMND¶

Type: int

The maximum number of nodes that a face can have. From the first table in Section Entities, its value should be 4.

Block.CLMND¶

Type: int

The maximum number of nodes that a cell can have. From the first table in Section Entities, its value should be 8.

Block.CLMFC¶

Type: int

The maximum number of faces that a cell can have. From the first table in Section Entities, its value should be 6.

Low-Level Interface to C Defined in `solvcon.mesh`¶

Although it is convenient to have data structure defined in the Python module solvcon.block, kernel of numerical methods are usually implemented in C. To bridge Python and C, we use Cython to write an interfacing module solvcon.mesh. This module enables C code to use the mesh data held by a solvcon.block.Block object, and allows Python to use those C functions.

A header file mesh.h contains the essential declarations to use the mesh data:

sc_mesh_t¶

This struct is the counterpart of the Python class solvcon.block.Block in C. It contains four sections of fields in order.

The first field section is for shape. These fields correspond to the instance properties (attributes) in solvcon.block.Block of the same names:

int ndim¶

int nnode¶

int nface¶

int ncell¶

int nbound¶

int ngstnode¶

int ngstface¶

int ngstcell¶

The second field section is for geometry arrays. These fields correspond to the instance variables (attributes) in solvcon.block.Block of the same names:

Note

All arrays in sc_mesh_t are shared arrays but the pointers point to the start of their interior portion. In this way, access to ghost information can be efficiently done by using negative indices of nodes, faces, and cells in the first dimension of these arrays. But negative indices in higher dimensions of the arrays is meaningless.

double* ndcrd¶

double* fccnd¶

double* fcnml¶

double* fcara¶

double* clcnd¶

double* clvol¶

The third field section is for type/meta arrays. These fields correspond to the instance variables (attributes) in solvcon.block.Block of the same names:

int* fctpn¶

int* cltpn¶

int* clgrp¶

The fourth and final field section is for connectivity arrays. These fields correspond to the instance variables (attributes) in solvcon.block.Block of the same names:

int* fcnds¶

int* fccls¶

int* clnds¶

int* clfcs¶

The SOLVCON C library (libsolvcon.a) contains five mesh-related functions that are used internally in Mesh. These functions are not meant to be part of the interface, but can be a reference about the usage of sc_mesh_t:

int sc_mesh_extract_faces_from_cells(sc_mesh_t *msd, int mface, int *pnface, int *clfcs, int *fctpn, int *fcnds, int *fccls)¶

This function extracts interior faces from the node lists of the cells given in the first argument msd. The second argument mface is also an input, which sets the maximum value of possible number of faces to be extracted.

The rest of the arguments is outputs. The arrays pointed by the last four arguments need to be pre-allocated with appropriate size or the memory will be corrupted.

int sc_mesh_calc_metric(sc_mesh_t *msd, int use_incenter)¶: This function calculates the geometry information and stores the calculated values into the arrays specified in msd. The second argument use_incenter is a flag. When it is set to 1, the function calculates and stores the incenter of the cells. Otherwise, the function calculates and stores the centroids of the cells.

void sc_mesh_build_ghost(sc_mesh_t *msd, int *bndfcs)¶

Build all information for ghost cells by mirroring information from interior cells. The arrays in the first argument msd will be altered, but data in the second argument bndfcs will remain intact. The action includes:

Define indices and build connectivities for ghost nodes, faces, and cells. In the same loop, mirror the coordinates of interior nodes to ghost nodes.
Compute center coordinates for faces for ghost cells.
Compute normal vectors and areas for faces for ghost cells.
Compute center coordinates for ghost cells.
Compute volume for ghost cells.

It should be noted that all the geometry, type/meta and connectivity data used in this function are SHARED arrays rather than interior arrays. The indices for ghost information should be carefully treated. All the ghost indices are negative in shared arrays.

int sc_mesh_build_rcells(sc_mesh_t *msd, int *rcells, int *rcellno)¶: This is a utility function used by Mesh.create_csr(). The first argument msd is input and will not be changed, and the output will be write to the second and third arguments, rcells and rcellno. Sufficient memory must be pre-allocated for the output arrays before calling or memory can be corrupted.

int sc_mesh_build_csr(sc_mesh_t *msd, int *rcells, int *adjncy)¶: This is a utility function used by Mesh.create_csr(). The first argument msd and the second argument rcells are input and will not be changed, while the third argument adjncy is output. Sufficient memory must be pre-allocated for the output array before calling or memory can be corrupted.

A Python class Mesh is written by using Cython to convert a Python-space solvcon.block.Block object into a sc_mesh_t struct variable for use in C. This class is meant to be subclassed to implement the core number-crunching algorithm of a numerical method. In addition, this class also provides functionalities that need the C utility functions listed above.

class solvcon.mesh.Mesh¶

This class associates the C functions for mesh operations to the mesh data and exposes the functions to Python.

_msd¶: This attribute holds a C struct sc_mesh_t for internal use.

setup_mesh(blk)¶

Parameters:	blk (`solvcon.block.Block`) – The block object that supplies the information.

extract_faces_from_cells(max_nfc)¶

Parameters:	max_nfc (C `int`) – Maximum value of possible number of faces to be extracted.
Returns:	Four interior `numpy.ndarray` for `solvcon.block.Block.clfcs`, `solvcon.block.Block.fctpn`, `solvcon.block.Block.fcnds`, and `solvcon.block.Block.fccls`.

Internally calls sc_mesh_extract_face_from_cells().

calc_metric()¶

Returns:	Nothing.

Calculates geometry information including normal vector and area of faces, and centroid/incenter coordinates and volume of cells. Internally calls sc_mesh_calc_metric().

build_ghost()¶

Returns:	Nothing.

Builds data for ghost cells. Internally calls sc_mesh_build_ghost().

create_csr()¶

Returns:	xadj, adjncy
Return type:	tuple of `numpy.ndarray`

Builds the connectivity graph in the CSR (compressed storage format) used by SCOTCH/METIS. Internally calls sc_mesh_build_rcells() and sc_mesh_build_csr().

partition(npart, vwgtarr=None)¶

Parameters:	npart (C `int`) – Number of parts to be partitioned to. vwgtarr (`numpy.ndarray`) – vwgt weighting settings. Default is `None`.
Returns:	A 2-tuple of (i) number of cut edges for the partitioning and (ii) a `numpy.ndarray` of shape (`solvcon.block.Block.ncell`,) and type `int32` that indicates the partition number of each cell in the mesh.
Return type:	`int`, `numpy.ndarray`

Internally calls METIS_PartGraphKway() of the SCOTCH library for mesh partitioning.

Numerical Code¶

The numerical calculations in SOLVCON rely on exploiting a two-level loop structure, i.e., the temporal loop and the spatial loops. For time-accurate solvers, there is always an outer loop that coordinates the time-marching. The outer loop is called the temporal loop, and it should be implemented in subclasses of MeshCase. Inside the temporal loop, there can be one or many inner loops that calculate the new values of the fields. The inner loops are called the spatial loops, and they should be implemented in subclasses of MeshSolver.

The outer temporal loop is more responsible for coordinating, while the inner spatial loops is closer to numerical algorithms. These two levels allow us to segregate code. An object of MeshCase can be seen as the realization of a simulation case in SOLVCON (as a convention the object’s name should contain or just be cse). Code in MeshCase is mainly about obtaining settings, input and output, and provision of the execution environment. On the other hand, we implement the numerical algorithm in MeshSolver to manipulate the field data (as a convention the object’s name should contain or just be svr). Its code shouldn’t involve input nor output (excepting that for debugging) but needs to take parallelism into account.

Code for data processing should go to MeshHook (as a convention the objects should be named with hok), which is the companion of MeshCase. Code that processes data close to numerical methods should go to MeshAnchor (as a convention the objects should be named with ank), which is the companion of MeshSolver.

In this section, for conciseness, the terms solver, anchor, case, and hook are used to denote the classes MeshSolver, MeshAnchor, MeshCase, and MehsHook or their instances, respectively. In the issue tracking system, solver, anchor, case, and hook form a component “sach”.

`solvcon.case`¶

Module solvcon.case contains code for making a simulation case (subclasses of solvcon.case.MeshCase). Because a case coordinates the whole process of a simulation run, for parallel execution, there can be only one MeshCase object residing in the controller (head) node.

By the design, MeshCase itself cannot be directly used. It must be subclassed to implement control logic for a specific application. The application can be a concrete model for a certain physical process, or an abstraction of a group of related physical processes, which can be further subclassed.

class solvcon.case.MeshCase(**kw)¶

Base class for simulation cases based on solvcon.mesh.Mesh.

Inheritance diagram of MeshCase

init() and run() are the two primary methods responsible for the execution of the simulation case object. Both methods accept a keyword parameter “level”:

run level 0: fresh run (default),
run level 1: restart run,
run level 2: initialization only.

cleanup(signum=None, frame=None)¶

Parameters:	signum – Signal number. frame – Current stack frame.

A signal handler for cleaning up the simulation case on termination or when errors occur. This method can be overridden in subclasses. The base implementation is trivial, but usually doesn’t need to be overridden.

An example to connect this method to a signal:

>>> from .testing import create_trivial_2d_blk
>>> from .solver import MeshSolver
>>> blk = create_trivial_2d_blk()
>>> cse = MeshCase(basefn='meshcase', mesher=lambda *arg: blk,
...                domaintype=domain.Domain, solvertype=MeshSolver)
>>> cse.info.muted = True
>>> signal.signal(signal.SIGTERM, cse.cleanup)
0

An example to call this method explicitly:

>>> cse.init()
>>> cse.run()
>>> cse.cleanup()

defer(delayed, replace=False, **kw)¶

Insert (append or replace) hooks.

Parameters:	delayed (`solvcon.MeshHook` or `solvcon.MeshAnchor`.) – The delayed construct. replace (bool) – True if existing object should be replaced.
Returns:	Nothing.

>>> import solvcon as sc
>>> from solvcon.testing import create_trivial_2d_blk
>>> cse = MeshCase() # No arguments because of demonstration.
>>> len(cse.runhooks)
0
>>> # Insert a hook.
>>> cse.defer(sc.MeshHook, dummy='name1')
>>> cse.runhooks[0].kws['dummy']
'name1'
>>> # Insert the second hook to replace the first one.
>>> cse.defer(sc.MeshHook, replace=True, dummy='name2')
>>> cse.runhooks[0].kws['dummy'] # Got replaced.
'name2'
>>> len(cse.runhooks) # Still only one hook in the list.
1
>>> # Insert the third hook without replace.
>>> cse.defer(sc.MeshHook, dummy='name3')
>>> cse.runhooks[1].kws['dummy'] # Got replaced.
'name3'

Initialize¶

MeshCase.init(level=0)¶

Parameters:	level (`int`) – Run level; higher level does less work.
Returns:	Nothing

Load a block and initialize the solver from the geometry information in the block and conditions in the self case. If parallel run is specified (through domaintype), split the domain and perform corresponding tasks.

For a MeshCase to be initialized, some information needs to be supplied to the constructor:

>>> cse = MeshCase()
>>> cse.info.muted = True
>>> cse.init()
Traceback (most recent call last):
    ...
TypeError: coercing to Unicode: need string or buffer, NoneType found

Mesh information. We can provide meshfn that specifying the path of a valid mesh file, or provide mesher, which is a function that generates the mesh and returns the solvcon.block.Block object, like the following code:

>>> from solvcon.testing import create_trivial_2d_blk
>>> blk = create_trivial_2d_blk()
>>> cse = MeshCase(mesher=lambda *arg: blk)
>>> cse.info.muted = True
>>> cse.init() 
Traceback (most recent call last):
    ...
TypeError: isinstance() arg 2 must be a class, type, or tuple of
classes and types

Type of the spatial domain. This information is used for detemining sequential or parallel execution, and performing related operations:

>>> cse = MeshCase(mesher=lambda *arg: blk, domaintype=domain.Domain)
>>> cse.info.muted = True
>>> cse.init()
Traceback (most recent call last):
    ...
TypeError: 'NoneType' object is not callable

The type of solver. It is used to specify the underlying numerical method:

>>> from solvcon.solver import MeshSolver
>>> cse = MeshCase(mesher=lambda *arg: blk,
...                domaintype=domain.Domain,
...                solvertype=MeshSolver)
>>> cse.info.muted = True
>>> cse.init()
Traceback (most recent call last):
    ...
TypeError: cannot concatenate 'str' and 'NoneType' objects

The base name. It is used to name its output files:

>>> cse = MeshCase(
...     mesher=lambda *arg: blk, domaintype=domain.Domain,
...     solvertype=MeshSolver, basefn='meshcase')
>>> cse.info.muted = True
>>> cse.init()

Time-March¶

MeshCase.run(level=0)¶

Parameters:	level (`int`) – Run level; higher level does less work.
Returns:	Nothing

Temporal loop for the incorporated solver. A simple example:

>>> from .testing import create_trivial_2d_blk
>>> from .solver import MeshSolver
>>> blk = create_trivial_2d_blk()
>>> cse = MeshCase(basefn='meshcase', mesher=lambda *arg: blk,
...                domaintype=domain.Domain, solvertype=MeshSolver)
>>> cse.info.muted = True
>>> cse.init()
>>> cse.run()

Arrangement¶

solvcon.case.arrangements¶: The module-level registry for arrangements.

MeshCase.arrangements¶: The class-level registry for arrangements.

classmethod MeshCase.register_arrangement(func, casename=None)¶

Returns:	Simulation function.
Return type:	callable

This class method is a decorator that creates a closure (internal function) that turns the decorated function to an arrangement, and registers the arrangement into the module-level registry and the class-level registry. The decorator function should return a MeshCase object cse, and the closure performs a simulation run by the following code:

try:
    signal.signal(signal.SIGTERM, cse.cleanup)
    signal.signal(signal.SIGINT, cse.cleanup)
    cse.init(level=runlevel)
    cse.run(level=runlevel)
    cse.cleanup()
except:
    cse.cleanup()
    raise

The usage of this decorator can be exemplified by the following code, which creates four arrangements (although the first three are erroneous):

>>> @MeshCase.register_arrangement
... def arg1():
...     return None
>>> @MeshCase.register_arrangement
... def arg2(wrongname):
...     return None
>>> @MeshCase.register_arrangement
... def arg3(casename):
...     return None
>>> @MeshCase.register_arrangement
... def arg4(casename):
...     from .testing import create_trivial_2d_blk
...     from .solver import MeshSolver
...     blk = create_trivial_2d_blk()
...     cse = MeshCase(basefn='meshcase', mesher=lambda *arg: blk,
...                    domaintype=domain.Domain, solvertype=MeshSolver)
...     cse.info.muted = True
...     return cse

The created arrangements are collected to a class attribute arrangements, i.e., the class-level registry:

>>> sorted(MeshCase.arrangements.keys())
['arg1', 'arg2', 'arg3', 'arg4']

The arrangements in the class attribute arrangements are also put into a module-level attribute solvcon.case.arrangements:

>>> arrangements == MeshCase.arrangements
True

The first example arrangement is a bad one, because it allows no argument:

>>> arrangements.arg1()
Traceback (most recent call last):
  ...
TypeError: arg1() takes no arguments (1 given)

The second example arrangement is still a bad one, because although it has an argument, the name of the argument is incorrect:

>>> arrangements.arg2()
Traceback (most recent call last):
  ...
TypeError: arg2() got an unexpected keyword argument 'casename'

The third example arrangement is a bad one for another reason. It doesn’t return a MeshCase:

>>> arrangements.arg3()
Traceback (most recent call last):
  ...
AttributeError: 'NoneType' object has no attribute 'cleanup'

The fourth example arrangement is finally good:

>>> arrangements.arg4()

`solvcon.solver`¶

Module solvcon.solver provides the basic facilities for implementing numerical methods by subclassing MeshSolver. The base class is defined as:

class solvcon.solver.MeshSolver(blk, time=0.0, time_increment=0.0, enable_mesg=False, debug=False, **kw)¶

Base class for all solving code that take Mesh, which is usually needed to write efficient C/C++ code for implementing numerical methods.

Here’re some examples about using MeshSolver. The first example shows that we can’t directly use it. A vanilla MeshSolver can’t march:

>>> from . import testing
>>> svr = MeshSolver(testing.create_trivial_2d_blk())
>>> svr.march(0.0, 0.1, 1) 
Traceback (most recent call last):
    ...
TypeError: 'NoneType' object ...

At minimal we need to override the _MMNAMES class attribute:

>>> class DerivedSolver(MeshSolver):
...     _MMNAMES = MeshSolver.new_method_list()
>>> svr = DerivedSolver(testing.create_trivial_2d_blk())
>>> svr.march(0.0, 0.1, 1)
{}

Of course the above derived solver did nothing. Let’s see another example solver that does non-trivial things:

>>> class ExampleSolver(MeshSolver):
...     _MMNAMES = MeshSolver.new_method_list()
...     @_MMNAMES.register
...     def calcsomething(self, worker=None):
...         self.marchret['key'] = 'value'
>>> svr = ExampleSolver(testing.create_trivial_2d_blk())
>>> svr.march(0.0, 0.1, 1)
{'key': 'value'}

Two instance attributes are used to record the temporal information:

time = None¶: The current time of the solver. By default, time is initialized to 0.0, which is usually desired value. The default value can be overridden from the constructor.

time_increment = None¶: The temporal interval between the current and the next time steps. It is usually referred to as $\Delta t$ in the numerical literature. By default, time_increment is initialized to 0.0, but the default should be overridden from the constructor.

Four instance attributes are used to track the status of time-marching:

step_current = None¶: It is an int that records the current step of the solver. It is reset to 0 on every invokation of march().

step_global = None¶: It is similar to step_current, but persists over restart. Without restarts, step_global should be identical to step_current.

substep_run = None¶: The number of sub-steps that a single time step should be split into. It is initialized to 1 and should be overidden in subclasses if needed.

substep_current = None¶: The current sub-step of the solver. It is initialized to 0.

Derived classes of MeshSolver should use the following method new_method_list() to make a new class variable _MMNAMES to define numerical methods:

static new_method_list()¶

Returns:	An object to be set to `_MMNAMES`.
Return type:	`_MethodList`

In subclasses of MeshSolver, implementors can use this utility method to creates an instance of _MethodList, which should be set to _MMNAMES.

_MMNAMES = None¶: This class attribute holds the names of the methods to be called in march(). It is of type _MethodList. The default value is None and must be reset in subclasses by calling new_method_list().

Useful entities are attached to the class MeshSolver:

MeshSolver.ALMOST_ZERO = 1e-200¶: A positive floating-point number close to zero. The value is not DBL_MIN, which can be accessed through sys.float_info.

Time-Marching¶

MeshSolver.march(time_current, time_increment, steps_run, worker=None)¶

Parameters:	time_current (float) – Starting time of this set of marching steps. time_increment (float) – Temporal interval $\Delta t$ of the time step. steps_run (int) – The count of time steps to run.
Returns:	`marchret`

This method performs time-marching. The parameters time_current and time_increment are used to reset the instance attributes time and time_increment, respectively. In each invokation step_current is reset to 0.

There is a nested two-level loop in this method for time-marching. The outer loop iterates for time steps, and the inner loop iterates for sub time steps. The outer loop runs steps_run times, while the inner loop runs substep_run times. In total, the inner loop runs steps_run * substep_run times. In each sub time step (in the inner loop), the increment of the attribute time is time_increment/substep_run. The temporal increment per time step is effectively time_increment, with a slight error because of round-off.

Before entering and after leaving the outer loop, premarch and postmarch anchors will be run (through the attribute runanchors). Similarly, before entering and after leaving the inner loop, prefull and postfull anchors will be run. Inside the inner loop of sub steps, before and after executing all the marching methods, presub and postsub anchors will be run. Lastly, before and after invoking every marching method, a pair of anchors will be run. The anchors for a marching method are related to the name of the marching method itself. For example, if a marching method is named “calcsome”, anchor precalcsome will be run before the invocation, and anchor postcalcsome will be run afterward.

Derived classes can set marchret dictionary, and march() will return the dictionary at the end of execution. The dictionary is reset to empty at the begninning of the execution.

MeshSolver.marchret¶: Values to be returned by this solver. It will be set to a dict in march().

MeshSolver.runanchors¶: This attribute is of type MeshAnchorList, and the foundation of the anchor mechanism of SOLVCON. An MeshAnchorList object like this collects a set of MeshAnchor objects, and is callable. When being called, runanchors iterates the contained MeshAnchor objects and invokes the corresponding methods of the individual MeshAnchor.

MeshSolver.der¶: Derived data container as a dict.

class solvcon.solver._MethodList¶

A custom list that provides a decorator for keeping names of functions.

>>> mmnames = _MethodList()
>>> @mmnames.register
... def func_of_a_name():
...     pass
>>> mmnames
['func_of_a_name']

This class is a private helper and should only be used in solvcon.solver.

Parallel Computing¶

For distributed-memory parallel computing (i.e., MPI runs), the member svrn indicates the serial number (0-based) the object is. The value of svrn comes from blk. Another member, nsvr, is the total number of collaborative solvers in the parallel run, and is initialized to None.

`solvcon.boundcond`¶

class solvcon.boundcond.BC(bc=None, fpdtype=None)¶

Generic boundary condition abstract class; the base class that all boundary condition classes should subclass.

FIXME: provide doctests as examples.

facn = None¶: An numpy.ndarray as a list of faces. First column is the face index in block. The second column is the face index in bndfcs. The third column is the face index of the related block (if exists).

value = None¶: An numpy.ndarray for attached (specified) value for each boundary face.

nvalue¶

Return the length of vnames as number of values per boundary face. It should be equivalent to the second shape element of value.

FIXME: provide doctests.

__len__()¶

Return the first shape element of facn.

FIXME: provide doctests.

cloneTo(another)¶

Parameters:	another (solvcon.boundcond.BC) – Another BC object.
Returns:	Nothing.

Clone self to another BC object.

create_bcd()¶

Returns:	An object contains the `sc_bound_t` variable for C interfacing.
Return type:	`solvcon.mesh.Bound`

The following code shows how and when to use this method:

>>> import numpy as np
>>> # craft some face numbers for testing.
>>> bndfcs = [0,1,2]
>>> # craft the BC object for testing.
>>> bc = BC()
>>> bc.name = 'some_name'
>>> bc.facn = np.empty((len(bndfcs), BC.BFREL), dtype='int32')
>>> bc.facn.fill(-1)
>>> bc.facn[:,0] = bndfcs
>>> bc.sern = 0
>>> bc.blk = None # should be set to a block.
>>> # test for this method.
>>> bcd = bc.create_bcd()

BC.vnames = []¶: Settable value names.

BC.vdefaults = {}¶: Default values.

Low-Level Interface to C for `BC`¶

sc_bound_t¶

This struct contains essential information for a solvcon.boundcond.BC object in C.

int nbound¶: Number of boundary faces. It’s equivalent to what __len__() returns.

int nvalue¶: Number of values per boundary face.

int *facn¶: Pointer to the data storage of BC.facn.

double *value¶: Pointer to the data storage of BC.value.

class solvcon.mesh.Bound¶

This class associates the C functions for mesh operations to the mesh data and exposes the functions to Python.

_bcd¶: This attribute holds a C struct sc_bound_t for internal use.

setup_bound(bc)¶

Parameters:	bc – The `BC` object that supplies information.

`solvcon.hook`¶

MeshHook performs custom operations at certain pre-defined stages.

class solvcon.hook.MeshHook(cse, **kw)¶: Base type for hooks needing a MeshCase.

`solvcon.anchor`¶

class solvcon.anchor.MeshAnchor(svr, **kw)¶

Callback class to be invoked by MeshSolver objects at various stages.

svr¶: The associated MeshSolver instance.

class solvcon.anchor.MeshAnchorList(svr, *args, **kw)¶

Sequence container for MeshAnchor instances.

svr¶: The associated MeshSolver instance.

References¶

ustmesh_2d_sample.geo¶

/*
 * A Gmsh template file for a rectangle domain.
 */
lc = 0.1;
// vertices.
Point(1) = {4,1,0,lc};
Point(2) = {2,2,0,lc};
Point(3) = {0,1,0,lc};
Point(4) = {0,0,0,lc};
Point(5) = {4,0,0,lc};
Point(6) = {3.5,  1,0,lc};
Point(7) = {  3,  1,0,lc};
Point(8) = {  3,0.5,0,lc};
Point(9) = {3.5,0.5,0,lc};
Point(10) = {  1,  1,0,lc};
Point(11) = {0.5,  1,0,lc};
Point(12) = {0.5,0.5,0,lc};
Point(13) = {  1,0.5,0,lc};
Point(14) = {  2,0.8,0,lc};
Point(15) = {1.5,0.2,0,lc};
Point(16) = {2.5,0.2,0,lc};
// lines.
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,5};
Line(5) = {5,1};
Line(6) = {6,7};
Line(7) = {7,8};
Line(8) = {8,9};
Line(9) = {9,6};
Line(10) = {10,11};
Line(11) = {11,12};
Line(12) = {12,13};
Line(13) = {13,10};
Line(14) = {14,15};
Line(15) = {15,16};
Line(16) = {16,14};
// surface.
Line Loop(1) = {1,2,3,4,5};
Line Loop(2) = {6,7,8,9};
Line Loop(3) = {10,11,12,13};
Line Loop(4) = {14,15,16};
Plane Surface(1) = {1,2,3,4};
// physics.
Physical Line("upper") = {1,2};
Physical Line("left") = {3};
Physical Line("lower") = {4};
Physical Line("right") = {5};
Physical Line("rwin") = {6,7,8,9};
Physical Line("lwin") = {10,11,12,13};
Physical Line("cwin") = {14,15,16};
Physical Surface("domain") = {1};
// vim: set ai et nu ff=unix ft=c:

The following command generate the mesh:

gmsh ustmesh_2d_sample.geo -3

The following command converts the mesh to a VTK file for ParaView:

scg mesh ustmesh_2d_sample.msh ustmesh_2d_sample.vtk

Footnotes

[1]	SOLVCON focuses on two- and three-dimensional meshes. But if we put an additional constraint on the mesh elements: Requiring them to be simplices, it wouldn’t be difficult to extend the data structure of SOLVCON meshes into higher-dimensional space.

Input and Output Facilities¶

`solvcon.io.gmsh`¶

class solvcon.io.gmsh.Gmsh(stream, load=False)¶

Gmsh mesh object. Indices nodes and elements in Gmsh is 1-based (Fortran convention), but 0-based (C convention) indices are used throughout SOLVCON. However, physics groups are using 0-based index.

__init__(stream, load=False)¶

>>> # sample data.
>>> import StringIO
>>> data = """$MeshFormat
... 2.2 0 8
... $EndMeshFormat
... $Nodes
... 3
... 1 -1 0 0
... 2 1 0 0
... 3 0 1 0
... $EndNodes
... $Elements
... 1
... 1 2 2 1 22 1 2 3
... $EndElements
... $PhysicalNames
... 1
... 2 1 "lower"
... $EndPhysicalNames
... $Periodic
... 1
... 0 1 3
... 1
... 1 3
... $EndPeriodic"""

Creation of the object doesn’t load data:

>>> gmsh = Gmsh(StringIO.StringIO(data))
>>> None is gmsh.ndim
True
>>> gmsh.load()
>>> gmsh.ndim
2
>>> gmsh.stream.close() # it's a good habit :-)

We can request to load data on creation by setting load=True. Note the stream will be closed after creation+loading. The default behavior is different to load().

>>> gmsh = Gmsh(StringIO.StringIO(data), load=True)
>>> gmsh.ndim
2
>>> gmsh.stream.closed
True

nnode¶: Number of nodes that is useful for SOLVCON.

ncell¶: Number of cells that is useful for SOLVCON and interior.

load(close=False)¶

Load mesh data from storage.

>>> # sample data.
>>> import StringIO
>>> data = """$MeshFormat
... 2.2 0 8
... $EndMeshFormat
... $Nodes
... 3
... 1 -1 0 0
... 2 1 0 0
... 3 0 1 0
... $EndNodes
... $Elements
... 1
... 1 2 2 1 22 1 2 3
... $EndElements
... $PhysicalNames
... 1
... 2 1 "lower"
... $EndPhysicalNames
... $Periodic
... 1
... 0 1 3
... 1
... 1 3
... $EndPeriodic"""

Load the mesh data after creation of the object. Note the stream is left opened after loading.

>>> stream = StringIO.StringIO(data)
>>> gmsh = Gmsh(stream)
>>> gmsh.load()
>>> stream.closed
False
>>> stream.close() # it's a good habit :-)

We can ask load() to close the stream after loading by using close=True:

>>> gmsh = Gmsh(StringIO.StringIO(data))
>>> gmsh.load(close=True)
>>> gmsh.stream.closed
True

Private Helpers for Loading and Parsing the Mesh File

These private methods are documented for demonstrating the data structure of the loaded meshes. Do not rely on their implementation.

static _check_meta(stream)¶

Load and check the meta data of the mesh. It doesn’t return anything to be stored.

>>> import StringIO
>>> stream = StringIO.StringIO("""$MeshFormat
... 2.2 0 8
... $EndMeshFormat""")
>>> stream.readline()
'$MeshFormat\n'
>>> Gmsh._check_meta(stream)
{}
>>> stream.readline()
''

static _load_nodes(stream)¶

Load node coordinates of the mesh data. Because of the internal data structure of Python, Numpy, and SOLVCON, the loaded nodes are using the 0-based index.

>>> import StringIO
>>> stream = StringIO.StringIO("""$Nodes
... 3
... 1 -1 0 0
... 2 1 0 0
... 3 0 1 0
... $EndNodes""") # a triangle.
>>> stream.readline()
'$Nodes\n'
>>> Gmsh._load_nodes(stream) 
{'nodes': array([[-1.,  0.,  0.], [ 1.,  0.,  0.], [ 0.,  1.,  0.]])}
>>> stream.readline()
''

classmethod _load_elements(stream, nodes)¶

Load element definition of the mesh data. The node indices defined for each element are still 1-based. It returns cltpn, eldim, elems, elgeo, elgrp, ndim, ndmap, and usnds for storage.

>>> from numpy import array
>>> nodes = array([[-1.,  0.,  0.], [ 1.,  0.,  0.], [ 0.,  1.,  0.]])
>>> import StringIO
>>> stream = StringIO.StringIO("""$Elements
... 1
... 1 2 2 1 22 1 2 3
... $EndElements""") # a triangle.
>>> stream.readline()
'$Elements\n'
>>> sorted(Gmsh._load_elements(
...     stream, nodes).items()) 
[('cltpn', array([3], dtype=int32)),
 ('eldim', array([2], dtype=int32)),
 ('elems', array([[ 3,  1,  2,  3, -1, -1, -1, -1, -1]], dtype=int32)),
 ('elgeo', array([22], dtype=int32)),
 ('elgrp', array([1], dtype=int32)),
 ('ndim', 2),
 ('ndmap', array([0, 1, 2], dtype=int32)),
 ('usnds', array([0, 1, 2], dtype=int32))]
>>> stream.readline()
''

static _load_physics(stream)¶

Load physics groups of the mesh data. Return physics for storage.

>>> import StringIO
>>> stream = StringIO.StringIO("""$PhysicalNames
... 1
... 2 1 "lower"
... $EndPhysicalNames""")
>>> stream.readline()
'$PhysicalNames\n'
>>> Gmsh._load_physics(stream)
{'physics': ['1', '2 1 "lower"']}
>>> stream.readline()
''

static _load_periodic(stream)¶

Load periodic definition of the mesh data. Return periodics for storage.

>>> import StringIO
>>> stream = StringIO.StringIO("""$Periodic
... 1
... 0 1 3
... 1
... 1 3
... $EndPeriodic""") # a triangle.
>>> stream.readline()
'$Periodic\n'
>>> Gmsh._load_periodic(stream) 
{'periodics': [{'ndim': 0,
                'stag': 1,
                'nodes': array([[1, 3]], dtype=int32),
                'mtag': 3}]}
>>> stream.readline()
''

_parse_physics()¶: Parse physics groups of the mesh data. Process physics and stores intels.

Mesh Definition and Data Attributes

ELMAP¶: Element definition map. The key is Gmsh element type ID. The value is a 4-tuple: (i) dimension, (ii) number of total nodes, (iii) SOLVCON cell type ID, and (iv) SOLVCON cell node ordering.

stream¶: Input stream (file) of the mesh data.

ndim¶: Number of dimension of this mesh (py:class:int). Stored by _load_elements().

nodes¶: Three-dimensional coordinates of all nodes of Gmsh nodes, 3). Note for even two-dimensional meshes the array still stores three-dimensional coordinates. Stored by _load_nodes().

usnds¶: Indices (0-based) of the nodes really useful for SOLVCON (numpy.ndarray). Stored by _load_elements().

ndmap¶: A mapping array from Gmsh node indices (0-based) to SOLVCON node indices (0-based) (numpy.ndarray). Stored by _load_elements().

cltpn¶: SOLVCON cell type ID for each Gmsh element (py:class:numpy.ndarray). Stored by _load_elements().

elgrp¶: Physics group number of each Gmsh element; the first tag (numpy.ndarray). Stored by _load_elements().

elgeo¶: Geometrical gropu number of each Gmsh element; the second tag (numpy.ndarray). Stored by _load_elements().

eldim¶: Dimension of each Gmsh element (numpy.ndarray). Stored by _load_elements().

elems¶: Gmsh node indices (1-based) of each Gmsh element (numpy.ndarray). Stored by _load_elements().

intels¶: Indices (0-based) of the elements inside the domain (numpy.ndarray). Stored by _parse_physics().

physics¶: Physics groups as a list of 3-tuples: (i) dimension, (ii) index (0-based), and (iii) name. If a physics group has the same dimension as the mesh, it is an interior group. Otherwise, the physics group must have one less dimension than the mesh, and it must be used as the boundary definition. Stored by _load_physics() and then processed by _parse_physics().

periodics¶: Periodic relation list. Each item is a dict:. Stored by _load_periodic().

class solvcon.io.gmsh.GmshIO¶

Proxy to Gmsh file format.

Inheritance diagram of GmshIO

load(stream, bcrej=None, bcmapper=None)¶: Load block from stream with BC mapper applied.

Gas Dynamics Parcel (Under Development)¶

The Euler Equations (Under Development)¶

This document shows the Euler equations of the first-order hyperbolic form.

Reflection of Oblique Shock Wave¶

This example solves a reflecting oblique shock wave, as shown in Figure 6. The system consists of two oblique shock waves, which separate the flow into three zones. The incident shock results from a wedge. The second reflects from a plane wall. Flow properties in all the three zones can be calculated with the following data:

The upstream (zone 1) Mach number $M_1$ and the flow properties density, pressure, and temperature.
The first oblique shock angle $\beta_1$ (between zone 1 and 2) or the flow deflection angle $\theta$ (across zone 1/2 and zone 2/3). Only one of the angle is needed. The other one can be calculated from the given one and $M_1$ . The calculation detail is in ObliqueShockRelation.calc_flow_angle() and ObliqueShockRelation.calc_shock_angle().

Figure 6: Oblique shock reflected from a wall

$M_{1,2,3}$ are the Mach number in the corresponding zone 1, 2, and 3. $\theta$ is the flow deflection angle. $\beta_{1,2}$ are the oblique shock angle behind the first and the second zone, respectively.

SOLVCON will be set up to solve this problem, and the simulated results will be compared with the analytical solution.

Figure 7: Simulated density (non-dimensionalized).

Driving Script¶

SOLVCON uses a driving script to control the numerical simulation. Its general layout is:

#!/usr/bin/env python2.7
# The shebang above directs the operating system to look for a correct
# program to run this script.
#
# We may provide additional information here.


# Import necessary modules.
import os # Python standard library
import numpy as np # http://www.numpy.org
import solvcon as sc # SOLVCON
from solvcon.parcel import gas # A specific SOLVCON solver package we'll use


# ...
# ... other code ...
# ...


# At the end of the file.
if __name__ == '__main__':
    sc.go()

Every driving script has the following lines at the end of the file:

if __name__ == '__main__':
    sc.go()

The if __name__ == '__main__': is a magical Python construct. It will detect that the file is run as a script, not imported as a library (module). Once the detection is evaluated as true, the script call a common execution flow defined in solvcon.go(), which uses the content of the driving script to perform the calculation.

Of course, the file has a lot of other code to set up and configure the calculation, as we’ll describe later. It’s important to note that a driving script is a valid Python program file. The Python language is good for specifying parameters the calculation needs, and as a platform to conduct useful operations much more complex than settings. Any Python module can be imported for use.

See Full Listing of the Driving Script for the driving script of this example: $SCSRC/examples/gas/obrf/go. SOLVCON separates apart the configuration and the execution of a simulation case. The separation is necessary for distributed-memory parallel computing (e.g., MPI). Everything run in the driving script is about the configuration. The execution is conducted by code hidden from users.

To run the simulation, go to the example directory and execute the driving script with the command run and the simulation arrangement name obrf:

$ ./go run obrf

The driving script will then run and print messages:

********************************************************************************
*** Start init (level 0) obrf ... 
*** ****************************************************************************
*** ****************************************************************************
*** *** Start build_domain ... 
*** *** ************************************************************************
*** *** mesh file: None
*** *** ************************************************************************
*** *** *** Start create_block ... 
*** *** *** ********************************************************************
*** *** *** ********************************************************************
*** *** *** End create_block . Elapsed time (sec) = 0.0925829
*** *** ************************************************************************
*** *** ************************************************************************
*** *** End build_domain . Elapsed time (sec) = 0.092721
*** ****************************************************************************
*** ****************************************************************************
*** End init obrf . Elapsed time (sec) = 0.0942369
********************************************************************************
********************************************************************************
*** Start run ... 
*** ****************************************************************************
*** ****************************************************************************
*** *** Start run_provide ... 
*** *** ************************************************************************
*** *** ************************************************************************
*** *** End run_provide . Elapsed time (sec) = 0.000499964
*** ****************************************************************************
*** ****************************************************************************
*** *** Start run_preloop ... 
*** *** ************************************************************************
*** *** Relation of reflected oblique shock:
*** *** - theta = 10.00 deg (flow angle)
*** *** - beta1 = 27.38 deg (shock angle)
*** *** - beta1 = 31.80 deg (shock angle)
*** *** - mach, rho, p, T, a (1) =  3.000  1.000  1.000  1.000  1.183
*** *** - mach, rho, p, T, a (2) =  2.505  1.655  2.054  1.242  1.318
*** *** - mach, rho, p, T, a (3) =  2.090  2.565  3.833  1.494  1.446
*** *** Steps 0/600
*** *** Block information:
*** ***   [Block (2D/centroid): 565 nodes, 1592 faces (100 BC), 1028 cells]
*** *** ************************************************************************
*** *** End run_preloop . Elapsed time (sec) = 0.014756
*** ****************************************************************************
*** 
*** ****************************************************************************
*** *** Start run_march ... 
*** *** ************************************************************************
*** *** ##################################################
*** *** Step 200/600, 0.7s elapsed, 1.4s left
*** *** CFL = 0.11/0.95 - 0.11/0.95 adjusted: 0/0/0
*** *** Performance of obrf:
*** ***   0.696783 seconds in marching solver.
*** ***   0.00348392 seconds/step.
*** ***   3.38902 microseconds/cell.
*** ***   0.29507 Mcells/seconds.
*** ***   1.18028 Mvariables/seconds.
*** *** ##################################################
*** *** Step 400/600, 1.4s elapsed, 0.7s left
*** *** CFL = 0.42/0.95 - 0.11/0.95 adjusted: 0/0/0
*** *** Performance of obrf:
*** ***   1.35721 seconds in marching solver.
*** ***   0.00339303 seconds/step.
*** ***   3.30061 microseconds/cell.
*** ***   0.302974 Mcells/seconds.
*** ***   1.2119 Mvariables/seconds.
*** *** ##################################################
*** *** Step 600/600, 2.1s elapsed, 0.0s left
*** *** CFL = 0.47/0.95 - 0.11/0.95 adjusted: 0/0/0
*** *** ************************************************************************
*** *** End run_march . Elapsed time (sec) = 2.06248
*** ****************************************************************************
*** 
*** ****************************************************************************
*** *** Start run_postloop ... 
*** *** ************************************************************************
*** *** Probe result at Pt/poi#611(3.79306,0.358565,0)601:
*** *** - mach3 = 2.074/2.090 (error=%0.79)
*** *** - rho3  = 2.543/2.565 (error=%0.86)
*** *** - p3    = 3.824/3.833 (error=%0.23)
*** *** Performance of obrf:
*** ***   2.02795 seconds in marching solver.
*** ***   0.00337992 seconds/step.
*** ***   3.28786 microseconds/cell.
*** ***   0.30415 Mcells/seconds.
*** ***   1.2166 Mvariables/seconds.
*** *** Averaged maximum CFL = 0.945858.
*** *** Relation of reflected oblique shock:
*** *** - theta = 10.00 deg (flow angle)
*** *** - beta1 = 27.38 deg (shock angle)
*** *** - beta1 = 31.80 deg (shock angle)
*** *** - mach, rho, p, T, a (1) =  3.000  1.000  1.000  1.000  1.183
*** *** - mach, rho, p, T, a (2) =  2.505  1.655  2.054  1.242  1.318
*** *** - mach, rho, p, T, a (3) =  2.090  2.565  3.833  1.494  1.446
*** *** ************************************************************************
*** *** End run_postloop . Elapsed time (sec) = 0.00133896
*** ****************************************************************************
*** ****************************************************************************
*** *** Start run_exhaust ... 
*** *** ************************************************************************
*** *** ************************************************************************
*** *** End run_exhaust . Elapsed time (sec) = 7.51019e-05
*** ****************************************************************************
*** ****************************************************************************
*** *** Start run_final ... 
*** *** ************************************************************************
*** *** ************************************************************************
*** *** End run_final . Elapsed time (sec) = 9.20296e-05
*** ****************************************************************************
*** ****************************************************************************
*** End run obrf . Elapsed time (sec) = 2.07972
********************************************************************************

Data will be output in directory result/.

Arrangement¶

An arrangement sits at the center of a driving script. It’s nothing more than a decorated Python function with a specific signature. The following function obrf() is the main arrangement we’ll use for the shock reflection problem:

def obrf(casename, **kw):
    return obrf_base(
        # Required positional argument for the name of the simulation case.
        casename,
        # Arguments to the base configuration.
        ssteps=200, psteps=4, edgelength=0.1,
        gamma=1.4, density=1.0, pressure=1.0, mach=3.0, theta=10.0/180*np.pi,
        # Arguments to GasCase.
        time_increment=7.e-3, steps_run=600, **kw)

It’s typical for the arrangement function obrf() to be a thin wrapper which calls another function (in this case, obrf_base()). It should be noted that an arrangement function must take one and only one positional argument: casename. All the other arguments need to be keyword.

To make the function obrf() discoverable by SOLVCON, it needs to be registered with the decorator gas.register_arrangement (gas was imported at the beginning of the driving script):

@gas.register_arrangement
def obrf(casename, **kw):
    # ... contents ...

The function obrf_base() does the real work of configuration:

def obrf_base(
    casename=None, psteps=None, ssteps=None, edgelength=None,
    gamma=None, density=None, pressure=None, mach=None, theta=None, **kw):
    """
    Base configuration of the simulation and return the case object.

    :return: The created Case object.
    :rtype: solvcon.parcel.gas.GasCase
    """

    ############################################################################
    # Step 1: Obtain the analytical solution.
    ############################################################################

    # Calculate the flow properties in all zones separated by the shock.
    relation = ObliqueShockReflection(gamma=gamma, theta=theta, mach1=mach,
                                      rho1=density, p1=pressure, T1=1)

    ############################################################################
    # Step 2: Instantiate the simulation case.
    ############################################################################

    # Create the mesh generator.  Keep it for later use.
    mesher = RectangleMesher(lowerleft=(0,0), upperright=(4,1),
                             edgelength=edgelength)
    # Set up case.
    cse = gas.GasCase(
        # Mesh generator.
        mesher=mesher,
        # Mapping boundary-condition treatments.
        bcmap=generate_bcmap(relation),
        # Use the case name to be the basename for all generated files.
        basefn=casename,
        # Use `cwd`/result to store all generated files.
        basedir=os.path.abspath(os.path.join(os.getcwd(), 'result')),
        # Debug and capture-all.
        debug=False, **kw)

    ############################################################################
    # Step 3: Set up delayed callbacks.
    ############################################################################

    # Field initialization and derived calculations.
    cse.defer(gas.FillAnchor, mappers={'soln': gas.GasSolver.ALMOST_ZERO,
                                       'dsoln': 0.0, 'amsca': gamma})
    cse.defer(gas.DensityInitAnchor, rho=density)
    cse.defer(gas.PhysicsAnchor, rsteps=ssteps)
    # Report information while calculating.
    cse.defer(relation.hookcls)
    cse.defer(gas.ProgressHook, linewidth=ssteps/psteps, psteps=psteps)
    cse.defer(gas.CflHook, fullstop=False, cflmax=10.0, psteps=ssteps)
    cse.defer(gas.MeshInfoHook, psteps=ssteps)
    cse.defer(ReflectionProbe, rect=mesher, relation=relation, psteps=ssteps)
    # Store data.
    cse.defer(gas.PMarchSave,
              anames=[('soln', False, -4),
                      ('rho', True, 0),
                      ('p', True, 0),
                      ('T', True, 0),
                      ('ke', True, 0),
                      ('M', True, 0),
                      ('sch', True, 0),
                      ('v', True, 0.5)],
              psteps=ssteps)

    ############################################################################
    # Final: Return the configured simulation case.
    ############################################################################
    return cse

There are three steps:

Obtain the Analytical Solution to set up all quantities for the simulation.
Instantiate the simulation case object (of type GasCase). The GasCase object needs to know how to set up the mesh (see Mesh Generation) and the boundary-condition (BC) treatment (see BC Treatment Mapping). Section Case Instantiation will explain the details.
Configure callbacks for delayed operations by calling defer() of the constructed simulation GasCase object. Section Callback Configuration will explain these callbacks.

At the end of the base function, the constructed and configured GasCase object is returned.

Although the example has only one arrangement, it’s actually encouraged to have multiple arrangements in a script. In this way one driving script can perform simulations of different parameters or different kinds. Conventionally we place the arrangement functions near the end of the driving script, and the decorated functions (e.g., obrf()) are placed after the base (e.g., obrf_base()). The ordering will make the file easier to read.

Analytical Solution¶

To set up the numerical simulation for the shock-reflection problem, we’ll use class ObliqueShockRelation to calculate necessary parameters by creating a subclass of it:

class ObliqueShockReflection(gas.ObliqueShockRelation):
    def __init__(self, gamma, theta, mach1, rho1, p1, T1):
        super(ObliqueShockReflection, self).__init__(gamma=gamma)
        # Angles and Mach numbers.
        self.theta = theta
        self.mach1 = mach1
        self.beta1 = beta1 = self.calc_shock_angle(mach1, theta)
        self.mach2 = mach2 = self.calc_dmach(mach1, beta1)
        self.beta2 = beta2 = self.calc_shock_angle(mach2, theta)
        self.mach3 = mach3 = self.calc_dmach(mach2, beta2)
        # Flow properties in the first zone.
        self.rho1 = rho1
        self.p1 = p1
        self.T1 = T1
        self.a1 = np.sqrt(gamma*p1/rho1)
        # Flow properties in the second zone.
        self.rho2 = rho2 = rho1 * self.calc_density_ratio(mach1, beta1)
        self.p2 = p2 = p1 * self.calc_pressure_ratio(mach1, beta1)
        self.T2 = T2 = T1 * self.calc_temperature_ratio(mach1, beta1)
        self.a2 = np.sqrt(gamma*p2/rho2)
        # Flow properties in the third zone.
        self.rho3 = rho3 = rho2 * self.calc_density_ratio(mach2, beta2)
        self.p3 = p3 = p2 * self.calc_pressure_ratio(mach2, beta2)
        self.T3 = T3 = T2 * self.calc_temperature_ratio(mach2, beta2)
        self.a3 = np.sqrt(gamma*p3/rho3)

    def __str__(self):
        msg = 'Relation of reflected oblique shock:\n'
        msg += '- theta = %5.2f deg (flow angle)\n' % (self.theta/np.pi*180)
        msg += '- beta1 = %5.2f deg (shock angle)\n' % (self.beta1/np.pi*180)
        msg += '- beta1 = %5.2f deg (shock angle)\n' % (self.beta2/np.pi*180)
        def property_string(zone):
            values = [getattr(self, '%s%d' % (key, zone))
                      for key in 'mach', 'rho', 'p', 'T', 'a']
            messages = [' %6.3f' % val for val in values]
            return ''.join(messages)
        msg += '- mach, rho, p, T, a (1) =' + property_string(1) + '\n'
        msg += '- mach, rho, p, T, a (2) =' + property_string(2) + '\n'
        msg += '- mach, rho, p, T, a (3) =' + property_string(3)
        return msg

    @property
    def hookcls(self):
        relation = self
        class _ShowRelation(sc.MeshHook):
            def preloop(self):
                for msg in str(relation).split('\n'):
                    self.info(msg + '\n')
            postloop = preloop
        return _ShowRelation

For the detail of ObliqueShockRelation, see Oblique Shock Relation.

Case Instantiation¶

An instance of GasCase represents a numerical simulation using the gas module. In addition to Mesh Generation and BC Treatment Mapping, other miscellaneous settings can be supplied through the GasCase constructor.

Mesh Generation¶

An unstructured mesh is required for a SOLVCON simulation. A mesh file can be created beforehand or on-the-fly with the simulation. The example uses the latter approach. The following is an example of mesh generating function that calls Gmsh:

class RectangleMesher(object):
    """
    Representation of a rectangle and the Gmsh meshing helper.

    :ivar lowerleft: (x0, y0) of the rectangle.
    :type lowerleft: tuple
    :ivar upperright: (x1, y1) of the rectangle.
    :type upperright: tuple
    :ivar edgelength: Length of each mesh cell edge.
    :type edgelength: float
    """

    GMSH_SCRIPT = """
    // vertices.
    Point(1) = {%(x1)g,%(y1)g,0,%(edgelength)g};
    Point(2) = {%(x0)g,%(y1)g,0,%(edgelength)g};
    Point(3) = {%(x0)g,%(y0)g,0,%(edgelength)g};
    Point(4) = {%(x1)g,%(y0)g,0,%(edgelength)g};
    // lines.
    Line(1) = {1,2};
    Line(2) = {2,3};
    Line(3) = {3,4};
    Line(4) = {4,1};
    // surface.
    Line Loop(1) = {1,2,3,4};
    Plane Surface(1) = {1};
    // physics.
    Physical Line("upper") = {1};
    Physical Line("left") = {2};
    Physical Line("lower") = {3};
    Physical Line("right") = {4};
    Physical Surface("domain") = {1};
    """.strip()

    def __init__(self, lowerleft, upperright, edgelength):
        assert 2 == len(lowerleft)
        self.lowerleft = tuple(float(val) for val in lowerleft)
        assert 2 == len(upperright)
        self.upperright = tuple(float(val) for val in upperright)
        self.edgelength = float(edgelength)

    def __call__(self, cse):
        x0, y0 = self.lowerleft
        x1, y1 = self.upperright
        script_arguments = dict(
            edgelength=self.edgelength, x0=x0, y0=y0, x1=x1, y1=y1)
        cmds = self.GMSH_SCRIPT % script_arguments
        cmds = [cmd.rstrip() for cmd in cmds.strip().split('\n')]
        gmh = sc.Gmsh(cmds)()
        return gmh.toblock(bcname_mapper=cse.condition.bcmap)

BC Treatment Mapping¶

Boundary-condition treatments are specified by creating a dict to map the name of the boundary to a specific BC class.

def generate_bcmap(relation):
    # Flow properties (fp).
    fpb = {
        'gamma': relation.gamma, 'rho': relation.rho1,
        'v2': 0.0, 'v3': 0.0, 'p': relation.p1,
    }
    fpb['v1'] = relation.mach1*np.sqrt(relation.gamma*fpb['p']/fpb['rho'])
    fpt = fpb.copy()
    fpt['rho'] = relation.rho2
    fpt['p'] = relation.p2
    V2 = relation.mach2 * relation.a2
    fpt['v1'] = V2 * np.cos(relation.theta)
    fpt['v2'] = -V2 * np.sin(relation.theta)
    fpi = fpb.copy()
    # BC map.
    bcmap = {
        'upper': (sc.bctregy.GasInlet, fpt,),
        'left': (sc.bctregy.GasInlet, fpb,),
        'right': (sc.bctregy.GasNonrefl, {},),
        'lower': (sc.bctregy.GasWall, {},),
    }
    return bcmap

Callback Configuration¶

SOLVCON provides general-purpose, application-agnostic solving facilities. To describe the problem to SOLVCON, we need to provide both data (numbers) and logic (computer code) in the driving script. The supplied code will be called back at proper points while the simulation is running.

Classes MeshHook and MeshAnchor are the fundamental constructs to make callbacks in the sequential and parallel runtime environment, respectively. The module gas includes useful callbacks, but we still need to write a couple of them in the driving script.

The shock reflection problem uses three categories of callbacks.

Initialization and calculation:

FillAnchor

DensityInitAnchor

PhysicsAnchor

Reporting:

ObliqueShockReflection.hookcls()

ProgressHook

CflHook

MeshInfoHook

ReflectionProbe

Output:

PMarchSave

The order of these callbacks is important. Dependency between callbacks is allowed.

View Results¶

After simulation, the results are stored in directory result/ as VTK unstructured data files that can be opened and processed by using ParaView. The result in Figure 7 was produced in this way. Other quantities can also be visualized, e.g., the Mach number shown in Figure 8.

Figure 8: Mach number at the final time step of the arrangement obrf_fine.

Both of Figures 7 and 8 are obtained with the arrangement obrf_fine.

Oblique Shock Relation¶

An oblique shock results from a sudden change of direction of supersonic flow. The relations of density ( $\rho$ ), pressure ( $p$ ), and temprature ( $T$ ) across the shock can be obtained analytically [Anderson03]. In addition, two angles are defined:

The angle of the oblique shock wave deflected from the upstream is $\beta$ ; the shock angle.
The angle of the flow behind the shock wave deflected from the upstream is $\theta$ ; the flow angle.

See Figure 9 for the illustration of the two angles.

Figure 9: Oblique shock wave by a wedge

$M$ is Mach number. $\theta$ is the flow deflection angle. $\beta$ is the oblique shock angle.

Methods of calculating the shock relations are organized in the class ObliqueShockRelation. To obtain the relations of density ( $\rho$ ), pressure ( $p$ ), and temprature ( $T$ ), the control volume across the shock is emplyed, as shown in Figure 10. In the figure and in ObliqueShockRelation, subscript 1 denotes upstream properties and subscript 2 denotes downstream properties. Derivation of the relation uses a rotated coordinate system $(n, t)$ defined by the oblique shock, where $\hat{n}$ is the unit vector normal to the shock, and $\hat{t}$ is the unit vector tangential to the shock. But in this document we won’t go into the detail.

Figure 10: Properties across an oblique shock

The flow properties in the upstream zone of the oblique shock are $v_1, M_1, \rho_1, p_1, T_1$ . Those in the downstream zone of the shock are $v_2, M_2, \rho_2, p_2, T_2$ .

class solvcon.parcel.gas.oblique_shock.ObliqueShockRelation(gamma)¶

Calculators of oblique shock relations.

The constructor must take the ratio of specific heat:

>>> ObliqueShockRelation()
Traceback (most recent call last):
    ...
TypeError: __init__() takes exactly 2 arguments (1 given)
>>> ob = ObliqueShockRelation(gamma=1.4)

The ratio of specific heat can be accessed through the gamma attribute:

>>> ob.gamma
1.4

The object can be used to calculate shock relations. For example, calc_density_ratio() returns the $\rho_2/\rho_1$ :

>>> round(ob.calc_density_ratio(mach1=3, beta=37.8/180*np.pi), 10)
2.4204302545

The solution changes as gamma changes:

>>> ob.gamma = 1.2
>>> round(ob.calc_density_ratio(mach1=3, beta=37.8/180*np.pi), 10)
2.7793244902

gamma¶: Ratio of specific heat $\gamma$ , dimensionless.

ObliqueShockRelation provides three methods to calculate the ratio of flow properties across the shock. $M_1$ and $\beta$ are required arguments:

With $M_1$ available, the shock angle $\beta$ can be calculated from the flow angle $\theta$ , or vice versa, by using the following two methods:

$\beta$ : ObliqueShockRelation.calc_shock_angle()
$\theta$ : ObliqueShockRelation.calc_flow_angle()

The following method calculates the downstream Mach number, with the upstream Mach number $M_1$ and either of $\beta$ or $\theta$ supplied:

$M_2$ : calc_dmach()

Listing of all methods:

ObliqueShockRelation.calc_density_ratio(mach1, beta)¶

Calculate the ratio of density $\rho$ across an oblique shock wave of which the angle deflected from the upstream flow is $\beta$ and the upstream Mach number is $M_1$ :

$\frac{\rho_2}{\rho_1} = \frac{(\gamma + 1) M_{n1}^2} {(\gamma - 1) M_{n1}^2 + 2}$

where $M_{n1} = M_1\sin\beta$ .

This method accepts scalar:

>>> ob = ObliqueShockRelation(gamma=1.4)
>>> round(ob.calc_density_ratio(mach1=3, beta=37.8/180*np.pi), 10)
2.4204302545

as well as numpy.ndarray:

>>> angle = 37.8/180*np.pi; angle = np.array([angle, angle])
>>> np.round(ob.calc_density_ratio(mach1=3, beta=angle), 10).tolist()
[2.4204302545, 2.4204302545]

Parameters:	mach1 – Upstream Mach number $M_1$ , dimensionless. beta – Oblique shock angle $\beta$ deflected from the upstream flow, in radian.

ObliqueShockRelation.calc_pressure_ratio(mach1, beta)¶

Calculate the ratio of pressure $p$ across an oblique shock wave of which the angle deflected from the upstream flow is $\beta$ and the upstream Mach number is $M_1$ :

$\frac{p_2}{p_1} = 1 + \frac{2\gamma}{\gamma+1}(M_{n1}^2 - 1)$

where $M_{n1} = M_1\sin\beta$ .

This method accepts scalar:

>>> ob = ObliqueShockRelation(gamma=1.4)
>>> round(ob.calc_pressure_ratio(mach1=3, beta=37.8/180*np.pi), 10)
3.7777114257

as well as numpy.ndarray:

>>> angle = 37.8/180*np.pi; angle = np.array([angle, angle])
>>> np.round(ob.calc_pressure_ratio(mach1=3, beta=angle), 10).tolist()
[3.7777114257, 3.7777114257]

Parameters:	mach1 – Upstream Mach number $M_1$ , dimensionless. beta – Oblique shock angle $\beta$ deflected from the upstream flow, in radian.

ObliqueShockRelation.calc_temperature_ratio(mach1, beta)¶

Calculate the ratio of temperature $T$ across an oblique shock wave of which the angle deflected from the upstream flow is $\beta$ and the upstream Mach number is $M_1$ :

$\frac{T_2}{T_1} = \frac{p_2}{p_1} \frac{\rho_1}{\rho_2}$

where both $p_2/p_1$ and $\rho_1/\rho_2$ are functions of $\gamma$ , $M_1$ , and $\beta$ . See also calc_pressure_ratio() and calc_density_ratio().

This method accepts scalar:

>>> ob = ObliqueShockRelation(gamma=1.4)
>>> round(ob.calc_temperature_ratio(mach1=3, beta=37.8/180*np.pi), 10)
1.5607602899

as well as numpy.ndarray:

>>> angle = 37.8/180*np.pi; angle = np.array([angle, angle])
>>> np.round(ob.calc_temperature_ratio(mach1=3, beta=angle), 10).tolist()
[1.5607602899, 1.5607602899]

Parameters:	mach1 – Upstream Mach number $M_1$ , dimensionless. beta – Oblique shock angle $\beta$ deflected from the upstream flow, in radian.

ObliqueShockRelation.calc_dmach(mach1, beta=None, theta=None, delta=1)¶

Calculate the downstream Mach number from the upstream Mach number $M_1$ and either of the shock or flow deflection angles:

$M_2 = \frac{M_{n2}}{\sin(\beta-\theta)}$

where $M_{n2}$ is calculated from calc_normal_dmach() with $M_{n1} = M_1\sin\beta$ .

The method can be invoked with only either $\beta$ or $\theta$ :

>>> ob = ObliqueShockRelation(gamma=1.4)
>>> ob.calc_dmach(3, beta=0.2, theta=0.1)
Traceback (most recent call last):
    ...
ValueError: got (beta=0.2, theta=0.1), but I need to take either beta or theta
>>> ob.calc_dmach(3)
Traceback (most recent call last):
    ...
ValueError: got (beta=None, theta=None), but I need to take either beta or theta

This method accepts scalar:

>>> round(ob.calc_dmach(3, beta=37.8/180*np.pi), 10)
1.9924827009
>>> round(ob.calc_dmach(3, theta=20./180*np.pi), 10)
1.9941316656

as well as numpy.ndarray:

>>> angle = 37.8/180*np.pi; angle = np.array([angle, angle])
>>> np.round(ob.calc_dmach(3, beta=angle), 10).tolist()
[1.9924827009, 1.9924827009]
>>> angle = 20./180*np.pi; angle = np.array([angle, angle])
>>> np.round(ob.calc_dmach(3, theta=angle), 10).tolist()
[1.9941316656, 1.9941316656]

Parameters:

mach1 – Upstream Mach number $M_1$ , dimensionless.
beta – Oblique shock angle $\beta$ deflected from the upstream flow, in radian.
theta – Downstream flow angle $\theta$ deflected from the upstream flow, in radian.
delta – A switching integer $\delta$ . For $\delta = 0$ , the function gives the solution of strong shock, while for $\delta = 1$ , it gives the solution of weak shock. This keyword argument is only valid when theta is given. The default value is 1.

ObliqueShockRelation.calc_normal_dmach(mach_n1)¶

Calculate the downstream Mach number from the given upstream Mach number $M_{n1}$ , in the direction normal to the shock:

$M_{n2} = \sqrt{\frac{(\gamma-1)M_{n1}^2 + 2} {2\gamma M_{n1}^2 - (\gamma-1)}}$

This method accepts scalar:

>>> ob = ObliqueShockRelation(gamma=1.4)
>>> round(ob.calc_normal_dmach(mach_n1=3), 10)
0.4751909633

as well as numpy.ndarray:

>>> np.round(ob.calc_normal_dmach(mach_n1=np.array([3, 3])), 10).tolist()
[0.4751909633, 0.4751909633]

Parameters:	mach_n1 – Upstream Mach number $M_{n1}$ normal to the shock wave, dimensionless.

ObliqueShockRelation.calc_flow_angle(mach1, beta)¶

Calculate the downstream flow angle $\theta$ deflected from the upstream flow by using calc_flow_tangent(), in radian.

This method accepts scalar:

>>> ob = ObliqueShockRelation(gamma=1.4)
>>> angle = 48.25848/180*np.pi
>>> round(ob.calc_flow_angle(mach1=4, beta=angle)/np.pi*180, 10)
32.0000000807

as well as numpy.ndarray:

>>> angle = 48.25848/180*np.pi; angle = np.array([angle, angle])
>>> np.round((ob.calc_flow_angle(mach1=4, beta=angle)/np.pi*180), 10).tolist()
[32.0000000807, 32.0000000807]

See Example 4.6 in [Anderson03] for the forward analysis. The above is the inverse analysis.

Parameters:	mach1 – Upstream Mach number $M_1$ , dimensionless. beta – Oblique shock angle $\beta$ deflected from the upstream flow, in radian.

ObliqueShockRelation.calc_flow_tangent(mach1, beta)¶

Calculate the trigonometric tangent function $\tan\beta$ of the downstream flow angle $\theta$ deflected from the upstream flow by using the $\theta$ - $\beta$ - $M$ relation:

$\tan\theta = 2\cot\beta \frac{M_1^2\sin^2\beta - 1} {M_1^2(\gamma+\cos2\beta) + 2}$

This method accepts scalar:

>>> ob = ObliqueShockRelation(gamma=1.4)
>>> angle = 48.25848/180*np.pi
>>> round(ob.calc_flow_tangent(mach1=4, beta=angle), 10)
0.6248693539

as well as numpy.ndarray:

>>> angle = 48.25848/180*np.pi; angle = np.array([angle, angle])
>>> np.round(ob.calc_flow_tangent(mach1=4, beta=angle), 10).tolist()
[0.6248693539, 0.6248693539]

See Example 4.6 in [Anderson03] for the forward analysis. The above is the inverse analysis.

Parameters:	mach1 – Upstream Mach number $M_1$ , dimensionless. beta – Oblique shock angle $\beta$ deflected from the upstream flow, in radian.

ObliqueShockRelation.calc_shock_angle(mach1, theta, delta=1)¶

Calculate the downstream shock angle $\beta$ deflected from the upstream flow by using calc_shock_tangent(), in radian.

This method accepts scalar:

>>> ob = ObliqueShockRelation(gamma=1.4)
>>> angle = 32./180*np.pi
>>> round(ob.calc_shock_angle(mach1=4, theta=angle, delta=1)/np.pi*180, 10)
48.2584798722

as well as numpy.ndarray:

>>> angle = np.array([angle, angle])
>>> np.round(ob.calc_shock_angle(mach1=4, theta=angle, delta=1)/np.pi*180,
...          10).tolist()
[48.2584798722, 48.2584798722]

See Example 4.6 in [Anderson03] for the analysis.

Parameters:	mach1 – Upstream Mach number $M_1$ , dimensionless. theta – Downstream flow angle $\theta$ deflected from the upstream flow, in radian. delta – A switching integer $\delta$ . For $\delta = 0$ , the function gives the solution of strong shock, while for $\delta = 1$ , it gives the solution of weak shock. The default value is 1.

ObliqueShockRelation.calc_shock_tangent(mach1, theta, delta)¶

Calculate the downstream shock angle $\beta$ deflected from the upstream flow by using the alternative $\beta$ - $\theta$ - $M$ relation:

$\tan\beta = \frac{M_1^2 - 1 + 2\lambda\cos\left(\frac{4\pi\delta + \cos^{-1}\chi}{3}\right)} {3\left(1 + \frac{\gamma-1}{2}M_1^2\right)\tan\theta}$

where $\lambda$ and $\chi$ are obtained internally by calling calc_shock_tangent_aux().

This method accepts scalar:

>>> ob = ObliqueShockRelation(gamma=1.4)
>>> angle = 32./180*np.pi
>>> round(ob.calc_shock_tangent(mach1=4, theta=angle, delta=1), 10)
1.1207391858

as well as numpy.ndarray:

>>> angle = np.array([angle, angle])
>>> np.round(ob.calc_shock_tangent(mach1=4, theta=angle, delta=1),
...          10).tolist()
[1.1207391858, 1.1207391858]

See Example 4.6 in [Anderson03] for the analysis.

Parameters:	mach1 – Upstream Mach number $M_1$ , dimensionless. theta – Downstream flow angle $\theta$ deflected from the upstream flow, in radian. delta – A switching integer $\delta$ . For $\delta = 0$ , the function gives the solution of strong shock, while for $\delta = 1$ , it gives the solution of weak shock.

ObliqueShockRelation.calc_shock_tangent_aux(mach1, theta)¶

Calculate the $\lambda$ and $\chi$ functions used by calc_shock_tangent():

$\lambda = \sqrt{(M_1^2-1)^2 - 3\left(1+\frac{\gamma-1}{2}M_1^2\right) \left(1+\frac{\gamma+1}{2}M_1^2\right)\tan^2\theta}$

$\chi = \frac{(M_1^2-1)^3 - 9\left(1+\frac{\gamma-1}{2}M_1^2\right) \left(1+\frac{\gamma-1}{2}M_1^2+\frac{\gamma+1}{4}M_1^4\right) \tan^2\theta} {\lambda^3}$

This method accepts scalar:

>>> ob = ObliqueShockRelation(gamma=1.4)
>>> lmbd, chi = ob.calc_shock_tangent_aux(mach1=4, theta=32./180*np.pi)
>>> round(lmbd, 10), round(chi, 10)
(11.2080188412, 0.7428957121)

as well as numpy.ndarray:

>>> angle = 32./180*np.pi; angle = np.array([angle, angle])
>>> lmbd, chi = ob.calc_shock_tangent_aux(mach1=4, theta=angle)
>>> np.round(lmbd, 10).tolist()
[11.2080188412, 11.2080188412]
>>> np.round(chi, 10).tolist()
[0.7428957121, 0.7428957121]

See Example 4.6 in [Anderson03] for the analysis.

Parameters:	mach1 – Upstream Mach number $M_1$ , dimensionless. theta – Downstream flow angle $\theta$ deflected from the upstream flow, in radian.

Full Listing of the Driving Script¶

#!/usr/bin/env python2.7
# -*- coding: UTF-8 -*-
#
# Copyright (c) 2014, Yung-Yu Chen <yyc@solvcon.net>
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
# - Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
# - Neither the name of the SOLVCON nor the names of its contributors may be
#   used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.


"""
This is an example for solving the problem of oblique shock reflection.
"""


import os # Python standard library
import numpy as np # http://www.numpy.org
import solvcon as sc # SOLVCON
from solvcon.parcel import gas # A specific SOLVCON solver package we'll use


class ObliqueShockReflection(gas.ObliqueShockRelation):
    def __init__(self, gamma, theta, mach1, rho1, p1, T1):
        super(ObliqueShockReflection, self).__init__(gamma=gamma)
        # Angles and Mach numbers.
        self.theta = theta
        self.mach1 = mach1
        self.beta1 = beta1 = self.calc_shock_angle(mach1, theta)
        self.mach2 = mach2 = self.calc_dmach(mach1, beta1)
        self.beta2 = beta2 = self.calc_shock_angle(mach2, theta)
        self.mach3 = mach3 = self.calc_dmach(mach2, beta2)
        # Flow properties in the first zone.
        self.rho1 = rho1
        self.p1 = p1
        self.T1 = T1
        self.a1 = np.sqrt(gamma*p1/rho1)
        # Flow properties in the second zone.
        self.rho2 = rho2 = rho1 * self.calc_density_ratio(mach1, beta1)
        self.p2 = p2 = p1 * self.calc_pressure_ratio(mach1, beta1)
        self.T2 = T2 = T1 * self.calc_temperature_ratio(mach1, beta1)
        self.a2 = np.sqrt(gamma*p2/rho2)
        # Flow properties in the third zone.
        self.rho3 = rho3 = rho2 * self.calc_density_ratio(mach2, beta2)
        self.p3 = p3 = p2 * self.calc_pressure_ratio(mach2, beta2)
        self.T3 = T3 = T2 * self.calc_temperature_ratio(mach2, beta2)
        self.a3 = np.sqrt(gamma*p3/rho3)

    def __str__(self):
        msg = 'Relation of reflected oblique shock:\n'
        msg += '- theta = %5.2f deg (flow angle)\n' % (self.theta/np.pi*180)
        msg += '- beta1 = %5.2f deg (shock angle)\n' % (self.beta1/np.pi*180)
        msg += '- beta1 = %5.2f deg (shock angle)\n' % (self.beta2/np.pi*180)
        def property_string(zone):
            values = [getattr(self, '%s%d' % (key, zone))
                      for key in 'mach', 'rho', 'p', 'T', 'a']
            messages = [' %6.3f' % val for val in values]
            return ''.join(messages)
        msg += '- mach, rho, p, T, a (1) =' + property_string(1) + '\n'
        msg += '- mach, rho, p, T, a (2) =' + property_string(2) + '\n'
        msg += '- mach, rho, p, T, a (3) =' + property_string(3)
        return msg

    @property
    def hookcls(self):
        relation = self
        class _ShowRelation(sc.MeshHook):
            def preloop(self):
                for msg in str(relation).split('\n'):
                    self.info(msg + '\n')
            postloop = preloop
        return _ShowRelation


class RectangleMesher(object):
    """
    Representation of a rectangle and the Gmsh meshing helper.

    :ivar lowerleft: (x0, y0) of the rectangle.
    :type lowerleft: tuple
    :ivar upperright: (x1, y1) of the rectangle.
    :type upperright: tuple
    :ivar edgelength: Length of each mesh cell edge.
    :type edgelength: float
    """

    GMSH_SCRIPT = """
    // vertices.
    Point(1) = {%(x1)g,%(y1)g,0,%(edgelength)g};
    Point(2) = {%(x0)g,%(y1)g,0,%(edgelength)g};
    Point(3) = {%(x0)g,%(y0)g,0,%(edgelength)g};
    Point(4) = {%(x1)g,%(y0)g,0,%(edgelength)g};
    // lines.
    Line(1) = {1,2};
    Line(2) = {2,3};
    Line(3) = {3,4};
    Line(4) = {4,1};
    // surface.
    Line Loop(1) = {1,2,3,4};
    Plane Surface(1) = {1};
    // physics.
    Physical Line("upper") = {1};
    Physical Line("left") = {2};
    Physical Line("lower") = {3};
    Physical Line("right") = {4};
    Physical Surface("domain") = {1};
    """.strip()

    def __init__(self, lowerleft, upperright, edgelength):
        assert 2 == len(lowerleft)
        self.lowerleft = tuple(float(val) for val in lowerleft)
        assert 2 == len(upperright)
        self.upperright = tuple(float(val) for val in upperright)
        self.edgelength = float(edgelength)

    def __call__(self, cse):
        x0, y0 = self.lowerleft
        x1, y1 = self.upperright
        script_arguments = dict(
            edgelength=self.edgelength, x0=x0, y0=y0, x1=x1, y1=y1)
        cmds = self.GMSH_SCRIPT % script_arguments
        cmds = [cmd.rstrip() for cmd in cmds.strip().split('\n')]
        gmh = sc.Gmsh(cmds)()
        return gmh.toblock(bcname_mapper=cse.condition.bcmap)


def generate_bcmap(relation):
    # Flow properties (fp).
    fpb = {
        'gamma': relation.gamma, 'rho': relation.rho1,
        'v2': 0.0, 'v3': 0.0, 'p': relation.p1,
    }
    fpb['v1'] = relation.mach1*np.sqrt(relation.gamma*fpb['p']/fpb['rho'])
    fpt = fpb.copy()
    fpt['rho'] = relation.rho2
    fpt['p'] = relation.p2
    V2 = relation.mach2 * relation.a2
    fpt['v1'] = V2 * np.cos(relation.theta)
    fpt['v2'] = -V2 * np.sin(relation.theta)
    fpi = fpb.copy()
    # BC map.
    bcmap = {
        'upper': (sc.bctregy.GasInlet, fpt,),
        'left': (sc.bctregy.GasInlet, fpb,),
        'right': (sc.bctregy.GasNonrefl, {},),
        'lower': (sc.bctregy.GasWall, {},),
    }
    return bcmap


class ReflectionProbe(gas.ProbeHook):
    """
    Place a probe for the flow properties in the reflected region.
    """

    def __init__(self, cse, **kw):
        """
        :param relation: Instruct shock relations.
        :type relation: ObliqueShockReflection
        :param rect: Instruct the mesh rectangle.
        :type rect: RectangleMesher
        """
        rect = kw.pop('rect')
        self.relation = relation = kw.pop('relation')
        factor = kw.pop('factor', 0.9)
        # Detemine location.
        theta = relation.theta
        beta1 = relation.beta1
        beta2 = relation.beta2
        x0, y0 = rect.lowerleft
        x1, y1 = rect.upperright
        lgh = (y1-y0) / np.tan(beta1)
        hgt = factor * (x1-x0-lgh) * np.tan((beta2-theta)/2)
        lgh = factor * (x1-x0-lgh) + lgh
        poi = ('poi', x0+lgh, y0+hgt, 0.0)
        # Call super.
        kw['coords'] = (poi,)
        kw['speclst'] = ['M', 'rho', 'p']
        super(ReflectionProbe, self).__init__(cse, **kw)

    def postloop(self):
        super(ReflectionProbe, self).postloop()
        rel = self.relation
        self.info('Probe result at %s:\n' % self.points[0])
        M, rho, p = self.points[0].vals[-1][1:]
        self.info('- mach3 = %.3f/%.3f (error=%%%.2f)\n' % (
            M, rel.mach3, abs((M-rel.mach3)/rel.mach3)*100))
        self.info('- rho3  = %.3f/%.3f (error=%%%.2f)\n' % (
            rho, rel.rho3, abs((rho-rel.rho3)/rel.rho3)*100))
        self.info('- p3    = %.3f/%.3f (error=%%%.2f)\n' % (
            p, rel.p3, abs((p-rel.p3)/rel.p3)*100))


def obrf_base(
    casename=None, psteps=None, ssteps=None, edgelength=None,
    gamma=None, density=None, pressure=None, mach=None, theta=None, **kw):
    """
    Base configuration of the simulation and return the case object.

    :return: The created Case object.
    :rtype: solvcon.parcel.gas.GasCase
    """

    ############################################################################
    # Step 1: Obtain the analytical solution.
    ############################################################################

    # Calculate the flow properties in all zones separated by the shock.
    relation = ObliqueShockReflection(gamma=gamma, theta=theta, mach1=mach,
                                      rho1=density, p1=pressure, T1=1)

    ############################################################################
    # Step 2: Instantiate the simulation case.
    ############################################################################

    # Create the mesh generator.  Keep it for later use.
    mesher = RectangleMesher(lowerleft=(0,0), upperright=(4,1),
                             edgelength=edgelength)
    # Set up case.
    cse = gas.GasCase(
        # Mesh generator.
        mesher=mesher,
        # Mapping boundary-condition treatments.
        bcmap=generate_bcmap(relation),
        # Use the case name to be the basename for all generated files.
        basefn=casename,
        # Use `cwd`/result to store all generated files.
        basedir=os.path.abspath(os.path.join(os.getcwd(), 'result')),
        # Debug and capture-all.
        debug=False, **kw)

    ############################################################################
    # Step 3: Set up delayed callbacks.
    ############################################################################

    # Field initialization and derived calculations.
    cse.defer(gas.FillAnchor, mappers={'soln': gas.GasSolver.ALMOST_ZERO,
                                       'dsoln': 0.0, 'amsca': gamma})
    cse.defer(gas.DensityInitAnchor, rho=density)
    cse.defer(gas.PhysicsAnchor, rsteps=ssteps)
    # Report information while calculating.
    cse.defer(relation.hookcls)
    cse.defer(gas.ProgressHook, linewidth=ssteps/psteps, psteps=psteps)
    cse.defer(gas.CflHook, fullstop=False, cflmax=10.0, psteps=ssteps)
    cse.defer(gas.MeshInfoHook, psteps=ssteps)
    cse.defer(ReflectionProbe, rect=mesher, relation=relation, psteps=ssteps)
    # Store data.
    cse.defer(gas.PMarchSave,
              anames=[('soln', False, -4),
                      ('rho', True, 0),
                      ('p', True, 0),
                      ('T', True, 0),
                      ('ke', True, 0),
                      ('M', True, 0),
                      ('sch', True, 0),
                      ('v', True, 0.5)],
              psteps=ssteps)

    ############################################################################
    # Final: Return the configured simulation case.
    ############################################################################
    return cse


@gas.register_arrangement
def obrf(casename, **kw):
    return obrf_base(
        # Required positional argument for the name of the simulation case.
        casename,
        # Arguments to the base configuration.
        ssteps=200, psteps=4, edgelength=0.1,
        gamma=1.4, density=1.0, pressure=1.0, mach=3.0, theta=10.0/180*np.pi,
        # Arguments to GasCase.
        time_increment=7.e-3, steps_run=600, **kw)


@gas.register_arrangement
def obrf_fine(casename, **kw):
    return obrf_base(
        casename,
        ssteps=200, psteps=4, edgelength=0.02,
        gamma=1.4, density=1.0, pressure=1.0, mach=3.0, theta=10.0/180*np.pi,
        time_increment=1.e-3, steps_run=4000, **kw)


if __name__ == '__main__':
    sc.go()

# vim: set ff=unix fenc=utf8 ft=python ai et sw=4 ts=4 tw=79:

Sod’s Shock Tube (Under Development)¶

A classic example to verify whether a CFD algorithm the Sod shock tube problem [Sod78]. We will introduce this problem in what follows.

In short, a shock tube problem is a Riemann problem with the Euler equations. This is a good benchmark to compare different CFD algorithm results.

The Euler equations consist of conservation of mass (Eq. (1)), of momentum (Eq. (2)), and of energy (Eq. (3)).

(1) $\dpd{\rho}{t} + \dpd{{\rho}{v}}{x} = 0$

(2) $\dpd{\rho{v}}{t} + \dpd{(p+\rho{v^2})}{x} = 0$

(3) $\dpd{}{t}\left(\frac{p}{\gamma-1} + \frac{\rho{v^2}}{2}\right) + \dpd{}{x}\left(\frac{\gamma}{\gamma-1}pv+\frac{1}{2}\rho{v^3}\right) = 0$

By defining

(4) $\bvec{u} = \left(\begin{array}{c} u_1 \\ u_2 \\ u_3 \end{array}\right) \defeq \left(\begin{array}{c} \rho \\ \rho v \\ \rho\left(\frac{1}{\gamma-1}\frac{p}{\rho} + \frac{v^2}{2}\right) \end{array}\right)$

(5) $\bvec{f} = \left(\begin{array}{c} f_1 \\ f_2 \\ f_3 \end{array}\right) \defeq \left(\begin{array}{c} u_2 \\ (\gamma-1)u_3 - \frac{\gamma-3}{2}\frac{u_2^2}{u_1} \\ \gamma\frac{u_2u_3}{u_1} - \frac{\gamma-1}{2}\frac{u_2^3}{u_1^2} \end{array}\right)$

we can rewrite Eqs. (1), (2), and (3) in a general form for nonlinear hyperbolic PDEs:

(6) $\dpd{\bvec{u}}{t} + \dpd{\bvec{f}(\bvec{u})}{x} = 0$

The initial condition of the Riemann problem is defined as:

(7) $\bvec{u} = \left(\begin{array}{c} \rho_L \\ u_L \\ p_L \end{array}\right) \text{ for } x <= 0 \text{ and } \bvec{u} = \left(\begin{array}{c} \rho_R \\ u_R \\ p_R \end{array}\right) \text{ for } x > 0$

By using Eq. (7), Sod’s initial conditions can be set as:

(8) $\bvec{u} = \left(\begin{array}{c} 1 \\ 0 \\ 1 \end{array}\right) \defeq \bvec{u}_L \text{ for } x <= 0 \text{ and } \bvec{u} = \left(\begin{array}{c} 0.125 \\ 0 \\ 0.1 \end{array}\right) \defeq \bvec{u}_R \text{ for } x > 0 \text{at } t=0$

We divide the solution of the problem in “5 zones”. From the left ( $x<0$ ) to the right ( $x>0$ ) of the diaphragm.

Region I
- There is no boundary of the tube. The status is always $\bvec{u}_L$ .
Region II
- Rarefaction wave. The status is continuous from the region 1 to the region 3.
Region III
- In the shock “pocket”, there is “no more shock” and the hyperbolic PDE (6) told us $u_{\mathrm{III}}=u_{\mathrm{IV}}$ are Riemann invariants. Together with Rankine-Hugoniot conditions, we know $p_{\mathrm{III}}=p_{\mathrm{IV}}$ and the density is not continuous.
Region IV
- Because of the expansion of the shock, there is shock discontinuity. The discontinuity status could be determined by Rankine-Hugoniot conditions [Wesselling01].
Region V
- There is no boundary of the tube, so the status is always $\bvec{u}_R$

To derive the analytic solution, we will begin from the region $\mathrm{IV}$ to get $\bvec{u}_{\mathrm{IV}}$ , then $\bvec{u}_{\mathrm{III}}$ and finally $\bvec{u}_{\mathrm{II}}$ .

Derivation of $\bvec{u}_{\mathrm{IV}}$ ¶

Rankine-Hugoniot conditions gives that the jump conditions must hold across a shock:

(9) $u_{shock}(\rho_{2} - \rho_{1}) = m_2 - m_1$

(10) $u_{shock}(m_2 - m_1) = \frac{{m_2}^2}{\rho_2} + p_2 - \frac{{m_1}^2}{\rho_1} - p_1$

(11) $u_{shock}(\rho_{2} E_2 - \rho_{1} E_1) = m_{2} H_{2} - m_{1} H_{1}$

If this is a stationary shock, $u_{shock} = 0$ . (9) tells us $m_2 = m_1$ .

Because $u_{shock} = 0$ and (9), from (10) we get:

$\frac{{m_2}^2}{\rho_2} + p_2 - \frac{{m_1}^2}{\rho_1} - p_1 = 0 \\ \text{divided by } m_1 \Rightarrow \frac{{m_2}^2}{\rho_{2}m_1} + \frac{p_2}{m_1} - \frac{{m_1}^2}{\rho_1{m_1}} - \frac{p_1}{m_1} = 0 \\ \text{please remember } m_1 = m_2 \Rightarrow \frac{m_{2}}{\rho_2} + \frac{p_2}{m_2} - \frac{m_{1}}{\rho_1} - \frac{p_1}{m_1}=0 \\ \text{use } m_1 = \rho_{1}{u_1}, m_2 = \rho_{2}{u_2} \Rightarrow u_2 + \frac{{\gamma}{p_2}}{\gamma{\rho_{2}{u_2}}} - u_1 - \frac{{\gamma}{p_1}}{\gamma{\rho_{1}{u_1}}} \\ \text{because } c_1 = \sqrt{\frac{{\gamma}{p_1}}{\rho_1}}, c_2 = \sqrt{\frac{{\gamma}{p_2}}{\rho_2}} \Rightarrow u_2 + \frac{{c_2}^2}{u_2} - u_1 - \frac{{c_1}^2}{u_1} = 0$

Thus, we get

(12) $u_1 - u_2 = \frac{{c_2}^2}{u_2} - \frac{{c_1}^2}{u_1}$

Since $u_{shock} = 0$ , $m_2 = m_1$ and (11), we get $H_1 = H_2$ . Use $H=h+\frac{{c}^2}{2}$ , namely $H_1=h_1+\frac{{c_1}^2}{2}$ and $H_2=h_2+\frac{{c_2}^2}{2}$ , and we could rewrite $H_1 = H_2$ as

$& H_1 = h_1+\frac{{u_1}^2}{2} = h_2+\frac{{u_2}^2}{2} = H_2 \\ & \text{Use } h = c_{p}T = \frac{c^2}{\gamma - 1} \\ & \text{that is } \quad h_1 = c_{p}T_1 = \frac{{c_1}^2}{\gamma - 1} \quad h_2 = c_{p}T_2 = \frac{{c_2}^2}{\gamma - 1} \\ \Rightarrow \quad & h_1 + \frac{{u_1}^2}{2} = \frac{{c_1}^2}{\gamma - 1} + \frac{{u_1}^2}{2} \\ \quad & h_2 + \frac{{u_2}^2}{2} = \frac{{c_2}^2}{\gamma - 1} + \frac{{u_2}^2}{2} \\ \Rightarrow \quad & \frac{{c_1}^2}{\gamma - 1} + \frac{{u_1}^2}{2} = \frac{{c_2}^2}{\gamma - 1} + \frac{{u_2}^2}{2}$

Assume $u_1 > \text{sonic speed} c_{*} > u_2$ . Because of continuity, there must be a point with the speed $u_{*}$ equal to the sound speed $c_{*}$ which satisfies:

$\frac{{c_1}^2}{\gamma - 1} + \frac{{u_{*}}^2}{2} = \frac{{c_1}^2}{\gamma - 1} + \frac{{c_{*}}^2}{2} = \frac{(\gamma-1)+2}{2(\gamma-1)}{c^2_{*}} = \frac{(\gamma+1)}{2(\gamma-1)}{c^2_{*}}$

And

(13) $\frac{{c_1}^2}{\gamma - 1} + \frac{{u_1}^2}{2} = \frac{{c_2}^2}{\gamma - 1} + \frac{{u_2}^2}{2} = \frac{(\gamma+1)}{2(\gamma-1)}{c^2_{*}}$

Now let’s try to get $c_{*}$ represented by $u_{1}$ and $u_{2}$ . Because of (13)

$& \frac{{c_1}^2}{\gamma - 1} + \frac{{u_1}^2}{2} = \frac{(\gamma+1)c_{*}}{2(\gamma-1)} \\ & \frac{{c_2}^2}{\gamma - 1} + \frac{{u_2}^2}{2} = \frac{(\gamma+1)c_{*}}{2(\gamma-1)} \\ \text{multipled by } \frac{(2\gamma-1)}{\gamma{u_1}} & \text{ and multipled by } \frac{(2\gamma-1)}{\gamma{u_2}} \text{ seperately} \\ \Rightarrow & \frac{2{c_1}^2}{\gamma{u_1}} + \frac{{u_1}(\gamma-1)}{\gamma} = \frac{(\gamma+1)c_{*}}{\gamma{u_1}} \\ & \frac{2{c_2}^2}{\gamma{u_2}} + \frac{{u_2}(\gamma-1)}{\gamma} = \frac{(\gamma+1)c_{*}}{\gamma{u_2}} \\ \Rightarrow & \frac{2{c_1}^2}{\gamma{u_1}} = \frac{(\gamma+1)c_{*}}{\gamma{u_1}} - \frac{{u_1}(\gamma-1)}{\gamma} \\ & \frac{2{c_2}^2}{\gamma{u_2}} = \frac{(\gamma+1)c_{*}}{\gamma{u_2}} - \frac{{u_2}(\gamma-1)}{\gamma} \\ \Rightarrow & \frac{{c_1}^2}{\gamma{u_1}} = \frac{(\gamma+1)c_{*}}{2\gamma{u_1}} - \frac{{u_1}(\gamma-1)}{2\gamma} = [\frac{(\gamma+1)c_{*}}{\gamma-1}+{u_1}^2] (\frac{(\gamma-1)}{2\gamma{u_1})}) \\ & \frac{{c_2}^2}{\gamma{u_2}} = \frac{(\gamma+1)c_{*}}{2\gamma{u_2}} - \frac{{u_2}(\gamma-1)}{2\gamma} = [\frac{(\gamma+1)c_{*}}{\gamma-1}+{u_2}^2] (\frac{(\gamma-1)}{2\gamma{u_2})}) \\ \Rightarrow & \frac{{c_1}^2}{\gamma{u_1}} - \frac{{c_2}^2}{\gamma{u_2}} = \frac{(\gamma+1)c_{*}}{2\gamma{u_1}} - \frac{{u_1}(\gamma-1)}{2\gamma} - \frac{(\gamma+1)c_{*}}{2\gamma{u_2}} + \frac{{u_2}(\gamma-1)}{2\gamma}$

please recall (12), thus

$u_2 - u_1 & = \frac{(\gamma+1)c_{*}}{2\gamma{u_1}} - \frac{{u_1}(\gamma-1)}{2\gamma} - \frac{(\gamma+1)c_{*}}{2\gamma{u_2}} + \frac{{u_2}(\gamma-1)}{2\gamma} \\ \Rightarrow & c_{*}(\frac{(\gamma+1)}{2\gamma{u_1}} - \frac{(\gamma+1)}{2\gamma{u_2}}) = u_2\frac{\gamma+1}{2\gamma} + u_1\frac{\gamma+1}{2\gamma} \\ \Rightarrow & c_{*}(\frac{1}{u_1}-\frac{1}{u_2}) = u_2 - u_1 \\ \Rightarrow & c_{*} = {u_1}{u_2}$

The relation

(14) $c_{*} = {u_1}{u_2}$

is called Prandel-Meyer relation. It means the flow at one side of a shock must be supersonic, and the other side must be subsonic.

Now defining the Mach number $M_1 = \frac{u_1}{c_1}$ and $M^{*} = \frac{u}{c_*}$ , we get from (13):

$& \frac{{c_1}^2}{(\gamma - 1)} + \frac{{u_1}^2}{2} = \frac{(\gamma+1)}{2(\gamma-1)}{c^2_{*}} \\ \Rightarrow & \frac{{c_1}^2}{(u_1)^2}\frac{(u_1)^2}{\gamma - 1} + \frac{{u_1}^2}{2} = \frac{(\gamma+1)}{2(\gamma-1)}{c^2_{*}} \\ \text{dividied by } c^2_{*} \Rightarrow & \frac{{c_1}^2}{(u_1)^2}\frac{(u_1)^2}{c^2_{*}(\gamma - 1)} + \frac{{u_1}^2}{(c^2_{*})2} = \frac{(\gamma+1)}{2(\gamma-1)} \\ & \frac{M^*_1}{M_1(\gamma-1)} + \frac{M^{*2}_{1}}{2} = \frac{{u_1}^2}{(c^2_{*})2}$

This parcel solvcon.parcel.gas contains code that performs computational fluid dynamics (CFD) for gas flows by using the CESE method. The model equations currently are the Euler equations. This parcel will be updated to use the Navier-Stokes equations in the future. See The Euler Equations (Under Development) for detail.

This parcel has the following examples:

Simulation Settings¶

solvcon.parcel.gas.register_arrangement()¶: This is an alias to the instance method GasCase.register_arrangement(), which inherits from the class solvcon.MeshCase. See solvcon.MeshCase.register_arrangement() for implementation detail.

class solvcon.parcel.gas.GasCase(**kw)¶: See case.GasCase for implementation detail.

class solvcon.parcel.gas.GasSolver(blk, **kw)¶: See solver.GasSolver for implementation detail.

Boundary-Condition Treatments¶

class solvcon.parcel.gas.GasNonrefl¶: See boundcond.GasNonrefl for implementaion detail.

class solvcon.parcel.gas.GasWall¶: See boundcond.GasWall for implementation detail.

class solvcon.parcel.gas.GasInlet¶: See boundcond.GasInlet for implementation detail.

Callbacks¶

class solvcon.parcel.gas.ProbeHook¶: See probe.ProbeHook for implementation detail.

class solvcon.parcel.gas.DensityInitAnchor¶: See physics.DensityInitAnchor for implementation detail.

class solvcon.parcel.gas.PhysicsAnchor¶: See physics.PhysicsAnchor for implementation detail.

class solvcon.parcel.gas.MeshInfoHook¶: See inout.MeshInfoHook for implementation detail.

class solvcon.parcel.gas.ProgressHook¶: See inout.ProgressHook for implementation detail.

class solvcon.parcel.gas.FillAnchor¶: See inout.FillAnchor for implementation detail.

class solvcon.parcel.gas.CflHook¶: See inout.CflHook for implementation detail.

class solvcon.parcel.gas.PMarchSave¶: See inout.PMarchSave for implementation detail.

Internal Refenrence¶

Simulation¶

The modules defining the physics process and the numerical methods are:

`case`¶

class solvcon.parcel.gas.case.GasCase(**kw)¶

Temporal loop for the gas-dynamic solver.

Inheritance diagram of GasCase

`solver`¶

class solvcon.parcel.gas.solver.GasSolver(blk, **kw)¶

Spatial loops for the gas-dynamics solver.

Inheritance diagram of GasSolver

__init__(blk, **kw)¶

Create a _algorithm.GasAlgorithm object.

>>> # create a valid solver as the test fixture.
>>> from solvcon.testing import create_trivial_2d_blk
>>> blk = create_trivial_2d_blk()
>>> blk.clgrp.fill(0)
>>> blk.grpnames.append('blank')
>>> class SubSolver(GasSolver):
...     pass
>>> svr = SubSolver(blk)
>>> # number of equations.
>>> svr.neq
4
>>> # valid GasAlgorithm.
>>> svr.alg is not None
True

Boundary-Condition Treatment¶

`boundcond`¶

All these treatments are in solvcon.parcel.gas.boundcond module.

class solvcon.parcel.gas.boundcond.GasBC(**kw)¶

Base class for all boundary conditions of the gas solver.

Inheritance diagram of GasBC

The base class for all boundary-condition treatments in this module. It should not be used in simulations directly.

class solvcon.parcel.gas.boundcond.GasNonrefl(**kw)¶

class solvcon.parcel.gas.boundcond.GasWall(**kw)¶

class solvcon.parcel.gas.boundcond.GasInlet(**kw)¶

Callback¶

Useful callbacks are included in:

`probe`¶

class solvcon.parcel.gas.probe.ProbeHook(cse, **kw)¶

Point probe.

Inheritance diagram of ProbeHook

`physics`¶

class solvcon.parcel.gas.physics.DensityInitAnchor(svr, **kw)¶

Initialize only density.

FIXME: Give me doctests.

Inheritance diagram of DensityInitAnchor

class solvcon.parcel.gas.physics.PhysicsAnchor(svr, **kw)¶

Calculates physical quantities for output. Implements (i) provide() and (ii) postfull() methods.

FIXME: I should be more integrated with GasSolver.

Variables:	gasconst – gas constant.

Inheritance diagram of PhysicsAnchor

`inout`¶

class solvcon.parcel.gas.inout.MeshInfoHook(cse, show_bclist=False, perffn=None, **kw)¶

Print mesh information.

Inheritance diagram of MeshInfoHook

class solvcon.parcel.gas.inout.ProgressHook(cse, linewidth=50, **kw)¶

Print simulation progess.

Inheritance diagram of ProgressHook

class solvcon.parcel.gas.inout.FillAnchor(svr, mappers=None, **kw)¶

Fill the specified arrays of a GasSolver with corresponding value.

Inheritance diagram of FillAnchor

class solvcon.parcel.gas.inout.CflAnchor(svr, rsteps=None, **kw)¶

Counting CFL numbers. Use MeshSolver.marchret to return results. Overrides postmarch() method.

Pair with CflHook.

Inheritance diagram of CflAnchor

__init__(svr, rsteps=None, **kw)¶

>>> from solvcon.testing import create_trivial_2d_blk
>>> from solvcon.solver import MeshSolver
>>> svr = MeshSolver(create_trivial_2d_blk())
>>> ank = CflAnchor(svr)
Traceback (most recent call last):
    ...
TypeError: int() argument must be a string or a number, not 'NoneType'
>>> ank = CflAnchor(svr, 1)
>>> ank.rsteps
1

class solvcon.parcel.gas.inout.CflHook(cse, name='cfl', cflmin=0.0, cflmax=1.0, fullstop=True, rsteps=None, **kw)¶

Makes sure CFL number is bounded and print averaged CFL number over time. Reports CFL information per time step and on finishing. Overrides (i) postmarch() and (ii) postloop() methods.

Pair with CflAnchor.

Inheritance diagram of CflHook

class solvcon.parcel.gas.inout.MarchSaveAnchor(svr, anames=None, compressor=None, fpdtype=None, psteps=None, vtkfn_tmpl=None, **kw)¶

Save solution data into VTK XML format for a solver.

Inheritance diagram of MarchSaveAnchor

class solvcon.parcel.gas.inout.PMarchSave(cse, anames=None, compressor='gz', fpdtype=None, altdir='', altsym='', vtkfn_tmpl=None, **kw)¶

Save the geometry and variables in a case when time marching in parallel VTK XML format.

Inheritance diagram of PMarchSave

Development¶

Issue tracker: https://bitbucket.org/solvcon/solvcon/issues
Users’ mailing list: solvcon@googlegroups.com (or http://groups.google.com/group/solvcon)
Python Style Guide
Jenkins Slave Setup
Hidden Applications contain documents of the parcels that are not actively maintained at the time being.

References¶

Bibliography
History
Papers and presentations:
- Published Applications of SOLVCON
- Yung-Yu Chen, A Multi-Physics Software Framework on Hybrid Parallel Computing for High-Fidelity Solutions of Conservation Laws, Ph.D. Thesis, The Ohio State University, United States, Aug. 2011. (OhioLINK)
- PyCon US 2011 talk: slides and video
- Yung-Yu Chen, David Bilyeu, Lixiang Yang, and Sheng-Tao John Yu, “SOLVCON: A Python-Based CFD Software Framework for Hybrid Parallelization”, 49th AIAA Aerospace Sciences Meeting, January 4-7 2011, Orlando, Florida. AIAA Paper 2011-1065
The CESE method:
- The CE/SE working group: http://www.grc.nasa.gov/WWW/microbus/
- The CESE research group at OSU: http://cfd.solvcon.net/research.html
- Selected papers:
  - Sin-Chung Chang, “The Method of Space-Time Conservation Element and Solution Element – A New Approach for Solving the Navier-Stokes and Euler Equations”, Journal of Computational Physics, Volume 119, Issue 2, July 1995, Pages 295-324. doi: 10.1006/jcph.1995.1137
  - Xiao-Yen Wang, Sin-Chung Chang, “A 2D Non-Splitting Unstructured Triangular Mesh Euler Solver Based on the Space-Time Conservation Element and Solution Element Method”, Computational Fluid Dynamics Journal, Volume 8, Issue 2, 1999, Pages 309-325.
  - Zeng-Chan Zhang, S. T. John Yu, Sin-Chung Chang, “A Space-Time Conservation Element and Solution Element Method for Solving the Two- and Three-Dimensional Unsteady Euler Equations Using Quadrilateral and Hexahedral Meshes”, Journal of Computational Physics, Volume 175, Issue 1, Jan. 2002, Pages 168-199. doi: 10.1006/jcph.2001.6934
Related Links
Other Software for Solving PDEs

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
Neither the name of the SOLVCON nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Contributors

Yung-Yu Chen <yyc@solvcon.net>
Sheng-Tao John Yu <yu.274@osu.edu>
David Bilyeu <bilyeu.4@osu.edu>
Tzu-Wei Lin <evolutionlin@gmail.com>
Po-Hsien Lin <lin.880@buckeyemail.osu.edu>
Taihsiang Ho <tai271828@gmail.com>

Navigation

SOLVCON¶

Install¶

Documents¶

Solver Architecture¶

Unstructured Meshes¶

Entities¶

Data Structure Defined in solvcon.block¶

Low-Level Interface to C Defined in solvcon.mesh¶

Numerical Code¶

Initialize¶

Time-March¶

Arrangement¶

Time-Marching¶

Parallel Computing¶

Low-Level Interface to C for BC¶

References¶

ustmesh_2d_sample.geo¶

Input and Output Facilities¶

Gas Dynamics Parcel (Under Development)¶

The Euler Equations (Under Development)¶

Reflection of Oblique Shock Wave¶

Driving Script¶

Arrangement¶

Analytical Solution¶

Case Instantiation¶

Callback Configuration¶

View Results¶

Oblique Shock Relation¶

Full Listing of the Driving Script¶

Sod’s Shock Tube (Under Development)¶

Derivation of ¶

Simulation Settings¶

Boundary-Condition Treatments¶

Callbacks¶

Internal Refenrence¶

Simulation¶

Boundary-Condition Treatment¶

Callback¶

Development¶

References¶

Quick search

Navigation

Data Structure Defined in `solvcon.block`¶

Low-Level Interface to C Defined in `solvcon.mesh`¶

Low-Level Interface to C for `BC`¶

Derivation of $\bvec{u}_{\mathrm{IV}}$ ¶