Toner-Tu hydrodynamics
The Toner-Tu hydrodynamics of a polar flock on a periodic lattice, from a disordered velocity field.
- Name:
toner_tu - Field:
physics/active_matter - Formalism: stochastic differential equations
- References: Toner and Tu 1995; Marchetti et al. 2013
Factory
@register_model(
field='physics/active_matter',
tags=('field', 'flocking', 'hydrodynamics', 'pde'),
references=('Toner and Tu 1995', 'Marchetti et al. 2013'),
defaults={'until': 20.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.02)},
)
def toner_tu(
shape: Sequence[int] = (32, 32),
*,
rho0: float = 1.0,
rho_c: float = 0.5,
a: float = 1.0,
beta: float = 1.0,
lam: float = 1.0,
sigma: float = 0.5,
D: float = 0.5,
spacing: float = 1.0,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return the Toner-Tu hydrodynamics of a polar flock on a periodic lattice, from a disordered velocity field.
Above the critical density the velocity grows and aligns until the whole fluid moves together with the speed
`sqrt(a (rho0 - rho_c) / beta)`.
Args:
shape:
The number of sites along each axis.
rho0:
The mean density.
rho_c:
The critical density of the ordering transition.
a:
The rate of the ordering instability per unit density.
beta:
The saturation of the speed.
lam:
The strength of the self-advection.
sigma:
The pressure per unit density.
D:
The viscosity of the velocity.
spacing:
The distance between sites.
random_state:
The seed or generator of the initial velocities.
Returns:
system:
The model, with the fields `flock.rho`, `flock.vx`, and `flock.vy`.
Examples:
>>> from skcomplex.models.physics.active_matter import toner_tu
>>> toner_tu((8, 8)).elements[0].check().prepare().size_
64
"""
shape = tuple(shape)
rng = np.random.default_rng(random_state)
values = {
'rho': np.full(shape, rho0),
'vx': rng.normal(0.0, 0.1, shape),
'vy': rng.normal(0.0, 0.1, shape) + 0.05,
}
parameters = {'a': a, 'rho_c': rho_c, 'beta': beta, 'lam': lam, 'sigma': sigma, 'D': D}
space = LatticeSpace(shape, spacing=spacing)
return System(
elements=[FieldElement('flock', element=PolarFluidElement, space=space, values={**values, **parameters})],
)
Components
class PolarFluidElement(Element):
"""Site of the Toner-Tu polar active fluid, with a density `rho` and a velocity `(vx, vy)`.
The density is conserved, `drho/dt + div(rho v) = 0`, and the velocity follows
`dv/dt + lam (v . grad) v = a (rho - rho_c) v - beta |v|^2 v - sigma grad(rho) + D laplacian(v)` with additive
noise of intensity `noise`. Above the critical density the fluid moves with the speed
`sqrt(a (rho - rho_c) / beta)`.
Args:
label:
The label of the site or the field.
visual:
How the element is drawn, or `None` for `default_visual_`.
rho:
The initial density.
vx:
The initial velocity along the first axis.
vy:
The initial velocity along the second axis.
a:
The rate of the ordering instability per unit density.
rho_c:
The critical density of the ordering transition.
beta:
The saturation of the speed.
lam:
The strength of the self-advection.
sigma:
The pressure per unit density.
D:
The viscosity of the velocity.
noise:
The intensity of the velocity noise.
Examples:
>>> from skcomplex.models.physics.active_matter import PolarFluidElement
>>> PolarFluidElement('site').rho_c
0.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
rho: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
vx: Annotated[float, StateSpec()] = 0.0,
vy: Annotated[float, StateSpec()] = 0.0,
a: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
rho_c: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
beta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
lam: Annotated[float, ParameterSpec()] = 1.0,
sigma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
noise: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.rho = rho
self.vx = vx
self.vy = vy
self.a = a
self.rho_c = rho_c
self.beta = beta
self.lam = lam
self.sigma = sigma
self.D = D
self.noise = noise
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the continuity equation and the Toner-Tu equation of each velocity component.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
velocity = [s.vx, s.vy]
pressure = gradient(s.sigma * s.rho, 2)
speed_squared = s.vx**2 + s.vy**2
laws: list[BaseLaw] = [RateLaw(s.rho, -divergence([s.rho * s.vx, s.rho * s.vy]))]
for axis, component in enumerate(velocity):
transport = advection(component, [s.lam * item for item in velocity])
growth = (s.a * (s.rho - s.rho_c) - s.beta * speed_squared) * component
laws += [
RateLaw(component, -transport + growth - pressure[axis] + s.D * laplacian(component)),
NoiseLaw(component, sqrt(2 * s.noise)),
]
return laws
Equations
flock.rho[i]: diff(rho_flock(t), t) == -partial(rho_flock(t)*vx_flock(t), 0) - partial(rho_flock(t)*vy_flock(t), 1)
flock.vx[i]: diff(vx_flock(t), t) == D_flock*laplacian(vx_flock(t)) - ((vx_flock(t)^2 + vy_flock(t)^2)*beta_flock + a_flock*(rho_c_flock - rho_flock(t)))*vx_flock(t) - partial(sigma_flock*rho_flock(t), 0) - upwind(vx_flock(t), 1, lam_flock*vy_flock(t)) - upwind(vx_flock(t), 0, lam_flock*vx_flock(t))
flock.vy[i]: diff(vy_flock(t), t) == D_flock*laplacian(vy_flock(t)) - ((vx_flock(t)^2 + vy_flock(t)^2)*beta_flock + a_flock*(rho_c_flock - rho_flock(t)))*vy_flock(t) - partial(sigma_flock*rho_flock(t), 1) - upwind(vy_flock(t), 1, lam_flock*vy_flock(t)) - upwind(vy_flock(t), 0, lam_flock*vx_flock(t))