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 Simulation object
Optionally update the hydro scheme, equation of state, code units, and gravity using the functions setHydroScheme(), setEquationOfState(), setCodeUnits(), and setGravityMode().
Execute the setDomain() function, which sets the domain in pixels and physical coordinates and creates the variable arrays.
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.
If the physical setup contains gravity, execute setGravityPotentials().
Execute the timestep() function until the simulation is finished; create plots as desired.

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.

Hydrodynamics background

The fundamental purpose of any hydrodynamics code is to solve the Euler equations. We can understand these equations as a series of conservation laws, which we can most generally express as

\[\frac{\partial}{\partial t} (\mathrm{density\ of\ }Q) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\mathrm{flux\ of\ }Q) = 0 \,.\]

In words, if the flux of a conserved quantity \(Q\) (e.g., mass) diverges, the density of that quantity decreases. If the flux converges (negative divergence), the density increases. The Euler equations can be expressed as the following conservation laws for mass, momentum, and energy:

\[\begin{split}\begin{align} \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho {\pmb u}) &= 0 \\ \frac{\partial (\rho {\pmb u})}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho {\pmb u} \otimes {\pmb u} + {\pmb I} P) &= -\rho \boldsymbol{\nabla} \Phi \\ \frac{\partial E}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} ([E + P] {\pmb u}) &= \rho \frac{\partial \Phi}{\partial t} \\ \end{align}\end{split}\]

where \(\rho\) is density, \({\pmb u}\) the velocity vector, \(P\) pressure, and \(E\) is the total energy, defined as

\[E = \rho({\pmb u}^2 / 2 + \epsilon + \Phi) \,.\]

Here \(\epsilon\) is the internal energy per unit mass. How this translates into temperature and pressure is governed by the microscopic particle properties of the gas, which are summarized in an equation of state. In Ulula, we always assume an ideal gas equation of state,

\[P = \rho \epsilon (\gamma - 1)\]

where \(\gamma\) is typically \(5/3\) for an atomic gas but can be changed by the user. The non-zero right-hand-side of the conservation laws above means that some quantities are not always conserved. In particular, in the presence of a gravitational potential \(\Phi\), momentum and energy are added to the gas due to gravity. In terms of how the equations are solved in Ulula (and most other finite-volume hydro codes), it is easier to return to the abstract form of a conservation law above and write the three conservation laws as a single vector equation,

\[\frac{\partial {\pmb U}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\mathcal{F}({\pmb U})) = S({\pmb U}) \,,\]

where \({\pmb U}\) is the vector of conserved quantities (in 2D, in our case), \(\mathcal{F}\) is the flux vector, and \(S\) the source vector,

\[\begin{split}{\pmb U} \equiv \left(\begin{array}{c} \rho \\ \rho u_{\rm x} \\ \rho u_{\rm y} \\ E \end{array}\right) \qquad \mathcal{F}({\pmb U}) \equiv \left(\begin{array}{c} \rho \pmb{u} \\ \rho u_{\rm x} {\pmb u} + \delta_{xj} P \\ \rho u_{\rm y} {\pmb u} + \delta_{yj} P \\ (E + P) {\pmb u} \end{array}\right) \qquad {\pmb S} \equiv \left(\begin{array}{c} 0 \\ - \rho\ \partial \Phi / \partial x \\ - \rho\ \partial \Phi / \partial y \\ \rho\ \partial \Phi / \partial t \end{array}\right)\end{split}\]

The Kronecker \(\delta\) means that \(P\) is only added to the x-component of the flux vector for \(u_{\rm x}\) and so on. The basic MO of a code such as Ulula is to update the conserved quantities \({\pmb U}\) by estimating the time-integrated fluxes \(\mathcal{F}({\pmb U})\) and taking their differences across cells. As shown by the equations above, in the process of computing those updates from the flux and source vectors we do need to know pressure, velocity and so on as well as the conserved variables. We thus frequently convert to and from so-called primitive variables,

\[\begin{split}{\pmb V} \equiv \left(\begin{array}{c} \rho \\ u_{\rm x} \\ u_{\rm y} \\ P \end{array}\right)\end{split}\]

