Dry active nematics
A dry active nematic, a conserved density coupled to a symmetric traceless alignment tensor.
- Name:
dry_active_nematics - Field:
physics/active_matter - Formalism: ordinary differential equations
- References: Marchetti et al. 2013; arXiv:1207.2929
Factory
@register_model(
field='physics/active_matter',
tags=('field', 'active matter', 'nematic', 'tensor', 'liquid crystal', 'giant number fluctuations', 'pde'),
references=('Marchetti et al. 2013', 'arXiv:1207.2929'),
defaults={'until': 100.0, 'dt': 5.0, 'solver': RK4Solver(dt=0.02), 'backend': NumbaBackend()},
)
def dry_active_nematics(
shape: Sequence[int] = (64, 64),
*,
density: float = 1.5,
critical_density: float = 1.0,
a: float = 1.0,
c: float = 2.0,
K: float = 1.0,
gamma: float = 1.0,
D: float = 0.5,
current: float = 0.5,
beta: float = 0.5,
perturbation: float = 0.01,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return a dry active nematic, a conserved density coupled to a symmetric traceless alignment tensor.
On a periodic lattice the density and the tensor follow Section II.2 of Marchetti et al.,
rho_t = -div(J), J = -D grad(rho) - current div(Q),
Q_t = H / gamma + beta (partial_i partial_j rho - delta_ij laplacian(rho) / 2),
with the molecular field `H` of the Landau-de Gennes energy, whose coefficient `a (1 - rho / rho_c)` turns negative
above the critical density `rho_c`. The current along `div(Q)` is the active current, which no free energy gives.
Args:
shape:
The number of sites along every axis of the periodic lattice.
density:
The mean density, which the dynamics conserve.
critical_density:
The density at which the isotropic state loses stability.
a:
The scale of the coefficient of `tr(Q^2) / 2` of the energy.
c:
The coefficient of `tr(Q^2)^2 / 4` of the energy.
K:
The elastic constant.
gamma:
The rotational viscosity.
D:
The diffusion coefficient of the density.
current:
The coefficient of the active current along `div(Q)`.
beta:
The coupling of the tensor to the curvature of the density.
perturbation:
The standard deviation of the initial noise of the density and the tensor.
random_state:
The seed or generator of the initial noise.
Returns:
system:
The model, with the fields `nematic.rho` and `nematic.Q`.
Examples:
>>> from skcomplex.models.physics.active_matter import dry_active_nematics
>>> dry_active_nematics((8, 8)).elements[0].values['beta']
0.5
"""
shape = tuple(shape)
rng = np.random.default_rng(random_state)
sites = int(np.prod(shape))
values = {
'rho': density + perturbation * rng.standard_normal(sites),
'Q': perturbation * rng.standard_normal((sites, 2)),
'beta': beta,
}
nematic = FieldElement('nematic', element=NematicDensitySiteElement, space=LatticeSpace(shape), values=values)
energy = LandauDeGennesEnergyInteraction(
'energy',
target='nematic',
a=a,
c=c,
K=K,
gamma=gamma,
quadratic=lambda s: s.a * (1 - s.target.rho / critical_density),
)
def flux(s: SimpleNamespace) -> list[Expression]:
"""Return the current of the density, its diffusion and its active current."""
active = tensor_divergence(s.target.Q)
return [-D * item - current * other for item, other in zip(gradient(s.target.rho), active, strict=True)]
transport = ConservationInteraction(
'transport',
target='nematic',
state='rho',
flux=flux,
limiter=None,
riemann=None,
)
return System(elements=[nematic], interactions=[energy, transport])
Components
class NematicDensitySiteElement(Element):
"""Site of a dry active nematic, whose alignment tensor grows with the traceless curvature of the density.
The rate of the tensor from the density is `beta (partial_i partial_j rho - delta_ij laplacian(rho) / 2)`. The
Landau-de Gennes energy and the current of the density are interactions of the model.
Args:
label:
The label of the site or the field.
visual:
How the site is drawn, or `None` for `default_visual_`.
rho:
The initial density.
Q:
The initial alignment tensor, its independent components `Q_xx` and `Q_xy`.
beta:
The coupling of the tensor to the curvature of the density.
Attributes:
default_visual_:
The visual of the element when it sets none.
Examples:
>>> from skcomplex.models.physics.active_matter import NematicDensitySiteElement
>>> NematicDensitySiteElement('site').beta
0.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
rho: Annotated[npt.ArrayLike, StateSpec()] = 1.0,
Q: Annotated[npt.ArrayLike, StateSpec(tensor='traceless')] = 0.0,
beta: Annotated[float, ParameterSpec()] = 0.5,
) -> None:
super().__init__(label, visual=visual)
self.rho = rho
self.Q = Q
self.beta = beta
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the rate of the tensor from the curvature of the density.
Args:
namespace:
The namespace of the symbols the laws can use.
Returns:
laws:
The rate of `Q`.
"""
s = namespace
curvature = symmetric_part(vector_gradient(gradient(s.rho)))
traceless = curvature - curvature.trace() / 2 * identity_matrix(SR, 2)
return [RateLaw(s.Q, s.beta * traceless)]
class LandauDeGennesEnergyInteraction(GradientInteraction):
"""Gradient flow of a symmetric traceless tensor `Q` in the Landau-de Gennes free energy.
The free energy is `F = integral a tr(Q^2) / 2 + c tr(Q^2)^2 / 4 + K |grad Q|^2 / 2`. Its negative functional
derivative among symmetric traceless tensors is the molecular field `H = -a Q - c tr(Q^2) Q + K laplacian(Q)`, and
the tensor relaxes as `Q_t = H / gamma`. A uniform tensor `S (n n - I / 2)` in two dimensions with `a < 0` relaxes
to the scalar order `S = sqrt(-2 a / c)`.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field.
state:
The name of the state of the field that holds the tensor.
a:
The coefficient of `tr(Q^2) / 2`, negative in the ordered phase.
c:
The coefficient of `tr(Q^2)^2 / 4`.
K:
The elastic constant.
gamma:
The rotational viscosity, which divides the molecular field.
quadratic:
A function of the namespace that returns the coefficient of `tr(Q^2) / 2` instead of `a`, such as an
expression of a density, or `None` for `a`.
Examples:
>>> from sage.all__sagemath_symbolics import SR, matrix, var
>>> from skcomplex.models.physics.active_matter import LandauDeGennesEnergyInteraction
>>> q, a, c = var('q, a, c')
>>> tensor = matrix(SR, [[q, 0], [0, -q]])
>>> LandauDeGennesEnergyInteraction._molecular_field(tensor, a, c, 0)[0, 0].factor()
-(2*c*q^2 + a)*q
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'Q',
a: Annotated[float, ParameterSpec()] = -1.0,
c: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
K: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
quadratic: Callable[[SimpleNamespace], object] | None = None,
) -> None:
super().__init__(label, target=target, states=state, variational=True)
self.state = state
self.a = a
self.c = c
self.K = K
self.gamma = gamma
self.quadratic = quadratic
@staticmethod
def _density(tensor: Matrix, a: object, c: object, K: object) -> Expression:
"""Return the density of the free energy, with the gradient term by parts."""
square = (tensor * tensor).trace()
return a * square / 2 + c * square**2 / 4 - K * (tensor * tensor.apply_map(laplacian)).trace() / 2
@classmethod
def _molecular_field(cls, tensor: Matrix, a: object, c: object, K: object) -> Matrix:
"""Return the molecular field of a tensor, the negative functional derivative of the free energy."""
energy = cls('energy', target='field')
return -energy._get_matrix_derivative(cls._density(tensor, a, c, K), tensor)
def _get_potential(self, namespace: SimpleNamespace) -> Expression:
"""Return the density of the free energy over the rotational viscosity, whose gradient is `H / gamma`.
Args:
namespace:
The namespace of the symbols the density can use.
Returns:
density:
The density that the tensor moves down.
"""
quadratic = namespace.a if self.quadratic is None else self.quadratic(namespace)
tensor = getattr(namespace.target, self.state)
return self._density(tensor, quadratic, namespace.c, namespace.K) / namespace.gamma
class ConservationInteraction(Interaction):
"""Conservation law `du/dt + divergence(flux(ns)) = source(ns)` of a state `u` of a field.
The flux is one expression for a field of one dimension, or one expression per dimension, of the states of the
field, and the interaction lowers it to a finite-volume flux law with a slope limiter and a numerical flux, or,
without a numerical flux, to the rate `-divergence(flux)` that the field discretizes by its own operators, as a
spectral or finite-difference field does. Without a source the total of the state on a periodic field stays
constant. The conservation laws of several states of one field form one system, whose numerical flux reads the
wave speeds of all of them, and a flux may read parameters that vary across the cells. A flux or a source function
takes the namespace of the interaction, with the namespace `target` of the field, the time `t`, and the symbols of
the parameters a subclass declares.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field whose state is conserved.
state:
The name of the conserved state of the field.
flux:
A function of the namespace that returns the flux, or one flux per dimension.
source:
A function of the namespace that returns the source, or `None` for none.
limiter:
The slope limiter of the reconstruction, `'minmod'`, `'van-leer'`, or `None` for first order.
riemann:
The numerical flux of a finite-volume field, `'rusanov'` or `'godunov'`, the exact solver for one state, or
`None` for the divergence of the flux on a field of another discretization.
sources:
The paths of other elements the flux 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 ConservationInteraction
>>> 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': 1.0 + np.sin(2 * np.pi * np.arange(50) / 50)}
>>> volumes = 'finite-volume'
>>> field = FieldElement('wave', element=CellElement, space=space, values=values, discretization=volumes)
>>> burgers = ConservationInteraction('burgers', target='wave', state='u', flux=lambda s: s.target.u**2 / 2)
>>> system = System(elements=[field], interactions=[burgers])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=0.5).trajectory_
>>> round(float(trajectory.values('wave.u')[-1].sum()), 6)
50.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'u',
flux: 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)
self.target = target
self.state = state
self.flux = flux
self.source = source
self.limiter = limiter
self.riemann = riemann
self.sources = sources
def _check(self) -> None:
"""Raise an error for a state that is not a name, a function 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 isinstance(self.state, str) or not self.state.isidentifier():
error_msg = f'{type(self).__name__} {self.label!r} conserves the state {self.state!r}, which is no name.'
raise ValidationError(error_msg)
for name in ('flux', 'source'):
value = getattr(self, name)
if value is not None and not callable(value):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
if self.limiter not in LIMITERS:
error_msg = f'Unknown limiter {self.limiter!r} of {self.label!r}. Use minmod, van-leer, or None.'
raise ValidationError(error_msg)
if self.riemann not in RIEMANN_SOLVERS:
error_msg = f'Unknown numerical flux {self.riemann!r} of {self.label!r}. Use rusanov, godunov, or None.'
raise ValidationError(error_msg)
def _get_flux(self, namespace: SimpleNamespace) -> object:
"""Return the flux, or `None` for none, the hook a subclass overrides to derive it from its arguments."""
return None if self.flux is None else self.flux(namespace)
def _get_transport_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the laws of the divergence of the flux, the hook a subclass overrides to lower its own flux."""
flux = self._get_flux(namespace)
if flux is None:
return []
state = getattr(self._get_member(namespace, 'target'), self.state)
if self.riemann is None:
components = list(flux) if isinstance(flux, FreeModuleElement | list | tuple) else [flux]
return [RateLaw(state, -divergence(components))]
return [FluxLaw(state, flux, limiter=self.limiter, riemann=self.riemann)]
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the laws of the divergence of the flux and the rate of the source.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespace of the field.
Returns:
laws:
The flux law, and the rate law of the source when there is one.
"""
laws = self._get_transport_laws(namespace)
if self.source is not None:
laws.append(RateLaw(getattr(self._get_member(namespace, 'target'), self.state), self.source(namespace)))
return laws
Equations
nematic.rho[i]: diff(rho_nematic(t), t) == -partial(-0.5*partial(-Q_nematic_xx(t), 1) - 0.5*partial(Q_nematic_xy(t), 0) - 0.5*partial(rho_nematic(t), 1), 1) - partial(-0.5*partial(Q_nematic_xx(t), 0) - 0.5*partial(Q_nematic_xy(t), 1) - 0.5*partial(rho_nematic(t), 0), 0)
nematic.Q_xx[i]: diff(Q_nematic_xx(t), t) == -1/2*beta_nematic*(partial(partial(rho_nematic(t), 1), 1) - partial(partial(rho_nematic(t), 0), 0)) + 1/2*K_energy*laplacian(Q_nematic_xx(t))/gamma_energy - 1/2*(4*(Q_nematic_xx(t)^2 + Q_nematic_xy(t)^2)*c_energy*Q_nematic_xx(t) + 2*a_energy*(-1.0*rho_nematic(t) + 1)*Q_nematic_xx(t) - K_energy*laplacian(Q_nematic_xx(t)))/gamma_energy
nematic.Q_xy[i]: diff(Q_nematic_xy(t), t) == 1/2*beta_nematic*(partial(partial(rho_nematic(t), 1), 0) + partial(partial(rho_nematic(t), 0), 1)) + 1/2*K_energy*laplacian(Q_nematic_xy(t))/gamma_energy - 1/2*(4*(Q_nematic_xx(t)^2 + Q_nematic_xy(t)^2)*c_energy*Q_nematic_xy(t) + 2*a_energy*(-1.0*rho_nematic(t) + 1)*Q_nematic_xy(t) - K_energy*laplacian(Q_nematic_xy(t)))/gamma_energy