Boltzmann
The Boltzmann equation of Maxwell molecules over one position and two velocity axes.
- Name:
boltzmann_equation - Field:
physics/fluids - Formalism: ordinary differential equations
- References: Bellomo et al. 2017; Mouhot and Pareschi 2006
Factory
@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])
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 EncounterInteraction(Interaction):
"""Encounters of a candidate density with a field density over activity classes, which move, remove, or create.
The candidate and the field are groups of `DensityElement` on the sites of lattices of activity classes, one member
per class. A candidate at class `h` meets a field at class `k` at the rate `eta[h, k]`. In the outcome
`'conservative'` the candidate moves to class `i` with the probability `kernel[h, k, i]`, in the outcome
`'destructive'` it is removed, and in the outcome `'proliferative'` it creates a new candidate at class `i` with
the probability `kernel[h, k, i]`, or at its own class without a kernel. The rate and the kernel are functions of
the namespace, which holds `namespace.candidate`, `namespace.field`, and the parameters of the interaction.
On a `PhaseSpace` the classes are the points of the class axes, and encounters pair classes at the same position. A
rate that is one expression may read the `coordinate` of the class of the candidate. Without a field the encounter
is a linear redistribution: a conservative candidate at class `h` moves to class `i` at the rate `eta[h]` with the
probability `kernel[h, i]`.
With a `cross_section` the candidate is a field on a `PhaseSpace` with two interval class axes, and its rate at
every position is the binary collision integral over those axes, evaluated by the fast spectral method on a padded
grid.
The density is conserved exactly, and the momentum and the energy within the spectral accuracy of the method.
Args:
label:
The label of the interaction, unique among the parts of its system.
candidate:
The path of the group of the candidate density.
field:
The path of the group of the field density, or `None` for a linear redistribution.
rate:
The function `rate(namespace)` of the encounter rate, a number or an expression, or an array `eta[h, k]`
of them, or `eta[h]` without a field, or `None` with a cross-section.
kernel:
The function `kernel(namespace)` of the array `kernel[h, k, i]` of the probabilities of the class of the
moved or created candidate, or `kernel[h, i]` without a field, or `None` for a destructive encounter or the
candidate's own class.
outcome:
`'conservative'`, `'destructive'`, or `'proliferative'`.
state:
The name of the state of the densities.
cross_section:
The function `B(speed, cosine)` of arrays of the relative speed and the cosine of the deflection angle, or
`None` for encounters at a rate.
factor:
A function of the namespace that returns a factor of the collision frequency, such as a function of the
density moment, or `None` for none.
Raises:
ValidationError:
The outcome is unknown, the state is no name, a function is not callable, a conservative encounter has no
kernel, an encounter without a field is not conservative, an encounter has neither a rate nor a
cross-section, or a cross-section comes with a rate, a kernel, a field, or an outcome.
Examples:
>>> import numpy as np
>>> from skcomplex.core.interactions import EncounterInteraction
>>> kernel = np.full((3, 3, 3), 1 / 3)
>>> mixing = EncounterInteraction(
... 'mixing', candidate='population', field='population', rate=lambda s: 1.0, kernel=lambda s: kernel
... )
>>> mixing.check().outcome
'conservative'
>>> collide = EncounterInteraction('collide', candidate='gas', cross_section=lambda speed, cosine: 1.0)
>>> collide.check().rate is None
True
"""
def __init__(
self,
label: str,
*,
candidate: Annotated[str, RoleSpec()],
field: Annotated[str | None, RoleSpec(optional=True)] = None,
rate: Callable[[SimpleNamespace], object] | None = None,
kernel: Callable[[SimpleNamespace], object] | None = None,
outcome: str = 'conservative',
state: str = 'f',
cross_section: Callable[..., object] | None = None,
factor: Callable[[SimpleNamespace], object] | None = None,
) -> None:
super().__init__(label)
self.candidate = candidate
self.field = field
self.rate = rate
self.kernel = kernel
self.outcome = outcome
self.state = state
self.cross_section = cross_section
self.factor = factor
def _check_collision(self) -> None:
"""Raise an error for a cross-section or a factor that is no function, or a cross-section with others."""
if not callable(self.cross_section):
error_msg = f'The cross-section of {self.label!r} is {self.cross_section!r}. Give a function of two arrays.'
raise ValidationError(error_msg)
if self.factor is not None and not callable(self.factor):
error_msg = f'The factor of {self.label!r} is {self.factor!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
others = self.rate, self.kernel, self.field
if any(item is not None for item in others) or self.outcome != 'conservative':
error_msg = f'The cross-section of {self.label!r} sets the rate. Give no rate, kernel, field, or outcome.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for an unknown outcome, a function that is not callable, or a conservative kernel missing.
Raises:
ValidationError:
The outcome is unknown, the state is no name, a function is not callable, a conservative encounter has
no kernel, an encounter without a field is not conservative, an encounter has neither a rate nor a
cross-section, or a cross-section comes with a rate, a kernel, a field, or an outcome.
"""
super()._check()
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = f'{type(self).__name__} {self.label!r} reads the state {self.state!r}, which is no name.'
raise ValidationError(error_msg)
if self.cross_section is not None:
self._check_collision()
return
if self.factor is not None:
error_msg = f'The factor of {self.label!r} scales a cross-section. Give `cross_section`.'
raise ValidationError(error_msg)
if self.outcome not in OUTCOMES:
error_msg = (
f'The outcome of encounter {self.label!r} is {self.outcome!r}. Use one of {", ".join(OUTCOMES)}.'
)
raise ValidationError(error_msg)
if not callable(self.rate) or (self.kernel is not None and not callable(self.kernel)):
error_msg = f'The rate and the kernel of encounter {self.label!r} must be functions of the namespace.'
raise ValidationError(error_msg)
if self.outcome == 'conservative' and self.kernel is None:
error_msg = f'The conservative encounter {self.label!r} moves candidates by a kernel. Give `kernel`.'
raise ValidationError(error_msg)
if self.field is None and self.outcome != 'conservative':
error_msg = f'The {self.outcome} encounter {self.label!r} reads a density of others. Give `field`.'
raise ValidationError(error_msg)
def _linear_laws(self, candidate: Expression, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the encounter law of a redistribution without a field, one field class of value one."""
assert self.rate is not None
rate = np.asarray(self.rate(namespace), dtype=object)
rates = rate if rate.ndim == 0 else rate.reshape(-1, 1)
kernel = np.asarray(self.kernel(namespace) if self.kernel is not None else (), dtype=object)
return [EncounterLaw(candidate, None, rates[()], np.expand_dims(kernel, axis=1), label=self.label)]
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the encounter, bilinear, or collision law of the candidate density.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.candidate` and `namespace.field`.
Returns:
laws:
The law of the encounters.
"""
candidate = getattr(self._get_member(namespace, 'candidate'), self.state)
if self.cross_section is not None:
factor = SR(1 if self.factor is None else self.factor(namespace))
return [
RateLaw(candidate, factor * collision(candidate), cross_section=self.cross_section, label=self.label),
]
if self.field is None:
return self._linear_laws(candidate, namespace)
field = getattr(self._get_member(namespace, 'field'), self.state)
assert self.rate is not None
rate = self.rate(namespace)
kernel = None if self.kernel is None else self.kernel(namespace)
if self.outcome == 'destructive':
loss = -np.asarray(rate, dtype=object)
return [BilinearLaw(candidate, candidate, field, loss, label=self.label)]
if self.outcome == 'proliferative':
return [BilinearLaw(candidate, candidate, field, rate, kernel, label=self.label)]
return [EncounterLaw(candidate, field, rate, kernel, label=self.label)]
Equations
gas.f[i]: diff(f_gas(t), t) == 1.0*collision(f_gas(t))