plasma_plots.theory.exact

Exact solutions to verify simulations: the Riemann problem of gas dynamics (with vacuum), the dam break, diffusion, advection, the wave equation, and pressureless flow and its caustics.

Profiles take the coordinates first and the time second, f(x, t, ...), so that lambda x, t: f(x, t, ...) (or a field of the returned result) can be passed as the reference= of plasma-plots’ profile plots. Everything is plain numpy and vectorized: the coordinates and times broadcast against each other.

Where a state can be empty (a vacuum in gas dynamics, a dry bed in shallow water) its density or depth is zero and its velocity is nan (undefined, so plots leave a gap); momenta and energies there are zero.

Classes

Name	Description
LagrangianFlow	A pressureless (ballistic) flow followed along its fluid elements.
RiemannSolution	The primitive variables of a gas-dynamics solution.
ShallowWaterSolution	The depth and velocity of a shallow-water solution.
StarState	The star region of an Euler Riemann problem, between the left and the right wave.

Functions

Name	Description
advected	Compute a profile carried unchanged at constant velocity, u(x, t) = u0(x − v t).
caustic_time	Compute the time of the first caustic (shell crossing) of a pressureless flow.
dalembert	Compute d'Alembert's solution of the 1-D wave equation ∂²u/∂t² = c² ∂²u/∂x².
dam_break	Compute Ritter's exact dam-break solution of the 1-D shallow-water equations on a dry bed.
heat_kernel	Compute a spreading Gaussian, the solution of the diffusion equation ∂u/∂t = D ∇²u.
pressureless	Follow a 1-D pressureless flow (the Zel'dovich approximation) in Lagrangian coordinates.
pressureless_eulerian	Compute a 1-D pressureless flow as a function of the position, before shell crossing.
riemann_euler	Compute the exact solution of the Riemann problem of the 1-D Euler equations.
sod_shock_tube	Compute the exact solution of Sod's shock tube.
star_state	Solve for the star region of the Riemann problem of the 1-D Euler equations.

LagrangianFlowclassdataclass#

View source

class LagrangianFlow(position: object, velocity: object, density: object)

A pressureless (ballistic) flow followed along its fluid elements.

Attributes

Name	Type	Description
position	float or numpy.ndarray	The Eulerian position x = q + v0(q) t of each element.
velocity	float or numpy.ndarray	Its velocity v0(q), constant in time.
density	float or numpy.ndarray	The density at x, ρ0(q)/|1 + v0’(q) t|; inf at a caustic.

densityattributeinstance attribute#

View source

density: object

positionattributeinstance attribute#

View source

position: object

velocityattributeinstance attribute#

View source

velocity: object

RiemannSolutionclassdataclass#

View source

class RiemannSolution(density: object, velocity: object, pressure: object, internal_energy: object, gamma: float)

The primitive variables of a gas-dynamics solution.

Attributes

Name	Type	Description
density	float or numpy.ndarray	The mass density ρ.
velocity	float or numpy.ndarray	The velocity u; nan in a vacuum.
pressure	float or numpy.ndarray	The pressure p.
internal_energy	float or numpy.ndarray	The specific internal energy e = p/((γ − 1) ρ); nan in a vacuum.
gamma	float	The ratio of specific heats γ.

densityattributeinstance attribute#

View source

density: object

gammaattributeinstance attribute#

View source

gamma: float

internal_energyattributeinstance attribute#

View source

internal_energy: object

pressureattributeinstance attribute#

View source

pressure: object

velocityattributeinstance attribute#

View source

velocity: object

energyproperty#

View source

energy

The total energy density ρu²/2 + p/(γ − 1) (zero in a vacuum).

momentumproperty#

View source

momentum

The momentum density ρu (zero in a vacuum).

sound_speedproperty#

View source

sound_speed

The sound speed a = √(γ p/ρ); nan in a vacuum.

ShallowWaterSolutionclassdataclass#

View source

class ShallowWaterSolution(depth: object, velocity: object)

The depth and velocity of a shallow-water solution.

Attributes

Name	Type	Description
depth	float or numpy.ndarray	The water depth h (zero on a dry bed).
velocity	float or numpy.ndarray	The depth-averaged velocity u; nan on a dry bed.

depthattributeinstance attribute#

View source

depth: object

velocityattributeinstance attribute#

View source

velocity: object

dischargeproperty#

View source

discharge

The discharge (momentum per unit density) hu, zero on a dry bed.

StarStateclassdataclass#

View source

class StarState(pressure: float, velocity: float, density_left: float, density_right: float, left_wave: str, right_wave: str, vacuum: bool)

The star region of an Euler Riemann problem, between the left and the right wave.

Attributes

Name	Type	Description
pressure	float	The star pressure p* (zero if a vacuum forms).
velocity	float	The star (contact) velocity u*; nan if a vacuum forms or a side is a vacuum.
density_left	float	The density ρ*L left of the contact (zero for a vacuum).
density_right	float	The density ρ*R right of the contact (zero for a vacuum).
left_wave	str	"shock" or "rarefaction" ("none" if the left state is a vacuum).
right_wave	str	"shock" or "rarefaction" ("none" if the right state is a vacuum).
vacuum	bool	Whether the solution contains a vacuum (given, or generated between two rarefactions).

density_leftattributeinstance attribute#

View source

density_left: float

density_rightattributeinstance attribute#

View source

density_right: float

left_waveattributeinstance attribute#

View source

left_wave: str

pressureattributeinstance attribute#

View source

pressure: float

right_waveattributeinstance attribute#

View source

right_wave: str

vacuumattributeinstance attribute#

View source

vacuum: bool

velocityattributeinstance attribute#

View source

velocity: float

advectedfunction#

View source

def advected(profile, x, t, velocity, period=None)

Compute a profile carried unchanged at constant velocity, u(x, t) = u0(x − v t).

Parameters

Name	Type	Default	Description
profile	callable	required	The initial profile u0, a vectorized function of the position.
x	float or array_like	required	Positions.
t	float or array_like	required	Times, broadcast against x.
velocity	float	required	The advection velocity v.
period	float or (float, float)	None	For a periodic domain: its length L (the domain is [0, L)) or its bounds (a, b). The shifted position is wrapped into the domain before profile is called. Default: None (not periodic).

Returns

float or numpy.ndarray

u0(x − v t).

Examples

>>> # −0.25 wrapped into [0, 1)
>>> round(float(advected(lambda x: x, 0.1, 0.35, 1.0, period=1.0)), 6)
0.75
>>> advected(lambda x: np.exp(-(x**2)), [0.0, 1.0], 1.0, 1.0).round(4)
array([0.3679, 1.    ])

caustic_timefunction#

View source

def caustic_time(velocity, q)

Compute the time of the first caustic (shell crossing) of a pressureless flow.

Elements with v0’(q) < 0 meet at t = −1/v0’(q), so the first caustic forms at t_c = −1/min v0’(q), at the steepest decrease of the initial velocity.

Parameters

Name	Type	Description
velocity	callable	The initial velocity v0(q), a vectorized function.
q	array_like	The positions where v0’ is sampled (the minimum is taken over these, so resolve it).

Returns

float

t_c; inf if the velocity nowhere decreases.

See Also

pressureless() : The flow up to (and past) the caustic.

References

Ya. B. Zel’dovich, “Gravitational instability: an approximate theory for large density perturbations”, Astron. Astrophys. 5, 84 (1970).

Examples

>>> q = np.linspace(-np.pi, np.pi, 1001)
>>> round(caustic_time(lambda q: -0.5 * np.sin(q), q), 6)
2.0

dalembertfunction#

View source

def dalembert(x, t, initial, speed, initial_rate=None)

Compute d’Alembert’s solution of the 1-D wave equation ∂²u/∂t² = c² ∂²u/∂x².

u = [u0(x − ct) + u0(x + ct)]/2 + (1/(2c)) ∫_{x−ct}^{x+ct} v0(s) ds, for the initial profile u0 and the initial rate v0 = ∂u/∂t at t = 0, on the whole line. The integral is done by 64-point Gauss–Legendre quadrature, exact to rounding for smooth v0 over a few wavelengths. A pulse at rest splits into two halves moving at ±c.

Parameters

Name	Type	Default	Description
x	float or array_like	required	Positions.
t	float or array_like	required	Times, broadcast against x.
initial	callable	required	The initial profile u0, a vectorized function of the position.
speed	float	required	The wave speed c > 0.
initial_rate	callable	None	The initial rate v0, a vectorized function of the position. Default: None (at rest).

Returns

float or numpy.ndarray

u(x, t).

Raises

ValueError

If speed is not positive.

References

J. le Rond d’Alembert, “Recherches sur la courbe que forme une corde tenduë mise en vibration”, Hist. Acad. R. Sci. Berlin 3, 214 (1747); L. C. Evans, Partial Differential Equations, 2nd ed. (AMS, 2010), section 2.4.

Examples

>>> pulse = lambda x: np.exp(-(x**2) / 0.01)
>>> # the pulse split in two halves
>>> dalembert([0.0, 1.0], 1.0, pulse, 1.0).round(4)
array([0. , 0.5])
>>> # sin(ct)/c
>>> round(
...     float(dalembert(0.0, 0.5, np.zeros_like, 1.0, initial_rate=np.cos)), 6
... )
0.479426

dam_breakfunction#

View source

def dam_break(x, t, depth=1.0, gravity=1.0, x0=0.0)

Compute Ritter’s exact dam-break solution of the 1-D shallow-water equations on a dry bed.

Water at rest of depth h0 for x < x0 is released at t = 0 onto a dry bed. A rarefaction runs between x0 − c0 t and the front at x0 + 2 c0 t, where c0 = √(g h0), with h = (2c0 − ξ)²/(9g) and u = 2(c0 + ξ)/3 for ξ = (x − x0)/t.

Parameters

Name	Type	Default	Description
x	float or array_like	required	Positions.
t	float or array_like	required	Times ≥ 0, broadcast against x.
depth	float	1.0	The initial depth h0 behind the dam (x < x0). Default: 1.0.
gravity	float	1.0	The gravitational acceleration g. Default: 1.0.
x0	float	0.0	The position of the dam. Default: 0.0.

Returns

ShallowWaterSolution

Depth and velocity (nan on the dry bed), with the broadcast shape of x and t; also the discharge hu as a property.

Raises

ValueError

If a time is negative, or the depth not positive.

References

A. Ritter, “Die Fortpflanzung der Wasserwellen”, Z. Ver. Dtsch. Ing. 36, 947 (1892).

Examples

>>> w = dam_break([-2.0, 0.0, 1.0, 2.1], 1.0)
>>> w.depth.round(6), w.velocity.round(6)
(array([1.      , 0.444444, 0.111111, 0.      ]), array([0.      , 0.666667, 1.333333,      nan]))

heat_kernelfunction#

View source

def heat_kernel(x, t, diffusivity, width=0.0, center=0.0, mass=1.0)

Compute a spreading Gaussian, the solution of the diffusion equation ∂u/∂t = D ∇²u.

u = M (2π σ²)^(−d/2) exp(−|x − x_c|²/(2σ²)) with σ² = σ0² + 2 D t: the fundamental solution for σ0 = 0, and a Gaussian of initial standard deviation σ0 otherwise. The integral is M at all times; each coordinate’s variance grows as σ0² + 2Dt.

Parameters

Name	Type	Default	Description
x	float, array_like or tuple of array_like	required	Positions: an array in 1-D, or a tuple (x, y) or (x, y, z) of arrays that broadcast against each other in d dimensions.
t	float or array_like	required	Times ≥ 0, broadcast against the positions.
diffusivity	float	required	The diffusion coefficient D.
width	float	0.0	The initial standard deviation σ0. Default: 0.0 (a point source at t = 0).
center	float or tuple of float	0.0	The center x_c, one value per dimension (a scalar is used for all). Default: 0.0.
mass	float	1.0	The integral M of u. Default: 1.0.

Returns

float or numpy.ndarray

u; at t = 0 with σ0 = 0, inf at the center and 0 elsewhere.

Raises

ValueError

If a time or diffusivity or width is negative.

References

L. C. Evans, Partial Differential Equations, 2nd ed. (AMS, 2010), section 2.3.

Examples

>>> round(float(heat_kernel(0.0, 0.5, 1.0)), 6)  # 1/√(4π D t)
0.398942
>>> round(float(heat_kernel((1.0, 0.0), 0.25, 1.0, width=1.0)), 6)
0.076026

pressurelessfunction#

View source

def pressureless(q, t, velocity, density=None, velocity_derivative=None)

Follow a 1-D pressureless flow (the Zel’dovich approximation) in Lagrangian coordinates.

Without pressure every fluid element keeps its initial velocity, x(q, t) = q + v0(q) t, and mass conservation ρ dx = ρ0 dq gives ρ = ρ0(q)/|∂x/∂q| = ρ0(q)/|1 + v0’(q) t|. Where the velocity decreases, elements catch up with each other: the density becomes infinite (a caustic) at caustic_time(), after which the map is multivalued (shell crossing) and the density here is that of each stream.

Parameters

Name	Type	Default	Description
q	float or array_like	required	The initial (Lagrangian) positions of fluid elements.
t	float or array_like	required	Times, broadcast against q.
velocity	callable	required	The initial velocity v0(q), a vectorized function.
density	callable	None	The initial density ρ0(q), a vectorized function. Default: None (uniform, 1).
velocity_derivative	callable	None	v0’(q). Default: None (by finite differences).

Returns

LagrangianFlow

The positions, velocities and densities of the elements.

See Also

pressureless_eulerian() : The same flow as a function of the position x. caustic_time() : The time of the first shell crossing.

References

Ya. B. Zel’dovich, “Gravitational instability: an approximate theory for large density perturbations”, Astron. Astrophys. 5, 84 (1970).

Examples

A converging sine velocity, v0 = −0.5 sin(q), with its caustic at t = 2:

>>> flow = pressureless([0.0, np.pi / 2], 1.0, lambda q: -0.5 * np.sin(q))
>>> flow.position.round(6), flow.density.round(6)
(array([0.      , 1.070796]), array([2., 1.]))

pressureless_eulerianfunction#

View source

def pressureless_eulerian(x, t, velocity, density=None, velocity_derivative=None)

Compute a 1-D pressureless flow as a function of the position, before shell crossing.

Inverts x = q + v0(q) t for the Lagrangian position q (by Newton’s method safeguarded by bisection, to machine precision) and returns the flow there. The inverse is unique only before the first caustic (caustic_time()); later, one of the streams is returned.

Parameters

Name	Type	Default	Description
x	float or array_like	required	Positions.
t	float or array_like	required	Times, broadcast against x.
velocity	callable	required	The initial velocity v0(q), a vectorized function.
density	callable	None	The initial density ρ0(q), a vectorized function. Default: None (uniform, 1).
velocity_derivative	callable	None	v0’(q). Default: None (by finite differences).

Returns

LagrangianFlow

position is x itself; the velocity and density at x.

See Also

pressureless() : The flow along the fluid elements.

References

Ya. B. Zel’dovich, “Gravitational instability: an approximate theory for large density perturbations”, Astron. Astrophys. 5, 84 (1970).

Examples

>>> flow = pressureless_eulerian(
...     [-1.070796, 1.070796], 1.0, lambda q: -0.5 * np.sin(q)
... )
>>> flow.density.round(5), flow.velocity.round(5)
(array([1., 1.]), array([ 0.5, -0.5]))

riemann_eulerfunction#

View source

def riemann_euler(x, t, left, right, gamma=1.4, x0=0.0)

Compute the exact solution of the Riemann problem of the 1-D Euler equations.

The gas is ideal, p = (γ − 1) ρ e, and initially in the state left for x < x0 and right for x > x0. The solution is self-similar in (x − x0)/t: a left wave (shock or rarefaction), a contact discontinuity moving at u* and a right wave. All configurations are covered, including a vacuum on either side (gas expanding into vacuum, bounded by a front moving at u ± 2a/(γ − 1)) and the vacuum generated between two strong rarefactions.

Parameters

Name	Type	Default	Description
x	float or array_like	required	Positions.
t	float or array_like	required	Times ≥ 0, broadcast against x. At t = 0 the initial states are returned.
left	(float, float, float)	required	The left state (density, velocity, pressure); (0, u, 0) is a vacuum.
right	(float, float, float)	required	The right state (density, velocity, pressure); (0, u, 0) is a vacuum.
gamma	float	1.4	The ratio of specific heats γ > 1. Default: 1.4.
x0	float	0.0	The position of the initial discontinuity. Default: 0.0.

Returns

RiemannSolution

Density, velocity, pressure and specific internal energy, with the broadcast shape of x and t; also the momentum and energy densities as properties.

Raises

ValueError

If gamma ≤ 1, a time is negative, or a state is invalid (see star_state()).

See Also

star_state() : The star pressure and velocity alone. sod_shock_tube() : The classic test case.

References

E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed. (Springer, 2009), chapter 4 (4.9 for vacuum).

Examples

Toro’s test 1 at t = 0.25, left of, inside and right of the star region:

>>> w = riemann_euler(
...     [-0.4, 0.1, 0.5], 0.25, (1.0, 0.0, 1.0), (0.125, 0.0, 0.1)
... )
>>> w.density.round(5), w.pressure.round(5)
(array([1.     , 0.42632, 0.125  ]), array([1.     , 0.30313, 0.1    ]))

Gas expanding into vacuum: the front moves at 2a/(γ − 1) = 5.9161 here.

>>> w = riemann_euler([0.0, 3.0, 6.0], 1.0, (1.0, 0.0, 1.0), (0.0, 0.0, 0.0))
>>> w.density.round(5), w.velocity.round(5)
(array([0.40188, 0.01169, 0.     ]), array([0.98601, 3.48601,     nan]))

sod_shock_tubefunction#

View source

def sod_shock_tube(x, t, gamma=1.4, x0=0.5)

Compute the exact solution of Sod’s shock tube.

The Riemann problem with (ρ, u, p) = (1, 0, 1) left and (0.125, 0, 0.1) right of x0, usually on [0, 1] up to t = 0.2 (Toro’s test 1): a left rarefaction, a contact and a right shock.

Parameters

Name	Type	Default	Description
x	float or array_like	required	Positions.
t	float or array_like	required	Times ≥ 0, broadcast against x.
gamma	float	1.4	The ratio of specific heats. Default: 1.4.
x0	float	0.5	The position of the diaphragm. Default: 0.5.

Returns

RiemannSolution

Density, velocity, pressure and specific internal energy.

See Also

riemann_euler() : Any Riemann problem.

References

G. A. Sod, “A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws”, J. Comput. Phys. 27, 1 (1978).

Examples

>>> w = sod_shock_tube([0.1, 0.6, 0.8, 0.9], 0.2)
>>> w.density.round(5)
array([1.     , 0.42632, 0.26557, 0.125  ])
>>> float(w.velocity[2].round(5)), float(w.pressure[2].round(5))
(0.92745, 0.30313)

star_statefunction#

View source

def star_state(left, right, gamma=1.4)

Solve for the star region of the Riemann problem of the 1-D Euler equations.

The star pressure p* is the root of Toro’s pressure function f_L(p) + f_R(p) + u_R − u_L, found by Newton’s method safeguarded by bisection to machine precision. A wave is a shock if p* exceeds the pressure ahead of it and a rarefaction otherwise. If a side is a vacuum (density and pressure zero), or if 2(a_L + a_R)/(γ − 1) ≤ u_R − u_L so that the two rarefactions leave a vacuum between them, p* = 0.

Parameters

Name	Type	Default	Description
left	(float, float, float)	required	The left state (density, velocity, pressure).
right	(float, float, float)	required	The right state (density, velocity, pressure).
gamma	float	1.4	The ratio of specific heats γ > 1. Default: 1.4.

Returns

StarState

p*, u*, the densities on both sides of the contact and the types of the two waves.

Raises

ValueError

If gamma ≤ 1, or a state has a negative density or pressure, or only one of the two zero.

See Also

riemann_euler() : The whole solution.

References

E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed. (Springer, 2009), chapter 4.

Examples

Toro’s test 1 (Sod’s shock tube):

>>> s = star_state((1.0, 0.0, 1.0), (0.125, 0.0, 0.1))
>>> (
...     round(s.pressure, 5),
...     round(s.velocity, 5),
...     round(s.density_left, 5),
...     round(s.density_right, 5),
... )
(0.30313, 0.92745, 0.42632, 0.26557)
>>> s.left_wave, s.right_wave
('rarefaction', 'shock')