Simulation framework¶

This page documents the Ulula hydro solver and the implemented algorithms.

Overview¶

The Ulula simulation framework is basically implemented within a single class, called Simulation. The HydroScheme class contains information about the chosen algorithm. The user does not need to interact directly with the Simulation class because the simulation workflow is already implemented in the ulula.run.run() runtime function (see Running Ulula). If this workflow is not general enough, however, the user can run a simulation by implementing and/or modifying the following steps:

Create a HydroScheme object that contains all user-defined algorithmic settings
Create the Simulation (passing the HydroScheme to the constructor)
Set up the domain and variable arrays with the setDomain() function, which sets the size of the domain in pixels and physical coordinates, as well as the boundary conditions
Set the fluid properties (such as the adiabatic index) using the setFluidProperties() function
Set the initial conditions in primitive variables (into the simulation object’s V array); this operation is implemented in each problem setup class (see Hydro problem setups). The primitive variable fields are listed below.
Execute the timestep() function until the simulation is finished; create plots as desired

The options for hydro algorithms are described in the HydroScheme class below. The documentation of Simulation contains further details about the inner workings of the Ulula hydro solver.

Ulula can save the current state of a simulation to an hdf5 file (see the save() function). The load() function loads such a file and recovers a Simulation object identical to the one saved. The simulation can then be plotted or continued, meaning that the files serve as both snapshot and restart files.

Fluid quantities¶

The following identifiers are used for fluid variables throughout the code:

Abbreviation	Quantity
Primitive variables
`DN`	Density
`VX`	X-velocity
`VY`	Y-velocity
`PR`	Pressure
Conserved variables
`MS`	Density (here called mass for consistency)
`MX`	X-momentum
`MY`	Y-momentum
`ET`	Total energy (thermal + kinetic + potentials)

These abbreviations are used to indicate quantities when plotting, although a few additional plottable quantities are defined (see Plotting). Similarly, the abbreviations are used when setting fluid variables in a Simulation object (see the q_prim and q_cons fields).

Hydro schemes¶

The following tables list the algorithmic choices implemented in Ulula. Broadly speaking, Ulula uses a Godunov solver with either piecewise-constant or piecewise-linear state reconstruction and a either Euler or 2nd-order time integration (MUSCL-Hancock).

Spatial reconstruction schemes¶
Identifier	Description
`const`	Piecewise-constant; the scheme is 1st-order in space because the state within each cell is not interpolated. This hydro scheme tends to be highly diffusive and is impelemented for test purposes.
`linear`	Piecewise-linear; linear interpolation within each cell, the scheme is 2nd-order in space. The slope limiter decides how aggressive the reconstruction is in terms of accuracy vs. stability.

Slope limiters¶
Identifier	Description
`none`	No limiter; results in a highly unstable scheme because the interpolation can lead to negative values and other artifacts. Implemented for demonstration purposes only
`minmod`	Minimum modulus; the most conservative limiter, takes the shallower of the left and right slopes
`vanleer`	The van Leer limiter; an intermediate between the minmod and mc limiters
`mc`	Monotonized central; the most aggressive limiter, takes the central slope unless the slopes are very different

Riemann solvers¶
Identifier	Description
`hll`	The Harten-Lax-van Leer (HLL) Riemann solver; this simple algorithm considers only the fastest waves to the left and right and constructs an intermediate state, ignoring contact discontinuities.

Time integration schemes¶
Identifier	Description
`euler`	First-order time integration; the fluid state is advanced by a full timestep without any attempt to time-average the Godunov fluxes.
`hancock`	Combined with linear reconstruction, this choice gives the MUSCL-Hancock scheme, where the states at the left and right cell edges are advanced by a half timestep. The change in time is computed to linear order from the primitive-variable form of the Euler equations, applied to the left and right states separately.
`hancock_cons`	The same as `hancock`, but using the conserved formulation where the fluxes corresponding to the left and right cell-edge states are computed and differenced to get an evolution in the fluid state. This formulation should be identical to the primitive formulation up to numerical errors, which provides a check of the algorithms. The conservative version is slower due to a number of primitive-conserved conversions; thus, the primitive version is preferred.

