Function Reference

DefaultSimulationParameters()

Return a set of default parameters for a simulation.

Returns

A Dict{String, Any} with default parameters for a simulation.

Notes

• Keys with a value of nothing must be set before using this dictionary to initialize a simulation. This is done instead of providing a default function to avoid the possibility of accidentally initializing a simulation with unintended initial conditions.

Parameters

• start_time: The initial time of the simulation. (Unit: s; Default: 0.0)
• start_position: The position of the left side of the domain. (Unit: m; Default: 0.0)
• end_position: The position of the right side of the domain. (Unit: m; Default: 1.0)
• number_of_zones: The number of zones to divide the domain into. (Unit: ⋅; Default: 1000)
• CFL: The CFL number to use. (Unit: ⋅; Default: 0.2)
• artificial_viscosity_coefficient: The coefficient to use for artificial viscosity. (Unit: ⋅; Default: 1e0)
• artificial_conductivity_coefficient: The coefficient to use for artificial conductivity. (Unit: ⋅; Default: 1e-2)
• minimum_timestep: The minimum allowable timestep size. Simulation will halt if timestep falls below this value. (Unit: s; Default: 1e-7)
• maximum_cycles: The maximum number of cycles to perform. Simulation will halt if more than this many timesteps are taken. (Unit: ⋅; Default: 1e6)
• init_density_function: A Function that returns the initial density as a function of position x. (Unit: kg/m³; Default: nothing, must be user-supplied)
• init_velocity_function: A Function that returns the initial velocity as a function of position x. (Unit: m/s; Default: nothing, must be user-supplied)
• init_pressure_function: A Function that returns the initial pressure as a function of position x. (Unit: N/m²; Default: nothing, must be user-supplied)
• init_gamma_function: A Function that returns the ratio of specific heats as a function of position x. (Unit: ⋅; Default: nothing, must be user-supplied)

InitializeSimulation( parameters::Dict{String, Any} )

Initialize a simulation with the parameters provided by parameters. A default set of parameters can be returned from the DefaultSimulationParameters() function.

Returns

A struct of type Simulation{T} that describes the current simulation state. The type T is determined by the type of the members of parameters.

Notes

• This function performs some basic sanity checking on input parameters such as ensuring that end_time > start_time or end_position > start_position, and will raise an error if any of these checks fail.
• A warning will be raised if any keys in the dictionary are unused, but this will not halt execution.
• This function will not warn about any potential issues with problem configuration or potential stability issues.

UpdateSimulationState!( state::Simulation{T}, gamma::Function, density::Function, velocity::Function, pressure::Function ) where { T <: AbstractFloat }

Update the simulation state using the supplied functions to define a new state. See Notes for information about the expected signature of the functions.

Returns

nothing. Updates the state input in-place.

Arguments

• state: A Simulation{T} representing the simulation state to be updated
• gamma: A Function that returns the new values of the ratio of specific heats, gamma
• density: A Function that returns the new values of the density
• velocity: A Function that returns the new values of velocity
• pressure: A Function that returns the new values of pressure

Notes

Unlike the functions used in initial problem setup, the functions supplied to UpdateSimulationState have a slightly different expected signature of ExampleFunction( x::T, oldValue::T ) where { T <: AbstractFloat } where oldValue will be the current value of the state variable at current position x (e.g., the previous value of pressure).

Side Effects

• Updates the values stored in the vectors for state.gamma, state.mass, state.velocity, state.intenergy, and fields derived from the equation of state in-place.

ConfigureSimulationCallbacks( simulation::Simulation{T} ) where { T <: AbstractFloat }

Configure a SimulationCallback structure for setting up callbacks as part of a simulation run.

Returns

• A SimulationCallback structure

Parameters

• simulation: A Simulation{T} describing the simulation that the callbacks are intended for.

Notes

• Currently, the simulation parameter has no effect in this function and the resulting SimulationCallback structure can be used by any simulation, not just the one passed as an argument. However, this may change in future version.
• Callback entries can be added using RegisterCycleCallback!(), RegisterTimeCallback!(), or RegisterTimeDeltaCallback!().