by using the definition of total energy and the equation of state (see above). All told, we need to keep track of quite a few quantities in Ulula. The following table shows these fluid variables and abbreviations used throughout the code, as well as their units in terms of length, time, and mass:

Symbol	Abbreviation	Quantity	Units
Primitive variables
\(\rho\)	`DN`	Density	\(m/l^3\)
\(u_{\rm x}\)	`VX`	X-velocity	\(l/t\)
\(u_{\rm y}\)	`VY`	Y-velocity	\(l/t\)
\(P\)	`PR`	Pressure	\(m/l/t^2\)
Conserved variables
\(\rho\)	`MS`	Mass density (same as density)	\(m/l^3\)
\(\rho u_{\rm x}\)	`MX`	X-momentum density	\(m/l^2/t\)
\(\rho u_{\rm y}\)	`MY`	Y-momentum density	\(m/l^2/t\)
\(E\)	`ET`	Total energy density	\(m/l/t^2\)
Gravity
\(\Phi\)	`GP`	Potential per unit mass	\(l^2/t^2\)
\(\partial \Phi / \partial x\)	`GX`	Potential gradient in the X-direction	\(l/t^2\)
\(\partial \Phi / \partial y\)	`GY`	Potential gradient in the Y-direction	\(l/t^2\)

The abbreviations are used to select variables in the Plotting module and when setting the initial conditions.

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='linear', limiter='mc', riemann='hll', time_integration='hancock', cfl=0.8, cfl_max=0.95, cfl_reduce_factor=1.2, cfl_max_attempts=3)

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), which determines the maximum allowable timestep such that the distance traveled by the maximum signal speed in the domain does not exceed the CFL number times the size of a cell.
cfl_max: float: Maximum CFL number (must be between 0 and 1, and greater than cfl). While each timestep is (typically) set to satisfy the CFL condition, each timestep actually consists of two sweeps in the two dimensions (x and y). Since the fluid state changes during the first sweep, satisfying the CFL condition at first does not guarantee that it is still satisfied during the second sweep. To avoid repeating the first sweep, we tolerate an actual, updated CFL factor that is larger than cfl, but it must still be smaller than cfl_max and smaller than unity, because exceeding unity will definitely break the hydro solver. Note, however, that setting cfl_max to a value close to unity (e.g., 0.99) may still lead to instabilities. On the other hand, choosing cfl and cfl_max to be close will mean that more timesteps need to be repeated, which slows down the code.
cfl_reduce_factor: float: If cfl_max is exceeded during the second sweep, we reduce the previous estimate of the timestep by this factor. Must be larger than unity.
cfl_max_attempts: int: If we still encounter a CFL violation after reducing the timestep, we keep doing so cfl_max_attempts times. After the last attempt, the simulation is aborted.

Simulation class

class ulula.simulation.Simulation

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
Domain and fluid variables
`is_2d`	Whether we are running a 1D or 2D simulation
`dx`	Width of cells (same in both x and y directions)
`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)
`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
`track_pressure`	Whether we need to explicitly track pressure
`nq_hydro`	Number of hydro fluid variables (typically 4 for rho, vx, vy, P)
`nq_all`	Total number of fluid variables (including gravitational pot.)
`bc_type`	Type of boundary condition (‘periodic’, ‘outflow’, ‘wall’)
`domain_set`	Once the domain is set, numerous settings cannot be changed any more
Hydro scheme and timestepping
`hs`	HydroScheme object
`use_source_term`	Whether any source terms are active (gravity etc.)
`t`	Current time of the simulation
`step`	Current step counter
`last_dir`	Direction of last sweep in previous timestep (x=0, y=1)
Equation of state
`eos_mode`	Type of equation of state (‘ideal’, ‘isothermal)
`eos_gamma`	Adiabatic index
`eos_gm1`	gamma - 1
`eos_gm1_inv`	1 / (gamma - 1)
`eos_eint_fixed`	Fixed internal energy per unit mass (if isothermal)
Code units
`unit_l`	Code unit length in centimeters
`unit_t`	Code unit time in seconds
`unit_m`	Code unit mass in grams
Gravity
`gravity_mode`	Type of gravity (‘none’, ‘fixed_acc’, ‘fixed_pot’)
`gravity_g`	Gravitational acceleration (for mode ‘fixed_acc’)
`gravity_dir`	Direction of fixed acceleration (0 for x, 1 for y)
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