The purpose of Ulula is to experiment with hydro algorithms, including their failure. Thus, the code allows sub-optimal (aka crazy) combinations of the parameters listed above. For example, it allows setting no limiter or combining piecewise-constant states with a Hancock-style time integration, which results in an unstable scheme. The user’s algorithmic choices are contained in the HydroScheme class:

class ulula.simulation.HydroScheme(reconstruction='const', limiter='minmod', riemann='hll', time_integration='euler', cfl=0.8)¶

Container class for hydro algorithms

Parameters

reconstruction: string: Reconstruction algorithm; see listing for valid choices
limiter: string: Slope limiter algorithm; see listing for valid choices
riemann: string: Riemann solver; see listing for valid choices
time_integration: string: Time integration scheme; see listing for valid choices
cfl: float: CFL number (must be between 0 and 1); determines the timestep as CFL number times cell size divided by the maximum signal speed in the domain

Simulation class¶

class ulula.simulation.Simulation(hydro_scheme)¶

Main class for the Ulula hydro solver

This class contains all simulation data and routines. The internal fields have the following meaning (after the hydro scheme, domain, fluid properties, and initial conditions have been set):

Field	Meaning
Scalars
`dx`	Width of cells (same in both x and y directions)
`nq`	Number of fluid variables (typically 4 for rho, vx, vy, P)
`nx`	Number of cells in the x-direction
`ny`	Number of cells in the y-direction
`ng`	Number of ghost cells around each edge
`xlo`	First index of physical grid in x-direction (without left ghost zone)
`xhi`	Last index of physical grid in x-direction (without right ghost zone)
`ylo`	First index of physical grid in y-direction (without bottom ghost zone)
`yhi`	Last index of physical grid in y-direction (without top ghost zone)
`t`	Current time of the simulation
`step`	Current step counter
`last_dir`	Direction of last sweep in previous timestep (x=0, y=1)
`gamma`	Adiabatic index
`gm1`	gamma - 1
`gm1_inv`	1 / (gamma - 1)
Settings
`bc_type`	Type of boundary condition (‘periodic’ or ‘outflow’)
`hs`	HydroScheme object
1D vectors
`x`	Array of x values at cell centers (dimensions [nx + 2 ng])
`y`	Array of y values at cell centers (dimensions [ny + 2 ng])
3D vectors
`U`	Vector of conserved fluid variables (dimensions [nq, nx + 2 ng, ny + 2 ng])
`V`	Vector of primitive fluid variables (dimensions [nq, nx + 2 ng, ny + 2 ng])
`V_im12`	Cell-edge states at left side (same dimensions as V)
`V_ip12`	Cell-edge states at right side (same dimensions as V)
Slices
`slc1dL`	1D slices for idir [0, 1], physical domain, shifted one cell left
`slc1dR`	1D slices for idir [0, 1], physical domain, shifted one cell right
`slc1dC`	1D slices for idir [0, 1], physical domain
`slc3dL`	3D slices for idir [0, 1], physical domain, shifted one cell left
`slc3dR`	3D slices for idir [0, 1], physical domain, shifted one cell right
`slc3dC`	3D slices for idir [0, 1], physical domain
`slc3aL`	3D slices for idir [0, 1], total domain, shifted one cell left
`slc3aR`	3D slices for idir [0, 1], total domain, shifted one cell right
`slc3aC`	3D slices for idir [0, 1], total domain
`slc3fL`	3D slice of flux vector from left interface
`slc3fR`	3D slice of flux vector from right interface
Dictionaries for fluid variables
`q_prim`	Dictionary containing the indices of the primitive fluid variables in the `V` array
`q_cons`	Dictionary containing the indices of the conserved fluid variables in the `U` array

The constructor takes the following parameters:

Parameters

hydro_scheme: HydroScheme: Container class for algorithmic choices

Methods