RegisterCycleCallback!( callbacks::SimulationCallback, func::Function, N::UInt; initial_cycle::UInt=0 )

Register a function func to be executed every N cycles of the simulation starting at cycle initial_cycle. Callbacks are stored in the callbacks input parameter and passed as an argument to the timestepping routines

Returns

nothing, mutates the callbacks parameter

Parameters

• callbacks: A SimulationCallback that describes the callbacks to be executed
• func: A Function that should be called by this callback. See Notes for information on expected function signature
• N: A UInt indicating how many cycles should occur between callbacks. Must be N ≥ 1
• initial_cycle: A UInt indicating the cycle number to start calling this callback at (Default: 0)

Notes

• func is expected to accept arguments of the form CallbackFunction( arg::Simulation{T} ) where { T <: AbstractFloat }. arg will be the simulation state at the end of the cycle that the callback is executed on.
• This type of callback is useful for cases where a callback should be called at a fixed number of cycles between each callback.

RegisterTimeCallback!( callbacks::SimulationCallback, func::Function, times::Vector{T} ) where { T <: AbstractFloat }

Register a function func to be executed at the list of times specified in times. Callbacks are stored in the callbacks input parameter and passed as an argument to the timestepping routines

Returns

nothing, mutates the callbacks parameter

Parameters

• callbacks: A SimulationCallback that describes the callbacks to be executed
• func: A Function that should be called by this callback. See Notes for information on expected function signature
• times: A Vector{T} indicating the times at which to execute the callback

Notes

• func is expected to accept arguments of the form CallbackFunction( arg::Simulation{T} ) where { T <: AbstractFloat }. arg will be the simulation state at the end of the cycle that the callback is executed on.
• Callbacks are executed at the end of the first cycle where the simulation time is greater than the next element in times. As a result, it is not guaranteed that the callback will be called at exactly the time specified in the times vector.
• times is sorted into ascending order before being stored.
• This type of callback is useful if the callback should be called at an irregular series of times. See CycleCallback if the callback should be called at a regular number of cycles, or TimeDeltaCallback if the callback should be called at a fixed temporal cadence.

RegisterTimeDeltaCallback!( callbacks::SimulationCallback, func::Function, delta::T, initial_time::T=T(0.0) ) where { T <: AbstractFloat }

Register a function func to be executed every delta seconds starting at initial_time. Callbacks are stored in the callbacks input parameter and passed as an argument to the timestepping routines

Returns

nothing, mutates the callbacks parameter

Parameters

• callbacks: A SimulationCallback that describes the callbacks to be executed
• func: A Function that should be called by this callback. See Notes for information on expected function signature
• delta: A scalar T indicating how often to execute the callback
• inital_time: A scalar T indicating the time of the first callback to execute. (Default: 0.0)

Notes

• func is expected to accept arguments of the form CallbackFunction( arg::Simulation{T} ) where { T <: AbstractFloat }. arg will be the simulation state at the end of the cycle that the callback is executed on.
• Callbacks are executed at the end of the first cycle where the simulation time is greater than the last time the callback was executed plus delta. As a result, it is not guaranteed that the callback will be called at exactly delta seconds since the last call.
• The time of the next callback is computed relative to the expected time of the current callback, not when the callback was actually called. In other words, if a callback that should be called at t₀ was actually called at t₀+ϵ due to the timestep size, the next callback will be scheduled for t₀+δ, not t₀+ϵ+δ.
• This type of callback is useful for cases where a callback should be called at a fixed temporal spacing between each callback.

AdvanceToTime( state::Simulation{T}, stoptime::T; timestep::Union{T,Nothing}=nothing, exact::Bool=false, callbacks::Union{SimulationCallback,Nothing}=nothing ) where { T <: AbstractFloat }

Advance the simulation with initial state state to a time specified by stoptime with a fixed timestep.

Returns

A Simulation{T} representing the state at the end of the final cycle.

Arguments

• input: A Simulation{T} representing the simulation state at the start of the first cycle.
• stoptime: The time to advance to. (Unit: s)
• timestep: The timestep size, or nothing to use adaptive timestepping. (Unit: s, Default: nothing)
• exact: If true, try to stop as close as possible to stoptime by adjusting the final timestep size (Default: false)
• callbacks: An optional SimulationCallback() structure containing information about callbacks to be performed during the simulation.