The constructor takes the following parameters:

Parameters:

hydro_scheme: HydroScheme: Container class for algorithmic choices

Methods

`addSourceTerms`(dt)	Add source terms to conserved quantities
`cflCondition`()	Compute the size of the next timestep
`conservedToPrimitive`(U, V)	Convert conserved to primitive variables
`emptyArray`([nq])	Get an empty array for fluid variables
`enforceBoundaryConditions`()	Enforce boundary conditions after changes
`fluxVector`(idir, V[, F])	Convert the flux vector F(V)
`limiterMC`(sL, sR, slim)	Monotonized-central limiter
`limiterMinMod`(sL, sR, slim)	Minimum-modulus limiter
`limiterNone`(sL, sR, slim)	Non-limiter (central derivative)
`limiterVanLeer`(sL, sR, slim)	The limiter of van Leer
`maxSpeedInDomain`()	Largest signal speed in domain
`primitiveEvolution`(idir, V, dV_dx)	Linear approximation of the Euler equations
`primitiveToConserved`(V, U)	Convert primitive to conserved variables
`primitiveToConservedRet`(V)	Convert primitive to new conserved array
`reconstructionConst`(idir, dt)	Piecewise-constant reconstruction
`reconstructionLinear`(idir, dt)	Piecewise-linear reconstruction
`riemannSolverHLL`(idir, VL, VR)	The HLL Riemann solver
`save`([filename])	Save the current state of a simulation
`setCodeUnits`([unit_l, unit_t, unit_m])	Define the meaning of the internal units
`setDomain`(nx, ny[, xmin, xmax, ymin, bc_type])	Set the physical and numerical size of the domain
`setEquationOfState`([eos_mode, gamma, eint_fixed])	Choose an equation of state
`setGravityMode`([gravity_mode, g, ...])	Add gravity to the simulation
`setGravityPotentials`()	Prepare gravitational potentials
`setHydroScheme`([hydro_scheme])	Set the hydro solver scheme
`soundSpeed`(V)	Sound speed
`timestep`([dt])	Advance the fluid state by a timestep dt
`xyGrid`()	Get a grid of the x and y cell center positions

setHydroScheme(hydro_scheme=None)

Set the hydro solver scheme

This function must be executed before the ulula.simulation.Simulation.setDomain() function.

Parameters:

hydro_scheme: HydroScheme: HydroScheme object that contains the settings for the hydro solver.

setEquationOfState(eos_mode='ideal', gamma=1.6666666666666667, eint_fixed=None)

Choose an equation of state

By default, an ideal gas EOS is chosen. If an isothermal EOS is chosen, the fixed temperature is given as the corresponding internal energy per unit mass. - conversion function - document parameters - write to file

This function must be executed before the ulula.simulation.Simulation.setDomain() function.

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.

setCodeUnits(unit_l=1.0, unit_t=1.0, unit_m=1.0)

Define the meaning of the internal units

The pure Euler equations (i.e., ignoring viscosity, cooling, etc.) are invariant under multiplications of up to three scale quantities, meaning that the solution of a hydro problem remains unchanged independent of what physical length, time, and mass scales the given numbers represent. One can alternatively think of rescalings in other, combined quantities such as density, pressure, and so on.

This function lets the user define the meaning of the internal length, time, and mass scales. The solution will not change unless the problem in question depends on physics beyond the pure Euler equations, such as gravity, cooling, and so on. However, plots of the solution will change (if a unit system other than code units is used).

The code units are given in cgs units. Some common units are defined in the utils module. For example, to set time units of years, unit_t = utils.units_t['yr']['in_cgs']. However, the code units can take on any positive, non-zero number chosen by the user.

This function must be executed before the ulula.simulation.Simulation.setDomain() function.

Parameters:

unit_l: float: Code unit for length in units of centimeters.
unit_t: float: Code unit for time in units of seconds.
unit_m: float: Code unit for mass in units of gram.

