Simulation framework

This page gives a brief introduction to the equations of hydrodynamics, lists the various hydro solver components that are implemented in Ulula, and documents the Ulula simulation class that constitutes the core of the code.

Overview

The Ulula simulation framework is implemented within a single class, called Simulation. The HydroScheme class contains information about the chosen algorithm. The user does not necessarily need to interact with the Simulation class directly because the simulation workflow is already implemented in the run() 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:
- Change the hydro scheme using setHydroScheme().
- Change the equation of state using setEquationOfState().
- Set the code units using setCodeUnits().
- Add gravity using setGravityMode().
- Add callback functions to manipulate the domain and boundary conditions at runtime using setUserUpdateFunction() and setUserBoundaryConditions().
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 (see Problem setups). The primitive variable fields are listed below.
If the physical setup contains gravity, execute the solveGravity() function with initial_setup=True.
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 v}) &= 0 \\ \frac{\partial (\rho {\pmb v})}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho {\pmb v} \otimes {\pmb v} + {\pmb I} P) &= -\rho \boldsymbol{\nabla} \Phi \\ \frac{\partial E}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} ([E + P] {\pmb v}) &= \rho \frac{\partial \Phi}{\partial t} \\ \end{align}\end{split}\]

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

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

Here \(\epsilon\) is the internal energy per unit mass. How this energy 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 fact that the right-hand-side of some of the conservation laws above is not always zero means that some quantities may not always be conserved. In particular, in the presence of a gravitational potential \(\Phi\), momentum and energy are added to the gas due to gravity, which itself depends on density via the Poisson equation,

\[\nabla^2 \Phi = 4 \pi G \rho \,.\]

However, we can imagine cases where we want to substitute a different potential for the Poisson equation, for example, if we want to mimic Earth’s essentially fixed gravitational acceleration without modeling Earth itself. Ulula can thus also accommodate fixed, user-defined potentials or accelerations.