Notes

• The current simulation time is determined by the Simulation field state.time. If state.time > stoptime, no steps will be taken.
• This function allocates two deepcopy()s of the input state and returns the copy corresponding to the final state.
• This function simply calls AdvanceOneCycle! repeatedly until the simulation time reaches stoptime. The primary advantage to using this function as opposed to AdvanceOneCycle() or AdvanceNCycles() is that various backing arrays are pre-allocated to improve speed.
• If exact=true, the timestep of the final cycle is adjusted so that the time of the final state is as close as possible to stoptime.
• If using adaptive timestepping, the timestep size is determined based on the minimum time for an acoustic wave to traverse a zone. See CalculateTimestepSize() for further details.

AdvanceOneCycle( state::Simulation{T}; timestep::Union{T,Nothing}=nothing, callbacks::Union{SimulationCallback,Nothing}=nothing ) where { T <: AbstractFloat }

Advance the simulation by one cycle.

Returns

A Simulation{T} representing the state at the end of the cycle

Arguments

• input: A Simulation{T} representing the simulation state at the start of the cycle.
• timestep: The size of the time step, or nothing to use adaptive timestepping. (Unit: s, Default: nothing)
• callbacks: An optional SimulationCallback() structure containing information about callbacks to be performed during the simulation.

Notes

• This function allocates a deepcopy() of the input state. The copy is modified and returned from this function.
• If using adaptive timestepping, the timestep size is determined based on the minimum time for an acoustic wave to traverse a zone. See CalculateTimestepSize() for further details.

AdvanceOneCycle!( output::Simulation{T}, input::Simulation{T}; timestep::Union{T,Nothing}=nothing, callbacks::Union{SimulationCallback,Nothing}=nothing ) where { T <: AbstractFloat }

Advance the simulation by one cycle.

Returns

nothing. Modifies output in-place.

Arguments

• output: A Simulation{T} that will represent the output state. This will be modified by the function to represent the simulation state after advancing one cycle.
• input: A Simulation{T} that represents the simulation state at the start of the cycle.
• timestep: The size of the time step, or nothing to use adaptive timestepping. (Unit: s, Default: nothing)
• callbacks: An optional SimulationCallback() structure containing information about callbacks to be performed during the simulation.

Side Effects

• All fields of output are modified in-place.
• If using adaptive timestepping, the timestep size is determined based on the minimum time for an acoustic wave to traverse a zone. See CalculateTimestepSize() for further details.

AdvanceNCycles( state::Simulation{T}, ncycles::UInt; timestep::Union{T,Nothing}=nothing, callbacks::Union{SimulationCallback,Nothing}=nothing ) where { T <: AbstractFloat }

Advance the simulation by ncycles cycles.

Returns

• A Simulation{T} representing the state at the end of the final cycle

Arguments

• input: A Simulation{T} representing the simulation state at the start of the first cycle.
• ncycles: The number of cycles to advance.
• timestep: The size of the time step, or nothing to use adaptive timestepping. (Unit: s, Default: nothing)
• callbacks: An optional SimulationCallback() structure containing information about callbacks to be performed during the simulation.

Notes

• This function allocates two deepcopy()s of the input state and returns the copy corresponding to the final state.
• This function calls AdvanceOneCycle!() a total of ncycles times to advance the simulation. The primary advantage to using this function as opposed to AdvanceOneCycle() if the number of cycles to advance is known is that various backing arrays are pre-allocated to improve speed.
• If using adaptive timestepping, the timestep size is determined based on the minimum time for an acoustic wave to traverse a zone. See CalculateTimestepSize() for further details.

CalculateTimestepSize( state::Simulation{T} ) where { T <: AbstractFloat }

Compute an automatic timestep size for the next simulation cycle based on the current simulation state. See the Notes section for details on how the timestep is determined.

Returns

• A scalar of type T representing the timestep size for the next cycle based on the current simulation state.

Arguments

• state: A Simulation{T} representing the problem state.

Notes

