BGK
The BGK equation of the Sod shock tube over one position and one velocity axis.
- Name:
bgk_equation - Field:
physics/fluids - Formalism: ordinary differential equations
- References: Bhatnagar, Gross, and Krook 1954; Mieussens 2000; Sod 1978
Factory
@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])
Components
class DensityElement(Element):
"""Value of a density at one class, a member of a distribution over the classes.
A distribution is a group of densities on the sites of a lattice of classes, one member per class, and encounter,
proliferation, and destruction laws change it. A density has no laws of its own.
Args:
label:
The label of the density, unique among the parts of its system.
visual:
How the element is drawn, or `None` for `default_visual_`.
f:
The initial density, a nonnegative number.
Attributes:
default_visual_:
The visual of the element when it sets none.
Examples:
>>> from skcomplex.core.elements import DensityElement, GroupElement
>>> from skcomplex.core.spaces import LatticeSpace
>>> classes = LatticeSpace(5, boundary='fixed')
>>> group = GroupElement('distribution', element=DensityElement, space=classes, on='sites', values={'f': 0.2})
>>> group.check().prepare().size_
5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
f: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.f = f
class AdvectionInteraction(ConservationInteraction):
"""Advection `du/dt + divergence(a u) = source(ns)` of a state `u` of a field by a drift field `a(ns)`.
The flux of the conservation law is the drift times the state, written in conservative form, so the total of the
state on a periodic field stays constant and, on a spectral field, the dealiasing of products applies to it. The
drift is one expression for a field of one dimension, or one per dimension, of the states and parameters of the
field, and may read other states than the one it carries, such as the flow of a fluid derived from a
streamfunction, or the drift of a probability density in a Fokker-Planck equation.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field whose state is carried.
state:
The name of the carried state of the field.
drift:
A function of the namespace that returns the drift, or one drift per dimension.
source:
A function of the namespace that returns the source, or `None` for none.
limiter:
The slope limiter of the reconstruction of a finite-volume field, `'minmod'`, `'van-leer'`, or `None`.
riemann:
The numerical flux of a finite-volume field, `'rusanov'` or `'godunov'`, or `None` for the divergence of
the flux on a field of another discretization.
sources:
The paths of other elements the drift and the source read, as the namespaces `sources`. A field over the
position lattice of a target on a phase space is read broadcast along the class axes.
Raises:
ValidationError:
The state is not a name, a function is not callable, or the limiter or the numerical flux is unknown.
Examples:
>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element, FieldElement
>>> from skcomplex.core.interactions import AdvectionInteraction
>>> from skcomplex.core.spaces import LatticeSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class CellElement(Element):
... def __init__(self, label, u: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.u = u
>>> space = LatticeSpace((50,), spacing=0.1, boundary='periodic')
>>> values = {'u': np.exp(-((np.arange(50) - 25) / 5.0) ** 2)}
>>> volumes = 'finite-volume'
>>> field = FieldElement('tracer', element=CellElement, space=space, values=values, discretization=volumes)
>>> drift = AdvectionInteraction('drift', target='tracer', state='u', drift=lambda s: 1.0)
>>> system = System(elements=[field], interactions=[drift])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=1.0).trajectory_
>>> total = trajectory.values('tracer.u')[-1].sum()
>>> round(float(total), 6) == round(float(values['u'].sum()), 6)
True
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'u',
drift: Callable[[SimpleNamespace], object] | None = None,
source: Callable[[SimpleNamespace], object] | None = None,
limiter: str | None = 'minmod',
riemann: str | None = 'rusanov',
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
) -> None:
super().__init__(
label,
target=target,
state=state,
source=source,
limiter=limiter,
riemann=riemann,
sources=sources,
)
self.drift = drift
def _check(self) -> None:
"""Raise an error for a state that is not a name, a drift that is not callable, or an unknown scheme.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is not callable, or the limiter or the
numerical flux is unknown.
"""
super()._check()
if not callable(self.drift):
error_msg = f'The drift of {self.label!r} is {self.drift!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
def _get_flux(self, namespace: SimpleNamespace) -> object:
"""Return the flux `a u`, one component per component of the drift.
Args:
namespace:
The namespace of the symbols the flux can use, with the namespace of the field.
Returns:
flux:
The flux, one expression or one per dimension.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
drift = self.drift(namespace) if callable(self.drift) else None
if isinstance(drift, FreeModuleElement | list | tuple):
return [SR(component) * state for component in drift]
return SR(drift) * state
class BhatnagarGrossKrookInteraction(RelaxationInteraction):
"""Relaxation of Bhatnagar, Gross, and Krook of a density towards the maximum-entropy density with its moments.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field on a phase space.
state:
The name of the state of the density.
axes:
The names of the velocity axes.
time:
A function of the namespace that returns the relaxation time, or a positive number.
Examples:
>>> from skcomplex.models.physics.fluids import BhatnagarGrossKrookInteraction
>>> BhatnagarGrossKrookInteraction('bgk', target='gas', axes=('v',), time=0.1).check().axes
('v',)
"""
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
def _get_goal(self, namespace: SimpleNamespace) -> object:
"""Return the maximum-entropy density with the moments of the state of the target."""
return maximum_entropy(getattr(namespace.target, self.state), axes=self.axes)
Equations
gas.f[i]: diff(f_gas(t), t) == -100.0*f_gas(t) + 100.0*maximum_entropy(f_gas(t), 2, _coordinate_v)