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)
• end_time: The final time of the simulation. (Unit: s; Default: 1.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)
• min_timestep: The minimum allowable timestep size. Simulation will halt if timestep falls below this value. (Unit: s; Default: 1e-7)
• max_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.

AdvanceToTime( state::Simulation{T}, stoptime::T, Δt::T; exact::Bool=false ) 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)
• Δt: The timestep size. (Unit: s)
• exact: If true, try to stop as close as possible to stoptime by adjusting the final timestep size (Default: false)

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.

AdvanceToTime( state::Simulation{T}, stoptime::T; exact::Bool=false ) where { T <: AbstractFloat }

Advance the simulation with initial state state to a time specified by stoptime with a variable 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)
• exact: If true, try to stop as close as possible to stoptime by adjusting the final timestep size (Default: false)

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.
• The timestep size, Δt, is determined for each cycle based on the minimum time for an acoustic wave to traverse a zone. See the documentation for CalculateTimestepSize() for further details.
• 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.

AdvanceOneCycle( state::Simulation{T}, Δt::T ) where { T <: AbstractFloat }

Advance the simulation by one cycle with a timestep of Δt.

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.
• Δt: The size of the time step. (Unit: s)

Notes

• This function allocates a deepcopy() of the input state. The copy is modified and returned from this function.

AdvanceOneCycle( state::Simulation{T} ) 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

Notes

• This function allocates a deepcopy() of the input state and returns the copy.
• The timestep size, Δt, 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}, Δt::T ) where { T <: AbstractFloat }

Advance the simulation by one cycle with a timestep of Δt.

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.
• Δt: The size of the time step. (Unit: s)

Side Effects

• All fields of output are modified in-place.

AdvanceOneCycle!( output::Simulation{T}, input::Simulation{T} ) 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} representing the simulation state at the start of the cycle

Notes

• The timestep size, Δt, is determined based on the minimum time for an acoustic wave to traverse a zone. See CalculateTimestepSize() for further details.

Side Effects

• All fields of output are modified in-place.

AdvanceNCycles( state::Simulation{T}, ncycles::UInt, Δt::T ) where { T <: AbstractFloat }

Advance the simulation by ncycles cycles 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.
• ncycles: The number of cycles to advance.
• Δt: The size of the time step. (Unit: s)

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.

AdvanceNCycles( state::Simulation{T}, ncycles::UInt ) where { T <: AbstractFloat }

Advance the simulation by ncycles cycles with a variable 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.
• ncycles: The number of cycles to advance. (Unit: ⋅)

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.
• The timestep size, Δt, 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³)

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