Fluids

Fluid dynamics: flows, convection cells, turbulence, shocks, shallow water, wave trains, and kinetic gases.

BhatnagarGrossKrookInteraction(label, *, target, state='f', axes=('v',), time=1.0)

Bases: RelaxationInteraction

Relaxation of Bhatnagar, Gross, and Krook of a density towards the maximum-entropy density with its moments.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the field on a phase space.	required
state	str	The name of the state of the density.	'f'
axes	Sequence[str]	The names of the velocity axes.	('v',)
time	Callable[[SimpleNamespace], object] | float	A function of the namespace that returns the relaxation time, or a positive number.	1.0

Examples:

>>> from skcomplex.models.physics.fluids import BhatnagarGrossKrookInteraction
>>> BhatnagarGrossKrookInteraction('bgk', target='gas', axes=('v',), time=0.1).check().axes
('v',)

Source code in src/skcomplex/models/physics/fluids/_bhatnagar_gross_krook.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    state: str = 'f',
    axes: Sequence[str] = ('v',),
    time: Callable[[SimpleNamespace], object] | float = 1.0,
) -> None:
    super().__init__(label, target=target, state=state, goal=0.0, time=time)
    self.axes = axes

BoussinesqElement(label, visual=None, *, omega=0.0, theta=0.0, Ra=1000.0, Pr=1.0)

Bases: Element

Vorticity and temperature of Boussinesq convection, in units of the thermal diffusion time.

The temperature theta is the departure from the linear profile between a hot bottom and a cold top, so dtheta/dt + u . grad(theta) = w + laplacian(theta) and domega/dt + u . grad(omega) = Pr laplacian(omega) + Ra Pr dtheta/dx. The advections, the diffusions, and the buoyancy are interactions of the model.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
omega	Annotated[float, StateSpec()]	The initial vorticity.	0.0
theta	Annotated[float, StateSpec()]	The initial temperature departure.	0.0
Ra	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The Rayleigh number.	1000.0
Pr	Annotated[float, ParameterSpec(domain=open(0, oo))]	The Prandtl number.	1.0

Examples:

>>> from skcomplex.models.physics.fluids import BoussinesqElement
>>> BoussinesqElement('site').Ra
1000.0

Source code in src/skcomplex/models/physics/fluids/_boussinesq.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    omega: Annotated[float, StateSpec()] = 0.0,
    theta: Annotated[float, StateSpec()] = 0.0,
    Ra: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1000.0,
    Pr: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.omega = omega
    self.theta = theta
    self.Ra = Ra
    self.Pr = Pr

BurgersElement(label, visual=None, *, u=0.0)

Bases: Element

Site of a field u that the Burgers equation du/dt + d(u^2 / 2)/dx = 0 carries.

Parameters:

Name	Type	Description	Default
label	str	The label of the site.	required
visual	VisualSpec | None	How the site is drawn, or None for none.	None
u	Annotated[float, StateSpec()]	The initial value.	0.0

Examples:

>>> from skcomplex.models.physics.fluids import BurgersElement
>>> BurgersElement('site', u=1.5).u
1.5

Source code in src/skcomplex/models/physics/fluids/_burgers.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    u: Annotated[float, StateSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.u = u

CarlemanGasElement(label, visual=None, *, f_plus=1.0, f_minus=1.0, c=1.0, kappa=1.0)

Bases: Element

Site of the Carleman discrete-velocity gas, densities f_plus and f_minus moving at the speeds c and -c.

Collisions turn two particles of one velocity into two of the other at the rate kappa, so df_plus/dt + c df_plus/dx = kappa (f_minus^2 - f_plus^2) and the opposite for f_minus, which conserves the mass and relaxes to the equilibrium f_plus = f_minus.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
f_plus	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial density of right-moving particles.	1.0
f_minus	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial density of left-moving particles.	1.0
c	Annotated[float, ParameterSpec(domain=open(0, oo))]	The particle speed.	1.0
kappa	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The collision rate.	1.0

Examples:

>>> from skcomplex.models.physics.fluids import CarlemanGasElement
>>> CarlemanGasElement('site').kappa
1.0

Source code in src/skcomplex/models/physics/fluids/_carleman_gas.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    f_plus: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    f_minus: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    c: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    kappa: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.f_plus = f_plus
    self.f_minus = f_minus
    self.c = c
    self.kappa = kappa

FHPSiteElement(label, visual=None, *, n0=0, n1=0, n2=0, n3=0, n4=0, n5=0, density=0)

Bases: Element

Site of the FHP-I lattice gas, with one boolean channel per direction of the triangular lattice.

The local rules of the gas collide the particles at every site, head-on pairs turning by a sixth of a turn to either side with equal probability and symmetric triples turning by a sixth of a turn, and then move every particle to the neighbouring site along its direction. Collisions conserve the number of particles and their momentum at every site. The state density holds the number of particles at the site.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
n0	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The initial occupation of the channel at 0 degrees.	0
n1	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The initial occupation of the channel at 60 degrees.	0
n2	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The initial occupation of the channel at 120 degrees.	0
n3	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The initial occupation of the channel at 180 degrees.	0
n4	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The initial occupation of the channel at 240 degrees.	0
n5	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The initial occupation of the channel at 300 degrees.	0
density	Annotated[float, StateSpec(domain=NN)]	The initial number of particles at the site.	0

Examples:

>>> from skcomplex.models.physics.fluids import FHPSiteElement
>>> FHPSiteElement('site').n0
0

Source code in src/skcomplex/models/physics/fluids/_lattice_gases.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    n0: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
    n1: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
    n2: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
    n3: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
    n4: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
    n5: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
    density: Annotated[float, StateSpec(domain=NN)] = 0,
) -> None:
    super().__init__(label, visual=visual)
    self.n0 = n0
    self.n1 = n1
    self.n2 = n2
    self.n3 = n3
    self.n4 = n4
    self.n5 = n5
    self.density = density

KuramotoSivashinskyElement(label, visual=None, *, u=0.0)

Bases: Element

Site of a flame front or a thin film, du/dt = -d(u^2 / 2)/dx - laplacian(u) - laplacian(laplacian(u)).

The conservation of the flux u^2 / 2 and the rate of the linear operator are interactions of the model.

Parameters:

Name	Type	Description	Default
label	str	The label of the site.	required
visual	VisualSpec | None	How the site is drawn, or None for none.	None
u	Annotated[float, StateSpec()]	The initial height of the front.	0.0

Examples:

>>> from skcomplex.models.physics.fluids import KuramotoSivashinskyElement
>>> KuramotoSivashinskyElement('site', u=0.5).u
0.5