`cflCondition`(self)	Compute the size of the next timestep
`conservedToPrimitive`(self, U, V)	Convert conserved to primitive variables
`emptyArray`(self[, nq])	Get an empty array for fluid variables
`enforceBoundaryConditions`(self)	Enforce boundary conditions after changes
`fluxVector`(self, idir, V)	Convert the flux vector F(V)
`limiterMC`(self, sL, sR, slim)	Monotonized-central limiter
`limiterMinMod`(self, sL, sR, slim)	Minimum-modulus limiter
`limiterNone`(self, sL, sR, slim)	Non-limiter (central derivative)
`limiterVanLeer`(self, sL, sR, slim)	The limiter of van Leer
`maxSpeedInDomain`(self)	Largest signal speed in domain
`primitiveEvolution`(self, idir, V, dV_dx)	Linear approximation of the Euler equations
`primitiveToConserved`(self, V, U)	Convert primitive to conserved variables
`primitiveToConservedRet`(self, V)	Convert primitive to new conserved array
`reconstructionConst`(self, idir, dt)	Piecewise-constant reconstruction
`reconstructionLinear`(self, idir, dt)	Piecewise-linear reconstruction
`riemannSolverHLL`(self, idir, VL, VR)	The HLL Riemann solver
`save`(self[, filename])	Save the current state of a simulation
`setDomain`(self, nx, ny[, xmin, xmax, ymin, …])	Set the physical and numerical size of the domain
`setFluidProperties`(self[, gamma])	Set the physical properties of the fluid
`soundSpeed`(self, V)	Sound speed
`timestep`(self[, dt])	Advance the fluid state by one timestep
`xyGrid`(self)	Get a grid of the x and y cell center positions

setDomain(self, nx, ny, xmin=0.0, xmax=1.0, ymin=0.0, bc_type='periodic')¶

Set the physical and numerical size of the domain

This function creates the memory structure for the simulation as well as pre-computed slices that index the arrays.

Parameters

nx: int: Number of grid points in x-direction
ny: int: Number of grid points in y-direction
xmin: float: Left edge in physical coordinates (code units)
xmax: float: Right edge in physical coordinates (code units)
ymin: float: Bottom edge in physical coordinates (code units)
ymax: float: Top edge in physical coordinates (code units)
bc_type: string: Type of boundary conditions; can be periodic or outflow

setFluidProperties(self, gamma=1.6666666666666667)¶

Set the physical properties of the fluid

Parameters

gamma: float: Adiabatic index of the ideal gas to be simulated; should be 5/3 for atomic gases or 7/5 for diatomic molecular gases.

emptyArray(self, nq=None)¶

Get an empty array for fluid variables

Parameters

nq: int: The number of quantities for which the array should contain space. If None, the number of fluid quantities is used (4 in two dimensions).

Returns

ret: array_like: Float array of size nq times the size of the domain including ghost cells. If nq == 1, the first dimension is omitted.

xyGrid(self)¶

Get a grid of the x and y cell center positions

This function returns two arrays with the x and y positions at each grid point. These arrays can be convenient when setting the initial conditions.

Returns

x: array_like: 2D array with x positions of all cells (including ghost cells)
y: array_like: 2D array with x positions of all cells (including ghost cells)

enforceBoundaryConditions(self)¶

Enforce boundary conditions after changes

This function fills the ghost cells with values from the physical domain. For periodic BCs, those originate from the other side; for outflow BCs, they are copied from the edge of the physical domain. This function must be executed at each timestep.

primitiveToConserved(self, V, U)¶

Convert primitive to conserved variables

This function takes the input and output arrays as parameters instead of assuming that it should use the main V and U arrays. In some cases, conversions need to be performed on other fluid states.

Parameters

V: array_like: Input array of primitive fluid variables with first dimension nq (rho, vx, vy, P…)
U: array_like: Output array of fluid variables with first dimension nq (rho, u * vx…)

primitiveToConservedRet(self, V)¶

Convert primitive to new conserved array

Same as primitiveToConserved(), but creating the conserved output array.

Parameters

V: array_like: Input array of primitive fluid variables with first dimension nq (rho, vx, vy, P…)