To understand how the equations above are solved in Ulula (and in 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 \({\pmb S}\) the source vector,

\[\begin{split}{\pmb U} \equiv \left(\begin{array}{c} \rho \\ \rho v_{\rm x} \\ \rho v_{\rm y} \\ E \end{array}\right) \qquad \mathcal{F}({\pmb U}) \equiv \left(\begin{array}{c} \rho \pmb{v} \\ \rho v_{\rm x} {\pmb v} + \delta_{xj} P \\ \rho v_{\rm y} {\pmb v} + \delta_{yj} P \\ (E + P) {\pmb v} \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_{ij}\) means that \(P\) is only added to the x-component of the flux vector for \(v_{\rm x}\) and to the y-component of the flux vector for \(v_{\rm y}\). 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. When computing those updates from the flux and source vectors, we need to know pressure and velocity in addition to the conserved variables. We thus frequently convert to and from so-called primitive variables,

\[\begin{split}{\pmb V} \equiv \left(\begin{array}{c} \rho \\ v_{\rm x} \\ v_{\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\)
\(v_{\rm x}\)	`VX`	X-velocity	\(l/t\)
\(v_{\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 v_{\rm x}\)	`MX`	X-momentum density	\(m/l^2/t\)
\(\rho v_{\rm y}\)	`MY`	Y-momentum density	\(m/l^2/t\)
\(E\)	`ET`	Total energy density	\(m/l/t^2\)
Gravity (depending on user choices)
\(\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\)
\(\Phi_{\rm user}\)	`GU`	User-defined fixed potential	\(l^2/t^2\)

The abbreviations are used when finding the position of each variable in the \({\pmb U}\) and \({\pmb V}\) arrays, when selecting variables in the Plotting & Analysis module, and when setting the initial conditions. The gravity-related fields exist only in the primitive variable array.

The equations solved by Ulula are invariant under unit transformations, meaning that the length, mass, and time units that make up the composite units listed above can be set to arbitrary values, so-called “code units” (using the setCodeUnits() function). Changing these units does not change the outcome of the simulation, but it will be reflected in plots (unless you choose to plot in code units). Code units can take on any value, but some common units are listed in the units module.

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 1st-order (“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 only.
`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 (MC). 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.
`hllc`	The HLLC Riemann solver, which corresponds to HLL plus contact discontinuities. This solver cannot be used with an isothermal EOS, which makes sense since an isothermal gas admits no 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 almost identical to the primitive formulation, 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 the default.

The purpose of Ulula is to experiment with hydro algorithms, including those that frequently fail. Thus, the code allows for sub-optimal (aka. crazy) combinations of the components 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.core.simulation.HydroScheme(reconstruction='linear', limiter='mc', riemann='hllc', 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 table for valid choices
limiter: string: Slope limiter algorithm; see table for valid choices
riemann: string: Riemann solver; see table for valid choices
time_integration: string: Time integration scheme; see table 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 in the first sweep 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. 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 may 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.core.simulation.Simulation

The Ulula hydro solver framework

This class contains all simulation data and routines. The internal fields have the following meaning:

Field	Meaning
Domain
`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
`nghost`	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)
`x`	Array of x values at cell centers (dimensions [nx + 2 ng])
`y`	Array of y values at cell centers (dimensions [ny + 2 ng])
`bc_type`	Type of boundary condition (‘periodic’, ‘outflow’, ‘wall’)
`domain_set`	Once the domain is set, numerous settings cannot be changed any more
Fluid variables
`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.)
`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
`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)
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_mu`	Mean particle mass in proton masses (can be None)
`eos_eint_fixed`	Fixed internal energy per unit mass (if isothermal)
`eos_pressure_floor`	Optional minimum pressure
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’, ‘poisson’, ‘fixed+poisson’)
`gravity_dir`	Direction of fixed acceleration (0 for x, 1 for y)
`gravity_g`	Gravitational acceleration (for mode ‘fixed_acc’)
`gravity_G`	Gravitational constant (for mode ‘poisson’ or ‘fixed+poisson’)
User-defined boundary and update functions
`do_user_updates`	True if the user has supplied a custom domain update function
`do_user_bcs`	True if the user has supplied a custom boundary condition function
`user_update_func`	A function to be called before the boundary conditions are enforced
`user_bc_func`	A function to be called after the boundary conditions are enforced

Methods

`addSourceTerms`(dt)	Add source terms to conserved quantities
`cflCondition`()	Compute the size of the next timestep
`checkPressure`(V)	Check for zero or negative pressures
`conservedToPrimitive`(U, V)	Convert conserved to primitive variables
`enforceBoundaryConditions`([a, slc_q])	Enforce boundary conditions after changes
`fluxVector`(idir, V)	Compute the flux vector F(V)
`limiterMC`(sL, sR, slim)	Monotonized-central limiter
`limiterMinMod`(sL, sR, slim)	Minimum-modulus limiter
`limiterNone`(sL, sR, slim)	No 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
`riemannSolverHLLC`(idir, VL, VR)	The HLLC 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 domain
`setEquationOfState`([eos_mode, gamma, mu, ...])	Choose an equation of state
`setGravityMode`([gravity_mode, gravity_dir, g, G])	Add gravity to the simulation
`setHydroScheme`([hydro_scheme])	Set the hydro solver scheme
`setUserBoundaryConditions`([user_bc_func])	Set user-defined boundary conditions
`setUserUpdateFunction`([user_update_func])	Set user-defined updates in the domain
`solveGravity`([initial_setup])	Compute gravitational potentials
`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 setDomain() function.

Parameters:

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

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

Choose an equation of state

The equation of state (EOS) captures the microphysics of the gas that is being simulated. This function must be executed before the setDomain() function.

The default is an ideal gas EOS with \(\gamma = 5/3\). If an isothermal EOS is chosen, the fixed temperature must be specified as an internal energy per unit mass, which can be computed from temperature using the internalEnergyFromTemperature() function.

Parameters:

eos_mode: str: Can be 'ideal' or 'isothermal'.
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.
mu: float: Mean particle mass in units of proton masses. Must be specified if temperature is to be plotted, regardless of the EOS.
eint_fixed: float: A fixed internal energy per unit mass in code units for an isothermal EOS, ignored otherwise.
pressure_floor: float: If not None, this pressure (in code units) is enforced as a minimum. Such a pressure floor can be helpful in simulations where the pressure becomes negative due to numerical errors. However, energy is no longer conserved since the pressure floor may effectively add thermal energy.

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. As a result, the values of the code units are not used at all in the simulation class. 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 units module. For example, to set time units of years, unit_t = units.units_t['yr']['in_cgs']. However, the code units can take on any positive, non-zero value chosen by the user.

This function must be executed before the 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', gravity_dir=1, g=1.0, G=1.0)

Add gravity to the simulation

If gravity is to be part of the setup, this function must be executed before the setDomain() function. Once the domain is set, the function solveGravity() must be called with 'initial_setup = True' to initialize potential and gradients. The Setup class performs this step automatically. The user has a choice between a number of gravity modes.

If 'fixed_acc' is chosen, an acceleration g is uniformly applied. If the simulation is 1D, we interpret a constant acceleration as pointing to the negative x-direction, otherwise in the negative y-direction.

If the chosen mode is 'fixed_pot', the user must subsequently set the initial potential at the same time as the other initial conditions. Ulula will take the gradients of this potential and apply them uniformly.

If 'poisson' is chosen, the Poisson equation is solved at each timestep so that the density field in the simulation generates a gravitational field with a gravitational constant G. This mode works only with periodic boundary conditions.

The 'fixed+poisson' mode is the same as 'poisson', but the user additionally provides a fixed potential as in the 'fixed_pot' mode. This potential is added to the Poisson potential at each timestep.

Parameters:

gravity_mode: str: The type of gravity to be added. Can be 'fixed_acc', 'fixed_pot', 'poisson', or 'fixed+poisson'.
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.
g: float: If gravity_mode is 'fixed_acc', then g gives the constant acceleration in code units.
G: float: If gravity_mode is 'poisson' or 'fixed+poisson', then G is the gravitational constant in code units.

setUserUpdateFunction(user_update_func=None)

Set user-defined updates in the domain

If a function is passed for user_update_func, that function will be called with the simulation object as an argument before boundary conditions are enforced. This mechanism allows the phyiscal setup to influence the simulation at runtime, while the boundary conditions are still automatically enforced. The Setup base class contains such a function, which becomes active if it is overwritten by a child class.

The code does not check that what a user update function does makes any sense! If the function implements unphysical behavior, an unphysical simulation will result. For example, the update function must ensure that primitive and conserved variables are consistent after the function call.

Parameters:

user_update_func: func: Function pointer that takes the simulation object as an argument.

setUserBoundaryConditions(user_bc_func=None)

Set user-defined boundary conditions

If a function is passed for user_bc_func, that function will be called with the simulation object as an argument every time after boundary conditions have been enforced. This mechanism allows for custom boundary conditions such as a flow entering the domain from a certain side. The Setup base class contains such a function, which becomes active if it is overwritten by a child class.

The code does not check that what a user update function does makes any sense! If the function implements unphysical behavior, an unphysical simulation will result. For example, the update function must ensure that primitive and conserved variables are consistent after the function call.

Parameters:

user_bc_func: func: Function pointer that takes the simulation object as an argument.

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

Set the physical and numerical domain

This function creates the memory structure for the simulation as well as pre-computed slices that index the variable 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 (in code units).
xmax: float: Right edge in physical coordinates (in code units).
ymin: float: Bottom edge in physical coordinates (in code units).
ymax: float: Top edge in physical coordinates (in code units).
bc_type: string: Type of boundary conditions, which can be 'periodic', 'outflow', or 'wall'. With periodic boundary conditions, the domain wraps around, meaning that flows out of the right edge re-enter on the left and so on. Outflow means that flows at the domain edge continue. Be careful with ouflow BCs! They do not conserved mass, energy, or momentum, and they can lead to weird effects if the flow is not actually out of the domain. For example, if there is a flow into the domain at an edge, this flow will continue adding mass indefinitely unless it is stopped by counter-acting pressures (which may or may not be desired). Finally, wall boundary conditions reflect the fluid flow. They conserve mass and energy but not momentum.

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.

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(a=None, slc_q=None)

Enforce boundary conditions after changes

This function fills the ghost cells with values from the physical domain to achieve certain behaviors. It is 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.

Parameters:

a: array_like: List of arrays to set boundary conditions in. If None, the function sets BCs for the hydro variables in the default primitive and conserved arrays. If not None, BCs are set in the given array. In that case, the user-defined BC and update functions are not executed.
slc_q: slice: If a is not None, this slice chooses which variables the BCs are set for.

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_all.
U: array_like: Output array of conserved fluid variables.

primitiveToConservedRet(V)

Convert primitive to new conserved array

Same as primitiveToConserved(), but creating a new array for the output conserved variables.

Parameters:

V: array_like: Input array of primitive fluid variables with first dimension nq_all.

Returns:

U: array_like: Array of fluid variables with first dimension nq_hydro and otherwise same dimensions as the 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, since in some cases, conversions need to be performed on other fluid variables. Note that gravitational potentials are stored only in the primitive variable vector, so that V must contain the correct potential even if U is technically the input.

Parameters:

U: array_like: Input array of conserved fluid variables with first dimension nq_hydro.
V: array_like: Output array of primitive fluid variables with first dimension nq_all, which must already contain the gravitational potential if gravity is on.

checkPressure(V)

Check for zero or negative pressures

Physically, pressure can never be negative. Numerically, however, negative pressures can when we compute pressure from conserved variables (by subtracting kinetic energy from total energy) or when we evolve it otherwise. After such operations, this function is called to either raise an error or enforce a user-defined pressure floor.

Parameters:

V: array_like: Input array of primitive fluid variables with first dimension nq_all.

fluxVector(idir, V)

Compute 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_all.

Returns:

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

primitiveEvolution(idir, V, dV_dx)

Linear approximation of the Euler equations

Instead of the usual conservation-law form, we can also think of the Euler equations in 1D as a change \(\partial {\pmb V}/ \partial t + A_x({\pmb V})\ \partial {\pmb V} / \partial x = {\pmb S}\). This function returns \(\Delta {\pmb V} / \Delta t\) given an input state \({\pmb V}\) and a vector of spatial derivatives \(\Delta {\pmb 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 at least nq_hydro.
dV_dx: array_like: Array of spatial derivatives of the fluid variables with first dimension at least nq_hydro.

Returns:

dV_dt: array_like: Array of linear approximations 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.

Returns:

cs: array_like: Array of sound speed with same dimensions as input array (for a single variable).

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.

Returns:

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

reconstructionConst(idir, dt)

Piecewise-constant reconstruction

The term “piecewise-constant” means that the fluid state at the cell center is interpreted as representing the entire cell, meaning that we perform no reconstruction. If this option is chosen, the left/right cell edge arrays already point to the cell-centered values so that this function does nothing at all. It merely 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. This interpolation leads to 2nd-order convergence in space.

If the Hancock time integration scheme is chosen, the reconstructed edge states are also advanced by half a timestep to get 2nd-order convergence in time. 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)

No 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. Must have same dimensions as sL.
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. Must have same dimensions as sL.
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. Must have same dimensions as sL.
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. Must have same dimensions as sL.
slim: array_like: Output array of limited slope. Must have same dimensions as sL and sR.

riemannSolverHLL(idir, VL, VR)

The HLL Riemann solver

A 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 to the left and to the right, but it takes into account only a single possible intermediate state, which ignores contact discontinuities.

Parameters:

idir: int: Direction of sweep (0 = x, 1 = y).
VL: array_like: Array of primitive fluid variables on the left sides of the interfaces.
VR: array_like: Array of primitive fluid variables on the right sides of the interfaces. Must have the same dimensions as VL.

Returns:

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

riemannSolverHLLC(idir, VL, VR)

The HLLC Riemann solver

Similar to the HLL Riemann solver, but with an additional distinction of whether the interfaces lie to the left or right of contact discontinuities. The implementation follows Chapter 10.4 in Toro 2009.

This Riemann solver explicitly uses pressure and total energy in its calculations and is thus not compatible with an isothermal EOS where we do not track those variables. This is, in some sense, by construction: in an isothermal gas, there are no contact discontinuities, and the HLLC solver yields no advantage over HLL.

Parameters:

idir: int: Direction of sweep (0 = x, 1 = y).
VL: array_like: Array of primitive fluid variables on the left sides of the interfaces.
VR: array_like: Array of primitive fluid variables on the right sides of the interfaces. Must have the same dimensions as VL.

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.

Returns:

dt: float: The next timestep.

solveGravity(initial_setup=False)

Compute gravitational potentials

This function computes the gravitational potential from density (if necessary) and the gradients of the potential. It must be executed once after the initial setup, i.e., after setDomain() but before the simulation is run. The Setup class performs this step automatically. In this first step, the function computes the potential that corresponds to a fixed acceleration (for 'fixed_acc' gravity), or the fixed gradients (for 'fixed_pot' gravity), or the initial potential and gradients from the density field (for 'poisson' gravity). In the latter case, the calculation needs to be repeated each time the gravity source term is added to the conserved fluid variables.

The Poisson equation is solved via Fourier transform because the Laplacian reduces to a multiplication by a “Green’s function” in Fourier space (which only works with periodic BCs due to the periodic nature of the sine/cosine basis functions). The Green’s function is specific to the discretization of the Laplacian derivative. We assume a standard discretized second derivative, \(\partial^2 \phi / \partial x^2 \approx (\phi_{i-1} - 2 \phi_i + \phi_{i+1}) / 2 \Delta x\). See e.g. Section 9.2 in the book by Mike Zingale for a derivation.

Parameters:

initial_setup: bool: Whether this is the first time the function is called.

addSourceTerms(dt)

Add source terms to conserved quantities

This function implements the source terms in the equations above, meaning that it adds \({\pmb U} \rightarrow {\pmb U} + {\pmb S} \Delta t\) to the conserved fluid variables. For fixed-acceleration gravity, we only set add the acceleration in the appropriate dimension. It might seem counter-intuitive that we are adding to the momentum but not to the energy. However, changes in momentum are balanced by “falling,” i.e., by changes in the gravitational potential. If the gravitational potential is changing (namely, for 'poisson' gravity), we also add the change in potential to the total energy.

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.

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).

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. Any source terms are added in two half-timesteps before and after the hydrodynamical sweeps.

Parameters:

dt: float: Size of the timestep to be taken. If None, the timestep is computed based on the CFL condition using the cflCondition() function. This timestep should be used in most circumstances, but sometimes we wish to manually set a smaller timestep, for example, to output a file or plot at a particular time. The user is responsible for ensuring that a manually chosen dt does not exceed the CFL criterion!

Returns:

dt: float: The timestep taken.

save(filename=None)

Save the current state of a simulation

Ulula uses the hdf5 format so save snapshots 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. Similarly, the attributes of the domain, eos, gravity groups and so on have the same meaning as in the Simulation object. The code group contains the version of Ulula that was used to create the file.

The grid group contains the 2D grid of the primitive fluid variables 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.core.simulation.load(filename, user_update_func=None, user_bc_func=None)

Load a snapshot file into a simulation object

A simulation can be exactly restored from a file, with the exception of user-defined function pointers which cannot be saved to file. If such pointers were given to the original simulation, the same functions must be passed to the loaded simulation object.

Parameters:

filename: str: Input filename.
user_update_func: func: Function pointer that takes the simulation object as an argument. See setUserUpdateFunction().
user_bc_func: func: Function pointer that takes the simulation object as an argument. See setUserBoundaryConditions().

Returns:

sim: Simulation: Object of type Simulation.