Source code in src/skcomplex/models/physics/fluids/_kuramoto_sivashinsky.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    u: Annotated[float, StateSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.u = u

MarangoniLayerElement(label, visual=None, *, u=0.0, r=0.1, g=1.0)

Bases: Element

Site of a convecting layer with a free surface, with the order parameter u.

Parameters:

Name	Type	Description	Default
label	str	The label of the site.	required
visual	VisualSpec | None	How the site is drawn, or None for none.	None
u	Annotated[float, StateSpec()]	The initial order parameter.	0.0
r	Annotated[float, ParameterSpec()]	The distance from the onset of convection.	0.1
g	Annotated[float, ParameterSpec()]	The strength of the quadratic term that the free surface adds.	1.0

Examples:

>>> from skcomplex.models.physics.fluids import MarangoniLayerElement
>>> MarangoniLayerElement('site', g=0.5).g
0.5

Source code in src/skcomplex/models/physics/fluids/_marangoni.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    u: Annotated[float, StateSpec()] = 0.0,
    r: Annotated[float, ParameterSpec()] = 0.1,
    g: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.u = u
    self.r = r
    self.g = g

ShallowWaterElement(label, visual=None, *, h=1.0, q=0.0, g=9.81)

Bases: Element

Site of a channel of shallow water, with the depth h and the discharge q per unit width.

The mass h and the momentum q are conserved by two conservation interactions, with the fluxes q and q^2 / h + g h^2 / 2, which form one hyperbolic system.

Parameters:

Name	Type	Description	Default
label	str	The label of the site.	required
visual	VisualSpec | None	How the site is drawn, or None for none.	None
h	Annotated[float, StateSpec(domain=open(0, oo))]	The initial depth.	1.0
q	Annotated[float, StateSpec()]	The initial discharge.	0.0
g	Annotated[float, ParameterSpec()]	The gravitational acceleration.	9.81

Examples:

>>> from skcomplex.models.physics.fluids import ShallowWaterElement
>>> ShallowWaterElement('site', h=2.0).h
2.0

Source code in src/skcomplex/models/physics/fluids/_shallow_water.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    h: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 1.0,
    q: Annotated[float, StateSpec()] = 0.0,
    g: Annotated[float, ParameterSpec()] = 9.81,
) -> None:
    super().__init__(label, visual=visual)
    self.h = h
    self.q = q
    self.g = g

VorticityElement(label, visual=None, *, omega=0.0, nu=0.001, f=0.0)

Bases: Element

Vorticity of an incompressible two-dimensional flow, domega/dt + u . grad(omega) = nu laplacian(omega) + f.

The velocity comes from the streamfunction, the solution of a Poisson equation, so the flow is divergence free, and the advection is written in conservative form so that dealiasing applies to it. The advection, the viscosity, and the forcing are interactions of the model.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
omega	Annotated[float, StateSpec()]	The initial vorticity.	0.0
nu	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The kinematic viscosity.	0.001
f	Annotated[float, ParameterSpec()]	The forcing of the vorticity.	0.0

Examples:

>>> from skcomplex.models.physics.fluids import VorticityElement
>>> VorticityElement('site').nu
0.001

Source code in src/skcomplex/models/physics/fluids/_vorticity.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    omega: Annotated[float, StateSpec()] = 0.0,
    nu: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.001,
    f: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.omega = omega
    self.nu = nu
    self.f = f

WaveNumberElement(label, visual=None, *, k=1.0, g=1.0)

Bases: Element

Local wave number of deep-water waves, conserved as dk/dt + d omega(k)/dx = 0 with omega = sqrt(g k).

The wave number travels at the group speed C(k) = sqrt(g / k) / 2, which falls as the waves shorten, so long waves catch up with short ones and the wave number steepens until it breaks.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
k	Annotated[float, StateSpec(domain=open(0, oo))]	The initial wave number.	1.0
g	Annotated[float, ParameterSpec(domain=open(0, oo))]	The gravitational acceleration.	1.0

Examples:

>>> from skcomplex.models.physics.fluids import WaveNumberElement
>>> WaveNumberElement('site').g
1.0

Source code in src/skcomplex/models/physics/fluids/_wave_number.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    k: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 1.0,
    g: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.k = k
    self.g = g

benard_marangoni(shape=(64, 64), *, length=16 * math.pi, r=0.1, g=1.0, random_state=0)

Return the order-parameter equation of Bénard-Marangoni convection, which forms hexagonal cells.

The field u is the vertical velocity at mid-depth of a heated layer whose open surface has a surface tension that falls with temperature. It follows u_t = r u - (1 + laplacian)^2 u + g u^2 - u^3, the Swift-Hohenberg equation of convection with the quadratic term that the thermocapillary surface adds by breaking the symmetry between rising and sinking fluid. With g > 0 and small r > 0, the cells are hexagons of rising fluid with a wavenumber near 1, and with g = 0, the layer of a fixed surface, they are rolls.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of sites along every axis of the periodic lattice.	(64, 64)
length	float	The side of the domain.	16 * pi
r	float	The distance above onset.	0.1
g	float	The strength of the quadratic term of the thermocapillary surface.	1.0
random_state	int | Generator | None	The seed or generator of the small random initial field.	0

Returns:

Name	Type	Description
system	System	The model, with the field layer.u.

Examples:

>>> from skcomplex.models.physics.fluids import benard_marangoni
>>> benard_marangoni((8, 8), g=0.5).elements[0].values['g']
0.5

Source code in src/skcomplex/models/physics/fluids/_marangoni.py

