Vlasov
The Vlasov-Poisson equation of electrons over a neutralizing background, on a periodic position axis.
- Name:
vlasov_equation - Field:
physics/fluids - Formalism: ordinary differential equations
- References: Bellomo et al. 2017; Cheng and Knorr 1976
Factory
@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)],
)
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
Equations
gas.f[i]: diff(f_gas(t), t) == 0