• For each zone, the local speed of sound is computed according to c = √( γ P / ρ ), where γ, P, and ρ are the ratio of specific heats, the pressure, and the density of the gas in that zone.
• The time for an acoustic wave to traverse a zone with length Δx is computed as t = Δx / c.
• The minimum traversal time for all zones is multiplied by the user-specified CFL number to obtain the timestep size.
• Sanity checking for negative zone sizes and small timesteps is performed to detect problem instability.

EOS_Density( mass::T, Δx::T ) where { T <: AbstractFloat } = mass / Δx

Compute the density of the fluid in a zone.

Returns

A scalar of type T representing the density of a zone.

Arguments

• mass: The total mass contained within the zone. (Unit: kg)
• Δx: The size of the zone. (Unit: m)

Notes

Density is calculated as:

ρ = mass / Δx

As this is the physical definition of density, this calculation does not assume any particular equation of state.

EOS_Pressure( γ::T, ρ::T, e::T ) where { T <: AbstractFloat }
EOS_Pressure( γ::T, mass::T, Δx::T, e::T ) where { T <: AbstractFloat }

Compute the pressure of a zone using an ideal gas equation of state.

Returns

A scalar of type T representing the pressure within the zone.

Arguments

• γ: The ratio of specific heats of the fluid in the zone. (Unit: ⋅)
• ρ: The density of the zone. (Unit: kg/m³)
• mass: The total mass contained within the zone. (Unit: kg)
• Δx: The size of the zone. (Unit: m)
• e: The internal energy per unit mass of the zone. (Unit: m²/s²)

Notes

The pressure is calculated as:

P = ( γ - 1 ) ⋅ ρ ⋅ e

The four-parameter version of this function computes density using EOS_Density(). See the documentation for that function for further details.

EOS_SpeedOfSound( γ::T, P::T, ρ::T ) where { T <: AbstractFloat }
EOS_SpeedOfSound( γ::T, e::T, mass::T, Δx::T ) where { T <: AbstractFloat }

Compute the speed of sound in a zone using an ideal gas equation of state.

Returns

A scalar of type T representing the speed of sound in the zone.

Arguments

• γ: The ratio of specific heats in the zone. (Unit: ⋅)
• e: The internal energy per unit mass in the zone. (Unit: m²/s²)
• P: The pressure in the zone. (Unit: kg/(m⋅s²))
• ρ: The density of the fluid in the zone. (Unit: kg/m³)
• mass: The mass contained in the zone. (Unit: kg)
• Δx: The length of the zone. (Unit: m)

Notes

The speed of sound is calculated as:

c = √( γ * P / ρ )

The four-parameter version of this function computes density using EOS_Density() and pressure using EOS_Pressure(). See the documentation of those functions for further details.

artificial_viscosity( Cᵥ::T, c::T, ρ::T, Δx::T, u₋::T, u₊::T ) where { T <: AbstractFloat }

Compute an artificial viscosity within a zone.

Returns

A value of type T representing the value of the artificial viscosity.

Arguments

• Cᵥ: An O(1) coefficient to control the strength of the artificial viscosity. (Unit: ⋅)
• c: The speed of sound in the zone. (Unit: m/s)
• ρ: The density of the zone. (Unit: kg/m³)
• Δx: The length of the zone. (Unit: m)
• u: The velocity of the zone boundaries, with superscripts - and + referring to the left and right boundaries of the zone, respectively. (Unit: m/s)

Notes

• This artificial viscosity is based on the method described by Wilkins (1980), which in turn relies upon the methods of Von Neumann and Richtmyer (1950) and Landschoff (1955). The values computed by this function should be added to the pressure field during governing equation updates.
• This function returns zero if ∂u/∂x > 0, which will be the case for regions where the flow is expanding. This is done to restrict artificial viscosity only to regions of compression.

artificial_conductivity( Cₖ::T, u::T, c₋::T, e₋::T, Δx₋::T, c₊::T, e₊::T, Δx₊::T ) where { T <: AbstractFloat }

Compute an artificial flux of internal energy across a zone interface.

Returns

A value of type T representing the artifical flux of internal energy across a zone boundary.

Arguments

