Bertin-Droz-Grégoire kinetic flocking
The kinetic theory of flocking, whose polar order grows above a transition density.
- Name:
bertin_droz_gregoire - Field:
physics/active_matter - Formalism: ordinary differential equations
- References: Bertin, Droz, and Grégoire 2006; Bertin, Droz, and Grégoire 2009
Factory
@register_model(
field='physics/active_matter',
tags=('field', 'kinetic', 'phase space', 'flocking', 'pde'),
references=('Bertin, Droz, and Grégoire 2006', 'Bertin, Droz, and Grégoire 2009'),
defaults={'until': 0.5, 'dt': 0.1, 'solver': RK4Solver(dt=0.05)},
)
def bertin_droz_gregoire(
positions: int = 4,
angles: int = 16,
*,
density: float = 1.0,
speed: float = 1.0,
interaction_range: float = 0.5,
tumbling: float = 1.0,
noise: float = 0.5,
perturbation: float = 0.01,
) -> System:
"""Return the kinetic theory of flocking, whose polar order grows above a transition density.
Particles stream along their orientation, tumble at a rate by a wrapped normal noise, and align in binary collisions
at their mean angle with the same noise. A homogeneous disordered state is unstable above the density
`pi tumbling (1 - exp(-noise^2 / 2)) / (8 interaction_range speed (exp(-noise^2 / 2) - 2 / 3))`.
Args:
positions:
The number of positions along each of the two periodic position axes of unit spacing.
angles:
The number of cells of the orientation axis.
density:
The density of the homogeneous state.
speed:
The speed of the particles.
interaction_range:
The interaction range.
tumbling:
The rate of the tumbling events.
noise:
The standard deviation of the noise of tumbling and alignment.
perturbation:
The relative amplitude of the initial polar perturbation.
Returns:
system:
The model, with the density `flock.f`.
Examples:
>>> from skcomplex.models.physics.active_matter import bertin_droz_gregoire
>>> bertin_droz_gregoire(2, 8).elements[0].check().prepare().size_
32
"""
axis = AxisSpec('theta', angles, 0.0, 2 * math.pi, periodic=True)
space = PhaseSpace(LatticeSpace((positions, positions)), [axis]).check().prepare()
theta = space.coordinates_['theta']
profile = density * (1.0 + perturbation * np.cos(theta)) / (2 * math.pi)
values = np.broadcast_to(profile, space.shape_).copy()
flock = FieldElement(
'flock',
element=DensityElement,
space=space,
values={'f': values},
discretization='finite-volume',
)
def drift(_: SimpleNamespace) -> list[Expression]:
return [speed * cos(coordinate('theta')), speed * sin(coordinate('theta')), 0]
return System(
elements=[flock],
interactions=[
AdvectionInteraction('stream', target='flock', state='f', drift=drift),
AngularDiffusionInteraction('tumble', candidate='flock', angles=angles, tumbling=tumbling, noise=noise),
AngularAlignmentInteraction(
'align',
candidate='flock',
field='flock',
angles=angles,
frequency=interaction_range * speed,
noise=noise,
),
],
)
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 AngularDiffusionInteraction(EncounterInteraction):
"""Self-propelled particles change their angle at a rate by a wrapped normal noise, a linear redistribution.
Args:
label:
The label of the interaction, unique among the parts of its system.
candidate:
The path of the field on a phase space with a periodic orientation axis.
angles:
The number of cells of the orientation axis on `[0, 2 pi)`.
tumbling:
The rate of the events.
noise:
The standard deviation of the noise of an event.
state:
The name of the state of the density.
Examples:
>>> from skcomplex.models.physics.active_matter import AngularDiffusionInteraction
>>> AngularDiffusionInteraction('tumble', candidate='flock', angles=8, tumbling=1.0, noise=0.5).check().field
"""
def __init__(
self,
label: str,
*,
candidate: Annotated[str, RoleSpec()],
angles: int = 32,
tumbling: float = 1.0,
noise: float = 0.5,
state: str = 'f',
) -> None:
super().__init__(label, candidate=candidate, rate=self._rate, kernel=self._kernel, state=state)
self.angles = angles
self.tumbling = tumbling
self.noise = noise
def _rate(self, _: SimpleNamespace) -> float:
"""Return the rate of the events."""
return self.tumbling
def _kernel(self, _: SimpleNamespace) -> npt.NDArray[np.float64]:
"""Return the probability of every angle after an event at every angle."""
centres = 2 * np.pi * (np.arange(self.angles) + 0.5) / self.angles
gaps = np.angle(np.exp(1j * np.subtract.outer(centres, centres)))
weights = np.exp(-(gaps**2) / (2 * self.noise**2))
return weights / weights.sum(axis=1, keepdims=True)
class AngularAlignmentInteraction(EncounterInteraction):
"""Pairs of self-propelled particles leave a collision at their mean angle plus a wrapped normal noise.
A pair of angles meets at the rate `4 d0 v0 |sin((a - b) / 2)|`, the frequency of Bertin, Droz, and Grégoire, with
`d0` the interaction range and `v0` the speed. A pair of opposite angles leaves at either of its two mean angles
with equal probability.
Args:
label:
The label of the interaction, unique among the parts of its system.
candidate:
The path of the field on a phase space with a periodic orientation axis.
field:
The path of the field the candidate meets, usually the candidate itself.
angles:
The number of cells of the orientation axis on `[0, 2 pi)`.
frequency:
The product `d0 v0` of the interaction range and the speed.
noise:
The standard deviation of the noise of the outcome.
state:
The name of the state of the density.
Examples:
>>> from skcomplex.models.physics.active_matter import AngularAlignmentInteraction
>>> align = AngularAlignmentInteraction('align', candidate='flock', field='flock', frequency=0.5, noise=0.5)
>>> align.check().outcome
'conservative'
"""
def __init__(
self,
label: str,
*,
candidate: Annotated[str, RoleSpec()],
field: Annotated[str, RoleSpec()],
angles: int = 32,
frequency: float = 1.0,
noise: float = 0.5,
state: str = 'f',
) -> None:
super().__init__(label, candidate=candidate, field=field, rate=self._rate, kernel=self._kernel, state=state)
self.angles = angles
self.frequency = frequency
self.noise = noise
def _centres(self) -> npt.NDArray[np.float64]:
"""Return the angle of every cell."""
return 2 * np.pi * (np.arange(self.angles) + 0.5) / self.angles
def _rate(self, _: SimpleNamespace) -> npt.NDArray[np.float64]:
"""Return the rate of every pair of angles, times the width of the cell of the partner."""
centres = self._centres()
return 8 * np.pi / self.angles * self.frequency * np.abs(np.sin(np.subtract.outer(centres, centres) / 2))
def _noise_weights(
self,
centres: npt.NDArray[np.float64],
means: npt.NDArray[np.float64],
) -> npt.NDArray[np.float64]:
"""Return the wrapped normal weight of every angle around the mean angle of every pair."""
gaps = np.angle(np.exp(1j * (centres[None, None, :] - means[:, :, None])))
return np.exp(-(gaps**2) / (2 * self.noise**2))
def _kernel(self, _: SimpleNamespace) -> npt.NDArray[np.float64]:
"""Return the probability of every outcome angle of every pair of angles, around the mean angle of the pair."""
centres = self._centres()
cosines, sines = np.cos(centres), np.sin(centres)
along, across = np.add.outer(cosines, cosines), np.add.outer(sines, sines)
means = np.arctan2(across, along)
opposite = np.hypot(along, across) < ANTIPODAL
means = np.where(opposite, centres[:, None] + np.pi / 2, means)
weights = self._noise_weights(centres, means)
weights = np.where(opposite[:, :, None], (weights + self._noise_weights(centres, means + np.pi)) / 2, weights)
return weights / weights.sum(axis=2, keepdims=True)
Equations
flock.f[i]: diff(f_flock(t), t) == 0