setGravityMode(gravity_mode='none', g=1.0, gravity_dir=1, compute_gradients=True)

Add gravity to the simulation

This function must be executed before the ulula.simulation.Simulation.setDomain() function. If the user chooses the fixed_acc mode, an acceleration g must be set, which acts in the negative x or y direction. The potential and gradients are computed automatically.

If the chosen mode is fixed_pot, the user must subsequently set the initial potential at the same time as the other initial conditions (in primitive variables).

Afterwards, the function setGravityPotentials() must be called to propagate the information. If the fixed_pot mode and compute_gradients are chosen, this function computes the spatial gradients of the user-defined potential. Otherwise, the user needs to set them manually. The latter can be more accurate if the analytical form of the gradients is known.

Parameters:

gravity_mode: str: The type of gravity to be added. Can be fixed_acc or fixed_pot.
g: float: If gravity_mode == 'fixed_acc', then g gives the constant acceleration in code units.
gravity_dir: int: The direction of a fixed acceleration, 0 meaning x and 1 meaning y. For a 1D simulation, the direction is forced to be x. For a 2D simulation, the direction is typically 1 (y) so that gravity points downwards.
compute_gradients: bool: If gravity_mode == 'fixed_pot', this parameter determines whether spatial gradients will be automatically computed or must be set by the user.

setDomain(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; must be at least 2.
ny: int: Number of grid points in y-direction; choosing ny = 1 leads to a 1D simulation.
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

emptyArray(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()

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)

setGravityPotentials()

Prepare gravitational potentials

This function must be executed after the ulula.simulation.Simulation.setDomain() function. If the gravity mode is fixed_acc, we automatically compute the potential. If it is fixed_pot, we expect that the user has set the potential and possibly the spatial gradients; if not, we compute them.

If the simulation is 1D, we interpret a constant acceleration as pointing to the negative x-direction, otherwise in the negative y-direction.

enforceBoundaryConditions()

Enforce boundary conditions after changes

This function fills the ghost cells with values from the physical domain to achieve certain behaviors. This function must be executed at each timestep. In particular:

Periodic: cells are rolled over from the other side of the domain so that it looks to the hydro solver as if the domain just continues on the other side.
Outflow: we take the value of the physical cells at the edge and copy them into the adjacent ghost cells, leading to flows that just continue across the edge.
Wall: the goal is to ensure that no mass or energy flux moves across the boundary. We achieve this condition by setting the ghost cells to a mirror image of the adjacent cells in the domain and inverting the perpendicular velocity, creating a counter-flow that balances any motion onto the edge.

primitiveToConserved(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(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(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(idir, V, F=None)

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(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(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()

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(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(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(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(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(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(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(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()

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

addSourceTerms(dt)

Add source terms to conserved quantities

This function implements adding the source terms S in dU / dt + div(F(U)) = S(U), namely we add the time-integrated source term to U. For a vector of conserved quantities

U = (rho, rho * ux, rho * uy, E)

the source term for gravity reads

S_grav = (0, -rho * dPhi/dx, -rho * dPhi/dy, rho * dPhi/dt).

For fixed-g gravity, we dPhi/dy = g and all other terms are zero. It might seem counter- intuitive that we are adding to the y-momentum and not to the energy. However, changes in momentum should be balanced by “falling,” i.e., by changes in the gravitational potential.

Parameters:

dt: float: The time over which the source term should be integrated.

timestep(dt=None)

Advance the fluid state by a timestep dt

This timestepping routine implements a dimensionally split scheme, meaning that we execute two sweeps in the two dimensions. We alternate the ordering of the sweeps with each timestep (xy-yx-xy and so on). This so-called Strang splitting maintains 2nd-order accuracy, but only if the timestep is the same for the two sweeps.

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.

The function internally handles the case of a CFL violation during the second sweep, which can occur even if the timestep was initially set to obey the CFL criterion. In this case, the returned timestep will be different from the input timestep (if given).

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. The user is responsible for ensuring that dt does not exceed the CFL criterion!

Returns:

dt: float: The timestep taken

save(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