• Cₖ: An O(1) coefficient to control the strength of the artificial conductivity. (Unit: ⋅)
• u: The velocity of the zone interface. (Unit: m/s)
• c: The speed of sound, where superscript - and + refer to the zones to the left and right of the zone interface, respectively. (Unit: m/s)
• e: The internal energy per unit mass. Superscript - and + refer to the zones to the left and right of the zone interface, respectively. (Unit: m²/s²)
• Δx: Length of the zone. Superscript - and + refer to the zones to the left and right of the zone interface, respectively. (Unit: m)

Notes

The artificial conductivity is modeled as a Fickian diffusivity. That is, the flux of energy across a zone boundary, fₑ, is described by

fₑ = -κ ∂e/∂x

where e is the internal energy per unit mass and κ is the (artificial) conductivity coefficient. The gradient of internal energy is treated with a simple forward finite difference. The artificial conductivity coefficient is modeled as

κ = Cₖ * cₘ * Δx

where

• Cₖ an O(1) coefficient.
• cₘ is a characteristic velocity taken to be max(c̄±u, c̄), where c̄=0.5*(c₋+c₊) is the average speed of sound of the two adjacent zones and u is the velocity of the zone interface.
• Δx is the distance between the zone centers

By convention, this leads to a positive flux if energy is diffusing in the positive x direction, and negative if it is diffusing in the negative x direction.

struct Simulation{T}

A structure containing all the internal variables and arrays used in the simulation.

Parameters

• nzones::Int: The number of zones in the simulation. (Unit: ⋅)
• nedges::Int: The number of zone edges in the simulation, equal to nzones + 1. (Unit: ⋅)
• CFL::Float64: The CFL number to be used when calculating timesteps. (Unit: ⋅)
• start_time::Float64: The initial time of the simulation. (Unit: s)
• viscosity_coefficient::Float64: The coefficient used to scale artificial viscosity. (Unit: ⋅)
• conductivity_coefficient::Float64: The coefficient used to scale artificial conductivity. (Unit: ⋅)
• time::Base.RefValue{T}: The current time of the simulation. (Unit: s)
• dt::Base.RefValue{T}: The size of the timestep taken in the last cycle. (Unit: s)
• cycles::Base.RefValue{UInt}: The number of cycles performed so far. (Unit: ⋅)
• min_dt::Float64: The minimum allowable timestep size. Simulation will halt if Δt falls below this value. (Unit: s)
• max_cycles::UInt: The maximum number of cycles to perform. Simulation will halt if cycles exceeds this value. (Unit: ⋅)
• zone_edge::Vector{T}: A vector of locations of zone edges. Increases monotonically. (Unit: m)
• zone_center::Vector{T}: A vector of locations of zone centers, defined as the midpoint between two zone edges. Increases monotonically. (Unit: m)
• zone_length::Vector{T}: A vector of the length of each zone. (Unit: m)
• gamma::Vector{T}: A vector of the ratio of specific heats inside each zone. (Unit: ⋅)
• mass::Vector{T}: A vector of the mass contained within each zone. Assumed constant. (Unit: kg)
• density::Vector{T}: A vector of the density of the fluid within each zone. Computed using the equation of state as mass/Δx. (Unit: kg/m³)
• velocity::Vector{T}: A vector of velocities of zone edges. (Unit: m/s)
• pressure::Vector{T}: A vector of pressures inside each zone. Computed from the equation of state. (Unit: kg/(m⋅s²))
• intenergy::Vector{T}: A vector of the internal energy per unit mass of each zone. (Unit: m²/s²)
• speedofsound::Vector{T}: A vector of the speed of sound within each zone. Computed from the equation of state. (Unit: m/s)
• viscosity::Vector{T}: A vector of the artificial viscosity within each zone, added to the pressure field. See ArtificialDissipation.jl for more information. (Unit: kg/(m⋅s²))
• energy_flux::Vector{T}: A vector of fluxes of internal energy per unit mass across each zone boundary. Added to the energy equation as a diffusion term. See ArtificialDissipation.jl for more information. (Unit: m³/s³)
• momentum_rhs::Vector{T}: A vector of the right hand side of the momentum equation at each zone edge. (Unit: m/s²)
• energy_rhs::Vector{T}: A vector of the right hand side of the energy equation within each zone. (Unit: m²/s³)