Returns

U: array_like: Array of fluid variables with first dimension nq (rho, u * vx…) and same dimensions as input array.

conservedToPrimitive(self, U, V)¶

Convert conserved to primitive variables

This function takes the input and output arrays as parameters instead of assuming that it should use the main U and V arrays. In some cases, conversions need to be performed on other fluid states.

Parameters

U: array_like: Input array of conserved fluid variables with first dimension nq (rho, u * vx…)
V: array_like: Output array of primitive fluid variables with first dimension nq (rho, vx, vy, P…)

fluxVector(self, idir, V)¶

Convert the flux vector F(V)

The flux of the conserved quantities density, momentum, and total energy as a function of a primitive fluid state.

Parameters

idir: int: Direction of sweep (0 = x, 1 = y)
V: array_like: Input array of primitive fluid variables with first dimension nq (rho, vx, vy, P…)

Returns

F: array_like: Array of fluxes with first dimension nq and same dimensions as input array.

primitiveEvolution(self, idir, V, dV_dx)¶

Linear approximation of the Euler equations

Instead of the conservation-law form, we can also think of the Euler equations as \(dV/dt + A(V) dV/dx = S\). This function returns \(\Delta V/ \Delta t\) given an input state and a vector of spatial derivatives \(\Delta V/ \Delta x\). The result is used in the Hancock step.

Parameters

idir: int: Direction of sweep (0 = x, 1 = y)
V: array_like: Array of primitive fluid variables with first dimension nq (rho, vx, vy, P…)
dV_dx: array_like: Array of derivative of fluid variables with first dimension nq

Returns

dV_dt: array_like: Array of linear approximation to time evolution of fluid variables, with same dimensions as input arrays.

soundSpeed(self, V)¶

Sound speed

Parameters

V: array_like: Input array of primitive fluid variables with first dimension nq (rho, vx, vy, P…)

Returns

cs: array_like: Array of sound speed with first dimension nq and same dimensions as input array.

maxSpeedInDomain(self)¶

Largest signal speed in domain

This function returns the largest possible signal speed anywhere in the domain. It evaluates the sound speed and adds it to the absolute x and y velocities. We do not need to add those velocities in quadrature since we are taking separate sweeps in the x and y directions. Thus, the largest allowed timestep is determined by the largest speed in either direction.

Parameters

V: array_like: Input array of primitive fluid variables with first dimension nq (rho, vx, vy, P…)

Returns

c_max: float: Largest possible signal speed in the domain.

reconstructionConst(self, idir, dt)¶

Piecewise-constant reconstruction

Piecewise-constant means no reconstruction. The left/right cell edge value arrays are already set to the cell-centered values so that this function does nothing at all. It serves as a placeholder to which the reconstruction function pointer can be set.

Parameters

idir: int: Direction of sweep (0 = x, 1 = y)
dt: float: Timestep

reconstructionLinear(self, idir, dt)¶

Piecewise-linear reconstruction

This function creates left and right cell-edge states based on the cell-centered states. It first computes the left and right slopes, uses a slope limiter to determine the limited slope to use, and interpolates linearly within each cell.

If the time integration scheme is Hancock, the reconstructed edge states are also advanced by half a timestep to get 2nd-order convergence in the flux calculation. There are two ways to perform the Hancock step. The more conventionally described way is to take the fluxes according to the L/R states as an approximation for the flux differential across the cell (the hancock_cons integration scheme). The differential is then used to updated the conserved cell-edge states. However, this method necessitates a primitive->conserved->primitive conversion and a flux calculation. By contrast, the so-called primitive Hancock method uses the Euler equations in primitive variables to estimate the change in time from the change across the cell (see primitiveEvolution()). The two methods should give almost identical results, but the primitive version is noticeably faster.

Parameters

idir: int: Direction of sweep (0 = x, 1 = y)
dt: float: Timestep

limiterNone(self, sL, sR, slim)¶

Non-limiter (central derivative)

This limiter is the absence thereof: it does not limit the left and right slopes but returns their average (the central derivative). This generally produces unstable schemes but is implemented for testing and demonstration purposes.