@register_model(
    field='physics/fluids',
    tags=('field', 'convection', 'hexagons', 'surface tension', 'spectral', 'pde'),
    references=('Satz 2022', 'Bestehorn 1993', 'Cross and Hohenberg 1993'),
    defaults={'until': 300.0, 'dt': 5.0, 'solver': ETDRK4Solver(dt=0.5)},
)
def benard_marangoni(
    shape: Sequence[int] = (64, 64),
    *,
    length: float = 16 * math.pi,
    r: float = 0.1,
    g: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the order-parameter equation of Bénard-Marangoni convection, which forms hexagonal cells.

    The field `u` is the vertical velocity at mid-depth of a heated layer whose open surface has a surface tension
    that falls with temperature. It follows `u_t = r u - (1 + laplacian)^2 u + g u^2 - u^3`, the Swift-Hohenberg
    equation of convection with the quadratic term that the thermocapillary surface adds by breaking the symmetry
    between rising and sinking fluid. With `g > 0` and small `r > 0`, the cells are hexagons of rising fluid with a
    wavenumber near 1, and with `g = 0`, the layer of a fixed surface, they are rolls.

    Args:
        shape:
            The number of sites along every axis of the periodic lattice.

        length:
            The side of the domain.

        r:
            The distance above onset.

        g:
            The strength of the quadratic term of the thermocapillary surface.

        random_state:
            The seed or generator of the small random initial field.

    Returns:
        system:
            The model, with the field `layer.u`.

    Examples:
        >>> from skcomplex.models.physics.fluids import benard_marangoni
        >>> benard_marangoni((8, 8), g=0.5).elements[0].values['g']
        0.5
    """
    shape = tuple(shape)
    rng = np.random.default_rng(random_state)
    values = {'u': 0.01 * rng.uniform(-1.0, 1.0, shape), 'r': r, 'g': g}
    space = LatticeSpace(shape, spacing=length / shape[0])
    return System(
        elements=[
            FieldElement('layer', element=MarangoniLayerElement, space=space, values=values, discretization='spectral'),
        ],
        interactions=[
            GradientInteraction(
                'relaxation',
                target='layer',
                states='u',
                potential=_get_marangoni_density,
                variational=True,
            ),
        ],
    )

bgk_equation(positions=100, velocities=32, *, knudsen=0.01, half_width=6.0)

Return the BGK equation of the Sod shock tube over one position and one velocity axis.

The density relaxes at every position towards the maximum-entropy density with its density, momentum, and energy, and the gas starts from the two states of the shock tube of Sod on either side of the middle of a unit interval. As the Knudsen number falls the density moment tends to the Euler solution with the adiabatic index 3 of one velocity dimension.

Parameters:

Name	Type	Description	Default
positions	int	The number of positions.	100
velocities	int	The number of points of the velocity axis.	32
knudsen	float	The Knudsen number, the relaxation time.	0.01
half_width	float	The half-width of the velocity axis.	6.0

Returns:

Name	Type	Description
system	System	The model, with the density gas.f.

Examples:

>>> from skcomplex.models.physics.fluids import bgk_equation
>>> bgk_equation(8, 16).elements[0].check().prepare().size_
128

Source code in src/skcomplex/models/physics/fluids/_kinetic_equations.py

@register_model(
    field='physics/fluids',
    tags=('field', 'kinetic', 'phase space', 'shocks', 'pde'),
    references=('Bhatnagar, Gross, and Krook 1954', 'Mieussens 2000', 'Sod 1978'),
    defaults={'until': 0.2, 'dt': 0.05, 'solver': StrangSplittingSolver(SSPRKSolver(), RK4Solver(0.005), 0.005)},
)
def bgk_equation(
    positions: int = 100,
    velocities: int = 32,
    *,
    knudsen: float = 1e-2,
    half_width: float = 6.0,
) -> System:
    """Return the BGK equation of the Sod shock tube over one position and one velocity axis.

    The density relaxes at every position towards the maximum-entropy density with its density, momentum, and energy,
    and the gas starts from the two states of the shock tube of Sod on either side of the middle of a unit interval. As
    the Knudsen number falls the density moment tends to the Euler solution with the adiabatic index 3 of one velocity
    dimension.

    Args:
        positions:
            The number of positions.

        velocities:
            The number of points of the velocity axis.

        knudsen:
            The Knudsen number, the relaxation time.

        half_width:
            The half-width of the velocity axis.

    Returns:
        system:
            The model, with the density `gas.f`.

    Examples:
        >>> from skcomplex.models.physics.fluids import bgk_equation
        >>> bgk_equation(8, 16).elements[0].check().prepare().size_
        128
    """
    position = LatticeSpace((positions,), spacing=1.0 / positions, boundary='reflecting')
    space = PhaseSpace(position, [AxisSpec('v', velocities, -half_width, half_width)]).check().prepare()
    left = np.arange(positions) < positions // 2
    density = np.where(left, 1.0, 0.125)
    temperature = np.where(left, 1.0, 0.8)
    values = _get_maxwellian(density, np.zeros(positions), temperature, space.coordinates_['v'])
    gas = FieldElement(
        'gas',
        element=DensityElement,
        space=space,
        values={'f': values},
        discretization='finite-volume',
        visual=VLASOV_VISUAL,
    )
    stream = AdvectionInteraction('stream', target='gas', state='f', drift=_get_streaming)
    relax = BhatnagarGrossKrookInteraction('relax', target='gas', axes=('v',), time=knudsen)
    return System(elements=[gas], interactions=[stream, relax])

boltzmann_equation(positions=8, velocities=12, *, knudsen=1.0, half_width=6.0)

Return the Boltzmann equation of Maxwell molecules over one position and two velocity axes.

Free streaming moves a density wave along the periodic position axis, and binary collisions redistribute the density over the velocities at every position, by the fast spectral method.

Parameters:

Name	Type	Description	Default
positions	int	The number of positions.	8
velocities	int	The number of points of each velocity axis.	12
knudsen	float	The Knudsen number, the inverse of the collision frequency.	1.0
half_width	float	The half-width of each velocity axis.	6.0

Returns:

Name	Type	Description
system	System	The model, with the density gas.f.

Examples:

>>> from skcomplex.models.physics.fluids import boltzmann_equation
>>> boltzmann_equation(4, 8).elements[0].check().prepare().size_
256

Source code in src/skcomplex/models/physics/fluids/_kinetic_equations.py

@register_model(
    field='physics/fluids',
    tags=('field', 'kinetic', 'phase space', 'collisions', 'pde'),
    references=('Bellomo et al. 2017', 'Mouhot and Pareschi 2006'),
    defaults={
        'until': 0.5,
        'dt': 0.1,
        'solver': StrangSplittingSolver(SemiLagrangianSolver(0.1), RK4Solver(0.05), 0.1),
    },
)
def boltzmann_equation(
    positions: int = 8,
    velocities: int = 12,
    *,
    knudsen: float = 1.0,
    half_width: float = 6.0,
) -> System:
    """Return the Boltzmann equation of Maxwell molecules over one position and two velocity axes.

    Free streaming moves a density wave along the periodic position axis, and binary collisions redistribute the density
    over the velocities at every position, by the fast spectral method.

    Args:
        positions:
            The number of positions.

        velocities:
            The number of points of each velocity axis.

        knudsen:
            The Knudsen number, the inverse of the collision frequency.

        half_width:
            The half-width of each velocity axis.

    Returns:
        system:
            The model, with the density `gas.f`.

    Examples:
        >>> from skcomplex.models.physics.fluids import boltzmann_equation
        >>> boltzmann_equation(4, 8).elements[0].check().prepare().size_
        256
    """
    gas, frequency = _get_uniform_gas(positions, velocities, knudsen, half_width)
    stream = AdvectionInteraction('stream', target='gas', state='f', drift=_get_plane_streaming)
    collide = EncounterInteraction(
        'collide',
        candidate='gas',
        cross_section=_get_maxwell_molecules,
        factor=lambda _: frequency,
    )
    return System(elements=[gas], interactions=[stream, collide])

burgers(n=200, *, mean=1.5, amplitude=1.0)

Return the inviscid Burgers equation on a periodic interval of length 2 pi, whose sine steepens into a shock.

The cells hold averages updated by MUSCL reconstruction with a minmod limiter and Rusanov fluxes, so the total is conserved and the shock forms no new extrema. The shock forms at time 1 / amplitude.

Parameters:

Name	Type	Description	Default
n	int	The number of cells.	200
mean	float	The mean of the initial wave.	1.5
amplitude	float	The amplitude of the initial wave.	1.0

Returns:

Name	Type	Description
system	System	The model, with the velocity line.u.

Examples:

>>> from skcomplex.models.physics.fluids import burgers
>>> burgers(16).elements[0].check().prepare().size_
16

Source code in src/skcomplex/models/physics/fluids/_flows.py

@register_model(
    field='physics/fluids',
    tags=('field', 'shocks', 'conservation law', 'finite volume', 'pde'),
    references=('Burgers 1948', 'LeVeque 2002'),
    defaults={'until': 2.0, 'dt': 0.1, 'solver': SSPRKSolver()},
)
def burgers(n: int = 200, *, mean: float = 1.5, amplitude: float = 1.0) -> System:
    """Return the inviscid Burgers equation on a periodic interval of length 2 pi, whose sine steepens into a shock.

    The cells hold averages updated by MUSCL reconstruction with a minmod limiter and Rusanov fluxes, so the total is
    conserved and the shock forms no new extrema. The shock forms at time `1 / amplitude`.

    Args:
        n:
            The number of cells.

        mean:
            The mean of the initial wave.

        amplitude:
            The amplitude of the initial wave.

    Returns:
        system:
            The model, with the velocity `line.u`.

    Examples:
        >>> from skcomplex.models.physics.fluids import burgers
        >>> burgers(16).elements[0].check().prepare().size_
        16
    """
    spacing = 2 * math.pi / n
    centers = (np.arange(n) + 0.5) * spacing
    values = {'u': mean + amplitude * np.sin(centers)}
    space = LatticeSpace(n, spacing=spacing)
    return System(
        elements=[
            FieldElement('line', element=BurgersElement, space=space, values=values, discretization='finite-volume'),
        ],
        interactions=[ConservationInteraction('transport', target='line', state='u', flux=lambda s: s.target.u**2 / 2)],
    )

decaying_turbulence(shape=(64, 64), *, nu=0.001, modes=6, random_state=0)

Return decaying two-dimensional turbulence from a random vorticity field on a periodic square of side 2 pi.

Energy flows to large scales and enstrophy to small scales, where viscosity removes it, so vortices merge while the enstrophy decays faster than the energy.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of sites along both axes.	(64, 64)
nu	float	The kinematic viscosity.	0.001
modes	int	The largest wavenumber of the initial vorticity.	6
random_state	int | Generator | None	The seed or generator of the initial vorticity.	0

Returns:

Name	Type	Description
system	System	The model, with the vorticity flow.omega.

Examples:

>>> from skcomplex.models.physics.fluids import decaying_turbulence
>>> decaying_turbulence((16, 16)).elements[0].check().prepare().size_
256

Source code in src/skcomplex/models/physics/fluids/_flows.py

@register_model(
    field='physics/fluids',
    tags=('field', 'turbulence', 'spectral', 'pde'),
    references=('Kraichnan 1967', 'Satz 2022'),
    defaults={'until': 20.0, 'dt': 1.0, 'solver': ETDRK4Solver(dt=0.05)},
)
def decaying_turbulence(
    shape: Sequence[int] = (64, 64),
    *,
    nu: float = 1e-3,
    modes: int = 6,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return decaying two-dimensional turbulence from a random vorticity field on a periodic square of side 2 pi.

    Energy flows to large scales and enstrophy to small scales, where viscosity removes it, so vortices merge while
    the enstrophy decays faster than the energy.

    Args:
        shape:
            The number of sites along both axes.

        nu:
            The kinematic viscosity.

        modes:
            The largest wavenumber of the initial vorticity.

        random_state:
            The seed or generator of the initial vorticity.

    Returns:
        system:
            The model, with the vorticity `flow.omega`.

    Examples:
        >>> from skcomplex.models.physics.fluids import decaying_turbulence
        >>> decaying_turbulence((16, 16)).elements[0].check().prepare().size_
        256
    """
    shape = tuple(shape)
    vorticity = 5 * _smooth_noise(shape, modes, random_state)
    return incompressible_flow(shape, nu=nu, vorticity=vorticity)

discrete_velocity_boltzmann(n=200, *, c=1.0, kappa=1.0, length=10.0)

Return the Carleman discrete-velocity Boltzmann equation on a periodic line, from a pulse of right movers.

A pulse of right-moving particles on a uniform background splits into two moving pulses while collisions relax every site towards equal densities, and the total mass stays constant.

Parameters:

Name	Type	Description	Default
n	int	The number of cells.	200
c	float	The particle speed.	1.0
kappa	float	The collision rate.	1.0
length	float	The length of the line.	10.0

Returns:

Name	Type	Description
system	System	The model, with the densities gas.f_plus and gas.f_minus.

Examples:

>>> from skcomplex.models.physics.fluids import discrete_velocity_boltzmann
>>> discrete_velocity_boltzmann(16).elements[0].check().prepare().size_
16

Source code in src/skcomplex/models/physics/fluids/_kinetics.py

@register_model(
    field='physics/fluids',
    tags=('field', 'kinetic theory', 'discrete velocities', 'finite volume', 'pde'),
    references=('Carleman 1957', 'Gatignol 1975', 'Bellomo et al. 2017'),
    defaults={'until': 5.0, 'dt': 0.1, 'solver': SSPRKSolver()},
)
def discrete_velocity_boltzmann(
    n: int = 200,
    *,
    c: float = 1.0,
    kappa: float = 1.0,
    length: float = 10.0,
) -> System:
    """Return the Carleman discrete-velocity Boltzmann equation on a periodic line, from a pulse of right movers.

    A pulse of right-moving particles on a uniform background splits into two moving pulses while collisions relax
    every site towards equal densities, and the total mass stays constant.

    Args:
        n:
            The number of cells.

        c:
            The particle speed.

        kappa:
            The collision rate.

        length:
            The length of the line.

    Returns:
        system:
            The model, with the densities `gas.f_plus` and `gas.f_minus`.

    Examples:
        >>> from skcomplex.models.physics.fluids import discrete_velocity_boltzmann
        >>> discrete_velocity_boltzmann(16).elements[0].check().prepare().size_
        16
    """
    spacing = length / n
    centers = (np.arange(n) + 0.5) * spacing
    pulse = np.exp(-((centers - length / 2) ** 2))
    values = {'f_plus': 0.5 + pulse, 'f_minus': np.full(n, 0.5), 'c': c, 'kappa': kappa}
    space = LatticeSpace(n, spacing=spacing)
    return System(
        elements=[
            FieldElement('gas', element=CarlemanGasElement, space=space, values=values, discretization='finite-volume'),
        ],
        interactions=[
            ConservationInteraction(
                'right_movers',
                target='gas',
                state='f_plus',
                flux=lambda s: s.target.c * s.target.f_plus,
                source=_get_collisions,
            ),
            ConservationInteraction(
                'left_movers',
                target='gas',
                state='f_minus',
                flux=lambda s: -s.target.c * s.target.f_minus,
                source=lambda s: -_get_collisions(s),
            ),
        ],
    )

enskog(positions=8, velocities=12, *, knudsen=1.0, diameter=0.2, half_width=6.0)

Return the local Enskog equation, the Boltzmann equation with the contact value of the pair correlation of disks.

The collision frequency is multiplied at every position by Henderson's contact value (1 - 7 c / 16) / (1 - c)^2 of the area fraction c = pi n d^2 / 4, with n the density moment and d the diameter of the disks. The collision partners sit at one position, so the shift of one diameter of the full Enskog collision is left out.

Parameters:

Name	Type	Description	Default
positions	int	The number of positions.	8
velocities	int	The number of points of each velocity axis.	12
knudsen	float	The Knudsen number, the inverse of the collision frequency of the dilute gas.	1.0
diameter	float	The diameter of the disks.	0.2
half_width	float	The half-width of each velocity axis.	6.0

Returns:

Name	Type	Description
system	System	The model, with the density gas.f.

Examples:

>>> from skcomplex.models.physics.fluids import enskog
>>> enskog(4, 8).elements[0].check().prepare().size_
256

Source code in src/skcomplex/models/physics/fluids/_kinetic_equations.py

@register_model(
    field='physics/fluids',
    tags=('field', 'kinetic', 'phase space', 'collisions', 'pde'),
    references=('Bellomo et al. 2017', 'Enskog 1922'),
    defaults={
        'until': 0.5,
        'dt': 0.1,
        'solver': StrangSplittingSolver(SemiLagrangianSolver(0.1), RK4Solver(0.05), 0.1),
    },
)
def enskog(
    positions: int = 8,
    velocities: int = 12,
    *,
    knudsen: float = 1.0,
    diameter: float = 0.2,
    half_width: float = 6.0,
) -> System:
    """Return the local Enskog equation, the Boltzmann equation with the contact value of the pair correlation of disks.

    The collision frequency is multiplied at every position by Henderson's contact value `(1 - 7 c / 16) / (1 - c)^2` of
    the area fraction `c = pi n d^2 / 4`, with `n` the density moment and `d` the diameter of the disks. The collision
    partners sit at one position, so the shift of one diameter of the full Enskog collision is left out.

    Args:
        positions:
            The number of positions.

        velocities:
            The number of points of each velocity axis.

        knudsen:
            The Knudsen number, the inverse of the collision frequency of the dilute gas.

        diameter:
            The diameter of the disks.

        half_width:
            The half-width of each velocity axis.

    Returns:
        system:
            The model, with the density `gas.f`.

    Examples:
        >>> from skcomplex.models.physics.fluids import enskog
        >>> enskog(4, 8).elements[0].check().prepare().size_
        256
    """
    gas, frequency = _get_uniform_gas(positions, velocities, knudsen, half_width)
    area = math.pi * diameter**2 / 4

    def contact(s: SimpleNamespace) -> Expression:
        fraction = area * moment(s.candidate.f, 0, axes=VELOCITY_AXES)
        return frequency * (1 - 7 * fraction / 16) / (1 - fraction) ** 2

    stream = AdvectionInteraction('stream', target='gas', state='f', drift=_get_plane_streaming)
    collide = EncounterInteraction('collide', candidate='gas', cross_section=_get_maxwell_molecules, factor=contact)
    return System(elements=[gas], interactions=[stream, collide])

fhp_lattice_gas(shape=(64, 64), *, density=0.2, bump=0.6, radius=None, random_state=0)

Return the FHP-I lattice gas on a periodic triangular lattice, with a dense disc that spreads into the gas.

Every channel of every site is occupied at random with probability density, or bump inside a disc at the centre of the lattice, so the excess density spreads as a sound wave and then diffuses.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of rows and columns of the triangular lattice.	(64, 64)
density	float	The occupation probability of a channel outside the disc.	0.2
bump	float	The occupation probability of a channel inside the disc.	0.6
radius	float | None	The radius of the disc in lattice spacings, or None for a sixth of the smaller side.	None
random_state	int | Generator | None	The seed or generator of the initial occupations.	0

Returns:

Name	Type	Description
system	System	The model, with the channels gas.n0 to gas.n5 and the particle numbers gas.density.

Examples:

>>> from skcomplex.models.physics.fluids import fhp_lattice_gas
>>> fhp_lattice_gas((8, 8)).elements[0].check().prepare().size_
64

Source code in src/skcomplex/models/physics/fluids/_lattice_gases.py

@register_model(
    field='physics/fluids',
    tags=('lattice', 'lattice gas', 'stochastic', 'conservation'),
    references=('Frisch, Hasslacher, and Pomeau 1986', 'Dilão 2023'),
    defaults={'until': 200, 'dt': 1},
)
def fhp_lattice_gas(
    shape: Sequence[int] = (64, 64),
    *,
    density: float = 0.2,
    bump: float = 0.6,
    radius: float | None = None,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the FHP-I lattice gas on a periodic triangular lattice, with a dense disc that spreads into the gas.

    Every channel of every site is occupied at random with probability `density`, or `bump` inside a disc at the centre
    of the lattice, so the excess density spreads as a sound wave and then diffuses.

    Args:
        shape:
            The number of rows and columns of the triangular lattice.

        density:
            The occupation probability of a channel outside the disc.

        bump:
            The occupation probability of a channel inside the disc.

        radius:
            The radius of the disc in lattice spacings, or `None` for a sixth of the smaller side.

        random_state:
            The seed or generator of the initial occupations.

    Returns:
        system:
            The model, with the channels `gas.n0` to `gas.n5` and the particle numbers `gas.density`.

    Examples:
        >>> from skcomplex.models.physics.fluids import fhp_lattice_gas
        >>> fhp_lattice_gas((8, 8)).elements[0].check().prepare().size_
        64
    """
    space = LatticeSpace(tuple(shape), geometry='triangular').check().prepare()
    positions = space.positions_
    centre = positions.mean(axis=0)
    size = min(space.shape_) / 6 if radius is None else radius
    inside = np.linalg.norm(positions - centre, axis=1) < size
    probability = np.where(inside, bump, density)[:, None]
    occupied = (np.random.default_rng(random_state).random((space.size_, DIRECTIONS)) < probability).astype(float)
    channels = {name: occupied[:, direction] for direction, name in enumerate(CHANNELS)}
    gas = LatticeElement(
        'gas',
        element=FHPSiteElement,
        space=space,
        values={'density': occupied.sum(axis=1), **channels},
        visual=DENSITY_VISUAL,
    )
    return System(
        elements=[gas],
        interactions=[
            LocalRuleInteraction(
                'collision',
                target='gas',
                state=CHANNELS,
                rule=_collide,
                neighborhood=SITE,
                stochastic=True,
            ),
            LocalRuleInteraction('propagation', target='gas', state=CHANNELS, rule=_propagate, neighborhood=TRIANGULAR),
            LocalRuleInteraction(
                'count',
                target='gas',
                state='density',
                rule=_count,
                neighborhood=SITE,
                reads=CHANNELS,
            ),
        ],
    )

incompressible_flow(shape=(64, 64), size=(2 * math.pi, 2 * math.pi), *, nu=0.001, forcing=None, buoyancy=None, walls=False, vorticity=None, temperature=None)

Return incompressible two-dimensional flow in vorticity-streamfunction form on a periodic spectral lattice.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of sites along both axes.	(64, 64)
size	Sequence[float]	The lengths of the domain along both axes.	(2 * pi, 2 * pi)
nu	float	The kinematic viscosity, unused with buoyancy, whose viscosity is the Prandtl number.	0.001
forcing	ArrayLike | None	The forcing of the vorticity at every site, or None.	None
buoyancy	Mapping[str, float] | None	The Rayleigh number Ra and the Prandtl number Pr of Boussinesq convection with a temperature field, or None for a flow without temperature.	None
walls	bool	Whether the second axis is twice the height of a layer between free-slip walls in its middle and at its ends, the vorticity and temperature being odd about them.	False
vorticity	ArrayLike | None	The initial vorticity, or None for rest.	None
temperature	ArrayLike | None	The initial temperature departure of convection, or None for zero.	None

Returns:

Name	Type	Description
system	System	The flow, with the vorticity flow.omega and, with buoyancy, the temperature flow.theta.

Examples:

>>> from skcomplex.models.physics.fluids import incompressible_flow
>>> incompressible_flow((16, 16), buoyancy={'Ra': 800.0, 'Pr': 1.0}).elements[0].element.__name__
'BoussinesqElement'

Source code in src/skcomplex/models/physics/fluids/_flows.py

def incompressible_flow(
    shape: Sequence[int] = (64, 64),
    size: Sequence[float] = (2 * math.pi, 2 * math.pi),
    *,
    nu: float = 1e-3,
    forcing: npt.ArrayLike | None = None,
    buoyancy: Mapping[str, float] | None = None,
    walls: bool = False,
    vorticity: npt.ArrayLike | None = None,
    temperature: npt.ArrayLike | None = None,
) -> System:
    """Return incompressible two-dimensional flow in vorticity-streamfunction form on a periodic spectral lattice.

    Args:
        shape:
            The number of sites along both axes.

        size:
            The lengths of the domain along both axes.

        nu:
            The kinematic viscosity, unused with buoyancy, whose viscosity is the Prandtl number.

        forcing:
            The forcing of the vorticity at every site, or `None`.

        buoyancy:
            The Rayleigh number `Ra` and the Prandtl number `Pr` of Boussinesq convection with a temperature field,
            or `None` for a flow without temperature.

        walls:
            Whether the second axis is twice the height of a layer between free-slip walls in its middle and at its
            ends, the vorticity and temperature being odd about them.

        vorticity:
            The initial vorticity, or `None` for rest.

        temperature:
            The initial temperature departure of convection, or `None` for zero.

    Returns:
        system:
            The flow, with the vorticity `flow.omega` and, with buoyancy, the temperature `flow.theta`.

    Examples:
        >>> from skcomplex.models.physics.fluids import incompressible_flow
        >>> incompressible_flow((16, 16), buoyancy={'Ra': 800.0, 'Pr': 1.0}).elements[0].element.__name__
        'BoussinesqElement'
    """
    shape = tuple(int(item) for item in shape)
    space = LatticeSpace(shape, spacing=tuple(length / count for length, count in zip(size, shape, strict=True)))
    values: dict[str, object] = {}
    if vorticity is not None:
        values['omega'] = np.asarray(vorticity, dtype=float)
    if buoyancy is None:
        element: type[Element] = VorticityElement
        values['nu'] = nu
        if forcing is not None:
            values['f'] = np.asarray(forcing, dtype=float)
    else:
        element = BoussinesqElement
        values.update({key: float(value) for key, value in buoyancy.items()})
        if temperature is not None:
            values['theta'] = np.asarray(temperature, dtype=float)
    states = ['omega'] if buoyancy is None else ['omega', 'theta']
    symmetry = dict.fromkeys(states, (None, 'odd')) if walls else None
    field = FieldElement(
        'flow',
        element=element,
        space=space,
        values=values,
        discretization='spectral',
        symmetry=symmetry,
    )
    return System(elements=[field], interactions=_get_flow_interactions(buoyancy is not None))

kuramoto_sivashinsky(modes=128, *, length=32 * math.pi)

Return the Kuramoto-Sivashinsky equation u_t = -u u_x - u_xx - u_xxxx, spatiotemporally chaotic.

The field is spectral on a periodic line, started from cos(x / 16) (1 + sin(x / 16)) as in Kassam and Trefethen, and the stiff fourth derivative suits exponential time differencing.

The derivative of the nonlinear term u²/2 is dealiased by the two-thirds rule: it keeps only the wavenumbers below a third of 2 π modes / length, 2.67 for the defaults, so the spectrum drops to round-off from q = 2.69 on. The cut costs little, since the modes there decay at the rate q⁴ - q² above 45, and with the defaults the power at the last kept mode is five decades below the peak near q = 0.6. More modes on the same length move the cut up in proportion.

Parameters:

Name	Type	Description	Default
modes	int	The number of sites, and of Fourier modes.	128
length	float	The length of the line.	32 * pi

Returns:

Name	Type	Description
system	System	The model, with the field line.u.

Examples:

>>> from skcomplex.models.physics.fluids import kuramoto_sivashinsky
>>> kuramoto_sivashinsky(32).elements[0].discretization
'spectral'

Source code in src/skcomplex/models/physics/fluids/_flows.py

@register_model(
    field='physics/fluids',
    tags=('field', 'chaos', 'spectral', 'stiff', 'pde'),
    references=('Kuramoto 1978', 'Sivashinsky 1977', 'Kassam and Trefethen 2005'),
    defaults={'until': 150.0, 'dt': 1.0, 'solver': ETDRK4Solver(dt=0.25)},
)
def kuramoto_sivashinsky(modes: int = 128, *, length: float = 32 * math.pi) -> System:
    """Return the Kuramoto-Sivashinsky equation `u_t = -u u_x - u_xx - u_xxxx`, spatiotemporally chaotic.

    The field is spectral on a periodic line, started from `cos(x / 16) (1 + sin(x / 16))` as in Kassam and Trefethen,
    and the stiff fourth derivative suits exponential time differencing.

    The derivative of the nonlinear term `u²/2` is dealiased by the two-thirds rule: it keeps only the wavenumbers
    below a third of `2 π modes / length`, 2.67 for the defaults, so the spectrum drops to round-off from `q = 2.69`
    on. The cut costs little, since the modes there decay at the rate `q⁴ - q²` above 45, and with the defaults the
    power at the last kept mode is five decades below the peak near `q = 0.6`. More modes on the same length move the
    cut up in proportion.

    Args:
        modes:
            The number of sites, and of Fourier modes.

        length:
            The length of the line.

    Returns:
        system:
            The model, with the field `line.u`.

    Examples:
        >>> from skcomplex.models.physics.fluids import kuramoto_sivashinsky
        >>> kuramoto_sivashinsky(32).elements[0].discretization
        'spectral'
    """
    spacing = length / modes
    positions = np.arange(1, modes + 1) * spacing
    values = {'u': np.cos(positions / 16) * (1 + np.sin(positions / 16))}
    space = LatticeSpace(modes, spacing=spacing)
    transport = ConservationInteraction(
        'transport',
        target='line',
        state='u',
        flux=lambda s: s.target.u**2 / 2,
        riemann=None,
    )
    linear = RateInteraction(
        'linear',
        target='line',
        state='u',
        rate=lambda s: -laplacian(s.target.u) - laplacian(laplacian(s.target.u)),
    )
    return System(
        elements=[
            FieldElement(
                'line',
                element=KuramotoSivashinskyElement,
                space=space,
                values=values,
                discretization='spectral',
            ),
        ],
        interactions=[transport, linear],
    )

rayleigh_benard(Ra=1000.0, *, Pr=1.0, shape=(32, 32), amplitude=0.001)

Return Rayleigh-Bénard convection between free-slip walls, from a small perturbation of the critical mode.

The layer has height 1 and width 2 sqrt(2), the wavelength of the critical mode, and convection sets in above Ra = 27 pi^4 / 4, about 657.5. The lattice holds twice the layer's height with odd states, which gives the free-slip walls. The registered step ETDRK4Solver(dt=0.005) is stable up to about Ra = 4000, and faster convection needs a smaller step, such as ETDRK4Solver(dt=0.0025) at Ra = 5000.

Parameters:

Name	Type	Description	Default
Ra	float	The Rayleigh number.	1000.0
Pr	float	The Prandtl number.	1.0
shape	Sequence[int]	The number of sites along the width and twice the height.	(32, 32)
amplitude	float	The amplitude of the initial temperature perturbation.	0.001

Returns:

Name	Type	Description
system	System	The model, with the vorticity flow.omega and the temperature departure flow.theta.

Examples:

>>> from skcomplex.models.physics.fluids import rayleigh_benard
>>> rayleigh_benard(500.0).elements[0].values['Ra']
500.0

Source code in src/skcomplex/models/physics/fluids/_flows.py

@register_model(
    field='physics/fluids',
    tags=('field', 'convection', 'instability', 'spectral', 'pde'),
    references=('Rayleigh 1916', 'Chandrasekhar 1961', 'Mihailović et al. 2023'),
    defaults={'until': 10.0, 'dt': 0.5, 'solver': ETDRK4Solver(dt=0.005)},
)
def rayleigh_benard(
    Ra: float = 1000.0,
    *,
    Pr: float = 1.0,
    shape: Sequence[int] = (32, 32),
    amplitude: float = 1e-3,
) -> System:
    """Return Rayleigh-Bénard convection between free-slip walls, from a small perturbation of the critical mode.

    The layer has height 1 and width `2 sqrt(2)`, the wavelength of the critical mode, and convection sets in above
    `Ra = 27 pi^4 / 4`, about 657.5. The lattice holds twice the layer's height with odd states, which gives the
    free-slip walls. The registered step `ETDRK4Solver(dt=0.005)` is stable up to about `Ra = 4000`, and faster
    convection needs a smaller step, such as `ETDRK4Solver(dt=0.0025)` at `Ra = 5000`.

    Args:
        Ra:
            The Rayleigh number.

        Pr:
            The Prandtl number.

        shape:
            The number of sites along the width and twice the height.

        amplitude:
            The amplitude of the initial temperature perturbation.

    Returns:
        system:
            The model, with the vorticity `flow.omega` and the temperature departure `flow.theta`.

    Examples:
        >>> from skcomplex.models.physics.fluids import rayleigh_benard
        >>> rayleigh_benard(500.0).elements[0].values['Ra']
        500.0
    """
    width, height = 2 * math.sqrt(2), 2.0
    x = (np.arange(shape[0]) * width / shape[0])[:, None]
    z = (np.arange(shape[1]) * height / shape[1])[None, :]
    temperature = amplitude * np.cos(math.pi / math.sqrt(2) * x) * np.sin(math.pi * z)
    return incompressible_flow(
        shape,
        (width, height),
        buoyancy={'Ra': Ra, 'Pr': Pr},
        walls=True,
        temperature=temperature,
    )

shallow_water(n=400, *, length=10.0, upstream=2.0, downstream=1.0, g=9.81)

Return the shallow-water dam break: still water at two depths that a removed dam separated at the middle.

A shock runs downstream and a rarefaction upstream, with the depth and the discharge conserved by finite volumes and open edges.

Parameters:

Name	Type	Description	Default
n	int	The number of cells.	400
length	float	The length of the channel.	10.0
upstream	float	The depth upstream of the dam.	2.0
downstream	float	The depth downstream of the dam.	1.0
g	float	The gravitational acceleration.	9.81

Returns:

Name	Type	Description
system	System	The model, with the depth channel.h and the discharge channel.q.

Examples:

>>> from skcomplex.models.physics.fluids import shallow_water
>>> float(shallow_water(10).elements[0].values['h'][0])
2.0

Source code in src/skcomplex/models/physics/fluids/_flows.py

@register_model(
    field='physics/fluids',
    tags=('field', 'waves', 'shocks', 'conservation law', 'finite volume', 'pde'),
    references=('Stoker 1957', 'LeVeque 2002'),
    defaults={'until': 0.5, 'dt': 0.05, 'solver': SSPRKSolver()},
)
def shallow_water(
    n: int = 400,
    *,
    length: float = 10.0,
    upstream: float = 2.0,
    downstream: float = 1.0,
    g: float = 9.81,
) -> System:
    """Return the shallow-water dam break: still water at two depths that a removed dam separated at the middle.

    A shock runs downstream and a rarefaction upstream, with the depth and the discharge conserved by finite volumes
    and open edges.

    Args:
        n:
            The number of cells.

        length:
            The length of the channel.

        upstream:
            The depth upstream of the dam.

        downstream:
            The depth downstream of the dam.

        g:
            The gravitational acceleration.

    Returns:
        system:
            The model, with the depth `channel.h` and the discharge `channel.q`.

    Examples:
        >>> from skcomplex.models.physics.fluids import shallow_water
        >>> float(shallow_water(10).elements[0].values['h'][0])
        2.0
    """
    depth = np.where(np.arange(n) < n // 2, upstream, downstream).astype(float)
    space = LatticeSpace(n, spacing=length / n)
    values = {'h': depth, 'q': np.zeros(n), 'g': g}
    boundary = {'h': 'zero-flux', 'q': 'zero-flux'}
    field = FieldElement(
        'channel',
        element=ShallowWaterElement,
        space=space,
        boundary=boundary,
        values=values,
        discretization='finite-volume',
    )
    mass = ConservationInteraction('mass', target='channel', state='h', flux=lambda s: s.target.q)
    momentum = ConservationInteraction(
        'momentum',
        target='channel',
        state='q',
        flux=lambda s: s.target.q**2 / s.target.h + s.target.g * s.target.h**2 / 2,
    )
    return System(elements=[field], interactions=[mass, momentum])

vlasov_equation(positions=32, velocities=64, *, setup='landau', amplitude=0.01, half_width=6.0)

Return the Vlasov-Poisson equation of electrons over a neutralizing background, on a periodic position axis.

The electric field is solved from the density moment at every step. The set-up 'landau' starts from a Maxwellian with a density wave of wave number 0.5, whose electric energy decays at twice the Landau damping rate 0.1533. The set-up 'two-stream' starts from two beams of temperature 0.25 at the velocities -1.5 and 1.5, which are unstable at the wave number 0.5.

Parameters:

Name	Type	Description	Default
positions	int	The number of positions.	32
velocities	int	The number of points of the velocity axis.	64
setup	str	'landau' or 'two-stream'.	'landau'
amplitude	float	The relative amplitude of the initial density wave.	0.01
half_width	float	The half-width of the velocity axis.	6.0

Returns:

Name	Type	Description
system	System	The model, with the density gas.f.

Raises:

Type	Description
ValidationError	The set-up is unknown.

Examples:

>>> from skcomplex.models.physics.fluids import vlasov_equation
>>> vlasov_equation(8, 16).elements[0].check().prepare().size_
128

Source code in src/skcomplex/models/physics/fluids/_kinetic_equations.py

@register_model(
    field='physics/fluids',
    tags=('field', 'kinetic', 'phase space', 'plasma', 'pde'),
    references=('Bellomo et al. 2017', 'Cheng and Knorr 1976'),
    defaults={'until': 30.0, 'dt': 0.25, 'solver': SemiLagrangianSolver(0.25)},
)
def vlasov_equation(
    positions: int = 32,
    velocities: int = 64,
    *,
    setup: str = 'landau',
    amplitude: float = 0.01,
    half_width: float = 6.0,
) -> System:
    """Return the Vlasov-Poisson equation of electrons over a neutralizing background, on a periodic position axis.

    The electric field is solved from the density moment at every step. The set-up `'landau'` starts from a Maxwellian
    with a density wave of wave number 0.5, whose electric energy decays at twice the Landau damping rate 0.1533. The
    set-up `'two-stream'` starts from two beams of temperature 0.25 at the velocities -1.5 and 1.5, which are unstable
    at the wave number 0.5.

    Args:
        positions:
            The number of positions.

        velocities:
            The number of points of the velocity axis.

        setup:
            `'landau'` or `'two-stream'`.

        amplitude:
            The relative amplitude of the initial density wave.

        half_width:
            The half-width of the velocity axis.

    Returns:
        system:
            The model, with the density `gas.f`.

    Raises:
        ValidationError:
            The set-up is unknown.

    Examples:
        >>> from skcomplex.models.physics.fluids import vlasov_equation
        >>> vlasov_equation(8, 16).elements[0].check().prepare().size_
        128
    """
    if setup not in SETUPS:
        error_msg = f'Unknown set-up {setup!r} of the Vlasov equation. Use one of {", ".join(SETUPS)}.'
        raise ValidationError(error_msg)
    length = 4 * math.pi
    space = PhaseSpace(
        LatticeSpace((positions,), spacing=length / positions),
        [AxisSpec('v', velocities, -half_width, half_width)],
    )
    space.check().prepare()
    x = np.arange(positions) * length / positions
    v = space.coordinates_['v']
    profile = np.exp(-(v**2) / 2) / np.sqrt(2 * np.pi)
    if setup == 'two-stream':
        profile = (np.exp(-2 * (v - 1.5) ** 2) + np.exp(-2 * (v + 1.5) ** 2)) * np.sqrt(2 / np.pi) / 2
    gas = FieldElement(
        'gas',
        element=DensityElement,
        space=space,
        values={'f': np.outer(1.0 + amplitude * np.cos(0.5 * x), profile)},
        discretization='finite-volume',
        visual=VLASOV_VISUAL,
    )
    return System(
        elements=[gas],
        interactions=[AdvectionInteraction('vlasov', target='gas', state='f', drift=_get_poisson_drift)],
    )

wave_breaking(n=200, *, k0=1.0, amplitude=0.5, g=1.0)

Return the kinematic wave equation of deep-water wave trains on a periodic interval of length 2 pi.

The wave number starts as k0 (1 + amplitude sin x) and steepens into a jump, a break of the wavelength, at the time 1 / max(-C'(k0(x)) k0'(x)) of the first crossing of characteristics. Finite volumes conserve the total number of waves.

Parameters:

Name	Type	Description	Default
n	int	The number of cells.	200
k0	float	The mean initial wave number.	1.0
amplitude	float	The relative amplitude of the initial modulation, below 1.	0.5
g	float	The gravitational acceleration.	1.0

Returns:

Name	Type	Description
system	System	The model, with the wave number train.k.

Examples:

>>> from skcomplex.models.physics.fluids import wave_breaking
>>> wave_breaking(16).elements[0].check().prepare().size_
16

Source code in src/skcomplex/models/physics/fluids/_kinetics.py

@register_model(
    field='physics/fluids',
    tags=('field', 'waves', 'breaking', 'conservation law', 'finite volume', 'pde'),
    references=('Whitham 1974', 'Forrest 2018'),
    defaults={'until': 20.0, 'dt': 0.5, 'solver': SSPRKSolver()},
)
def wave_breaking(n: int = 200, *, k0: float = 1.0, amplitude: float = 0.5, g: float = 1.0) -> System:
    """Return the kinematic wave equation of deep-water wave trains on a periodic interval of length 2 pi.

    The wave number starts as `k0 (1 + amplitude sin x)` and steepens into a jump, a break of the wavelength, at the
    time `1 / max(-C'(k0(x)) k0'(x))` of the first crossing of characteristics. Finite volumes conserve the total
    number of waves.

    Args:
        n:
            The number of cells.

        k0:
            The mean initial wave number.

        amplitude:
            The relative amplitude of the initial modulation, below 1.

        g:
            The gravitational acceleration.

    Returns:
        system:
            The model, with the wave number `train.k`.

    Examples:
        >>> from skcomplex.models.physics.fluids import wave_breaking
        >>> wave_breaking(16).elements[0].check().prepare().size_
        16
    """
    spacing = 2 * math.pi / n
    centers = (np.arange(n) + 0.5) * spacing
    values = {'k': k0 * (1 + amplitude * np.sin(centers)), 'g': g}
    space = LatticeSpace(n, spacing=spacing)
    return System(
        elements=[
            FieldElement(
                'train',
                element=WaveNumberElement,
                space=space,
                values=values,
                discretization='finite-volume',
            ),
        ],
        interactions=[
            ConservationInteraction('action', target='train', state='k', flux=lambda s: sqrt(s.target.g * s.target.k)),
        ],
    )