struct SimulationCallback

A structure containing callbacks to be called by a simulation.

Parameters

• callback_cycle: A Vector of CycleCallback's that are called based on the number of cycles (timesteps) that have been performed
• callback_time: A Vector of TimeCallback's that are called at a fixed set of times
• callback_dt: A Vector of TimeDeltaCallback's that are called at a fixed temporal frequency

Notes

• This structure is typically initialized using ConfigureSimulationCallbacks(). See the documentation of that function for further detail.

struct CycleCallback

A structure containing information about a callback intended to be executed every every cycles

Parameters

• func: The Function to be executed
• every: How often (in cycles) to execute this callback
• last_called: The last cycle that this callback was called

Notes

• This type of callback is useful for cases where a callback should be called at a fixed number of cycles between each callback.
• While this structure can be initialized directly and added to the callback_cycle entry of a SimulationCallback structure, it is recommended to call RegisterCycleCallback!() instead.

struct TimeCallback

A structure containing information about a callback intended to be executed at a fixed list of times given by times.

Parameters

• func: The Function to be executed
• times: The list of times at which to execute this callback
• next_index: An index into the times vector indicating the next time the callback should be executed.

Notes

• This type of callback is useful if the callback should be called at an irregular series of times. See CycleCallback if the callback should be called at a regular number of cycles, or TimeDeltaCallback if the callback should be called at a fixed temporal cadence.
• While this structure can be initialized directly and added to the callback_time entry of a SimulationCallback structure, it is recommended to call RegisterTimeCallback!() instead.
• The entries in times are assumed to be sorted in ascending order. RegisterTimeCallback!() will handle this automatically, but this will need to be handled manually if creating this structure directly.

struct TimeDeltaCallback

A structure containing information about a callback intended to be executed every every seconds

Parameters

• func: The Function to be executed
• every: How often (in seconds) to execute this callback
• last_called: The last time that this callback was called

Notes

• This type of callback is useful for cases where a callback should be called at a fixed temporal spacing between each callback.
• While this structure can be initialized directly and added to the callback_dt entry of a SimulationCallback structure, it is recommended to call RegisterTimeDeltaCallback!() instead.

Momentum( m₋::T, P₋::T, m₊::T, P₊::T ) where { T <: AbstractFloat }

Compute the right hand side of the momentum equation at a zone interface. See the Notes section for information on the equation being solved.

Returns

A scalar of type T representing the rate of change of velocity of a zone interface over time.

Arguments

• m: The mass of the zone. (Unit: kg)
• P: The pressure of the zone (Unit: kg/(m⋅s²))

For each of these parameters, a - subscript refers to the zone to the left of the zone interface and a + subscript refers to the zone to the right of the zone interface.

Notes

The governing equation solved in this function is: ∂u/∂t = (1/ρ₀) * ∂P/∂x Through the numerical disretization, this reduces to ∂u/∂t = 1/m̄ ( P₊ - P₋ ) If using artificial viscosity per the method of Von Neumann and Richtmyer (1950), the artificial viscosity term should be added to the pressure field.

Energy( ρ::T, P::T, Δx::T, u₋::T, u₊::T, q₋::T, q₊::T ) where { T <: AbstractFloat }

Compute the right hand side of the energy equation in a zone. See the Notes section for the specific equations that are solved.

Returns

A scalar of type T representing the rate of change of internal energy inside the zone.

Arguments

• m: The mass of the zone. (Unit: kg)
• P: The pressure of the zone. (Unit: kg/(m⋅s²))
• Δx: The length of the zone. (Unit: m)
• u: The velocity of the zone edges on the (-): left and (+): right of the zone. (Unit: m/s)
• q: The (artificial) flux of internal energy across the (-): left and (+): right zone edges. (Unit: m³/s³)

Notes:

The governing equation solved in this function is ∂e/∂t = - ( P / ρ₀ ) * ∂u/∂x + ∑q Through discretization, this becomes ∂eᵢ/∂t = - ( Pᵢ / mᵢ ) * ( u₊ - u₋ ) + ∑q