Parameters

sL: array_like: Array of left slopes
sR: array_like: Array of right slopes
slim: array_like: Output array of limited slope; must have same dimensions as sL and sR.

limiterMinMod(self, sL, sR, slim)¶

Minimum-modulus limiter

The most conservative limiter, which always chooses the shallower out of the left and right slopes.

Parameters

sL: array_like: Array of left slopes
sR: array_like: Array of right slopes
slim: array_like: Output array of limited slope; must have same dimensions as sL and sR.

limiterVanLeer(self, sL, sR, slim)¶

The limiter of van Leer

An intermediate limiter that is less conservative than minimum modulus but more conservative than monotonized central.

Parameters

sL: array_like: Array of left slopes
sR: array_like: Array of right slopes
slim: array_like: Output array of limited slope; must have same dimensions as sL and sR.

limiterMC(self, sL, sR, slim)¶

Monotonized-central limiter

As the name suggests, this limiter chooses the central derivative wherever possible, but reduces its slope where it would cause negative cell-edge values. This limiter leads to the sharpest solutions but is also the least stable.

Parameters

sL: array_like: Array of left slopes
sR: array_like: Array of right slopes
slim: array_like: Output array of limited slope; must have same dimensions as sL and sR.

riemannSolverHLL(self, idir, VL, VR)¶

The HLL Riemann solver

The Riemann solver computes the fluxes across cell interfaces given two discontinuous states on the left and right sides of each interface. The Harten-Lax-van Leer (HLL) Riemann solver is one of the simplest such algorithms. It takes into account the fastest waves traveling left and right, but it computes only one intermediate state that ignores contact discontinuities.

Parameters

idir: int: Direction of sweep (0 = x, 1 = y)
VL: array_like: Array of primitive state vectors on the left sides of the interfaces
VR: array_like: Array of primitive state vectors on the right sides of the interfaces

Returns

flux: array_like: Array of conservative fluxes across interfaces; has the same dimensions as VL and VR.

cflCondition(self)¶

Compute the size of the next timestep

This function computes the maximum signal speed anywhere in the domain and sets a timestep based on the CFL condition. This routine

Returns

dt: float: Size of the next timestep

timestep(self, dt=None)¶

Advance the fluid state by one timestep

This timestepping routine implements a dimensionally split scheme, where we alternate the order of execution with each timestep to maintain a 2nd-order scheme in time (xy-yx-xy and so on). In each direction, we reconstruct the cell-edge states, compute the conservative Godunov fluxes with the Riemann solver, and add the flux difference to the converved fluid variables. Technically, we should make sure that the same timestep is taken for both the xy and yx steps because otherwise the first-order error terms do not perfectly cancel. However, we omit this complication, given that the timestep should not change drastically between steps.

Parameters

dt: float: Size of the timestep to be taken; if None the timestep is computed from the CFL condition using the cflCondition() function. This timestep should be used in most circumstances, but sometimes we wish to take a manually set timestep, for example, to output a file or plot. Thus, the two functions are separated.

Returns

dt: float: The timestep taken

save(self, filename=None)¶

Save the current state of a simulation

Parameters

filename: str: Output filename; auto-generated if None

Input/Output¶

Files are saved using the save() function in a Simulation object. The file format is relatively self-explanatory. The attributes of the hydro_scheme group correspond to the parameters of the HydroScheme object. The attributes of the domain, physics, and run groups have the same meaning as in the Simulation object. The code group contains the version of Ulula that the filetype corresponds to.

The grid group contains the 2D grid data for the given simulation in code units. The field names correspond to the abbreviations listed in the “Fluid quantities” section above. Only the physical domain is included (without ghost cells).

A file can be loaded into a new simulation object with the following function. It checks the file version against the current version of Ulula; if the file version is too old for the file to be compatible, the loading process is aborted.

ulula.simulation.load(filename)¶

Load a snapshot file into a simulation object

Parameters

filename: str: Input filename

Returns

sim: Simulation: Object of type Simulation