HPP lattice gas
The HPP lattice gas on a periodic square lattice, whose particles move along four directions and collide.
- Name:
hpp_lattice_gas - Field:
physics/fluids - Formalism: lattice dynamics
- References: Hardy, de Pazzis, and Pomeau 1976; Lindgren 2024
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('hpp_lattice_gas')
entry = describe('hpp_lattice_gas')
trajectory = create_simulation('hpp_lattice_gas').prepare().run().trajectory_
Factory
@register_system(
field='physics/fluids',
tags=('lattice', 'lattice gas', 'deterministic', 'conservation', 'reversible'),
references=('Hardy, de Pazzis, and Pomeau 1976', 'Lindgren 2024'),
defaults={'until': 200, 'dt': 1},
)
def hpp_lattice_gas(
shape: Sequence[int] = (64, 64),
*,
density: float = 0.2,
bump: float = 0.6,
radius: float | None = None,
full: bool = False,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return the HPP lattice gas on a periodic square lattice, whose particles move along four directions and collide.
Every site holds up to four particles, one per direction. Every step the particles first collide, where a site with
exactly two particles of opposite directions turns them by a quarter turn and every other site keeps its particles,
and then every particle moves one site along its direction. The rules are deterministic and reversible, and they
conserve the number of particles and the momentum, and also the momentum along every row and every column, a
spurious invariant that keeps the HPP gas from isotropic hydrodynamics. Every channel is occupied at random with
probability `density`, or `bump` inside a disc at the centre. With `full`, as in Problem 4.3 of Lindgren, the
particles start only in full cells, which hold all four directions, placed at random with probability `density`.
Args:
shape:
The number of rows and columns.
density:
The occupation probability of a channel outside the disc, or of a full cell with `full`.
bump:
The occupation probability of a channel inside the disc, unused with `full`.
radius:
The radius of the disc in lattice spacings, or `None` for a sixth of the smaller side.
full:
Whether the particles start in full cells only.
random_state:
The seed or generator of the initial occupations.
Returns:
system:
The system, with the channels `gas.n0` to `gas.n3` and the particle numbers `gas.density`.
Examples:
>>> from skcomplex.systems.physics.fluids import hpp_lattice_gas
>>> hpp_lattice_gas((8, 8)).elements[0].prepare().size_
64
"""
space = LatticeSpace(tuple(shape), boundary='periodic').prepare()
rng = np.random.default_rng(random_state)
if full:
occupied = np.repeat((rng.random((space.size_, 1)) < density).astype(float), SQUARE_DIRECTIONS, axis=1)
else:
rows, columns = np.indices(tuple(shape))
positions = np.column_stack([rows.ravel(), columns.ravel()]).astype(float)
size = min(shape) / 6 if radius is None else radius
inside = np.linalg.norm(positions - positions.mean(axis=0), axis=1) < size
probability = np.where(inside, bump, density)[:, None]
occupied = (rng.random((space.size_, SQUARE_DIRECTIONS)) < probability).astype(float)
channels = {name: occupied[:, direction] for direction, name in enumerate(SQUARE_CHANNELS)}
gas = LatticeElement(
'gas',
element=HPPSiteElement,
space=space,
values={'density': occupied.sum(axis=1), **channels},
visual=DENSITY_VISUAL,
)
return System(
elements=[gas],
interactions=[
LocalRuleInteraction(
'collision',
target='gas',
state=SQUARE_CHANNELS,
rule=_collide_square,
neighborhood=SITE,
),
LocalRuleInteraction(
'propagation',
target='gas',
state=SQUARE_CHANNELS,
rule=_propagate_square,
neighborhood=ARRIVALS,
),
LocalRuleInteraction(
'count',
target='gas',
state='density',
rule=_count,
neighborhood=SITE,
reads=SQUARE_CHANNELS,
),
],
)
Components
Elements:
HPPSiteElement, fromElement
Interactions:
LocalRuleInteraction
class HPPSiteElement(Element):
"""Site of the HPP lattice gas, with one boolean channel per direction of the square lattice.
The local rules of the gas turn every head-on pair of particles at a site by a quarter turn, and then move every
particle to the neighbouring site along its direction. The state `density` holds the number of particles at the
site.
Args:
label:
The label of the site or the lattice group.
visual:
How the element is drawn, or `None` for `default_visual_`.
n0:
The initial occupation of the channel along the columns, towards higher column indices.
n1:
The initial occupation of the channel along the rows, towards higher row indices.
n2:
The initial occupation of the channel along the columns, towards lower column indices.
n3:
The initial occupation of the channel along the rows, towards lower row indices.
density:
The initial number of particles at the site.
Attributes:
default_visual_:
The visual of the site when it sets none.
Examples:
>>> from skcomplex.systems.physics.fluids import HPPSiteElement
>>> HPPSiteElement('site').n0
0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
n0: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
n1: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
n2: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
n3: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
density: Annotated[float, StateSpec(domain=NN)] = 0,
) -> None:
super().__init__(label, visual=visual)
self.n0 = n0
self.n1 = n1
self.n2 = n2
self.n3 = n3
self.density = density
class LocalRuleInteraction(Interaction):
"""Local rule of a cellular automaton, which sets the next value of every cell from its ordered neighbourhood.
The rule is a table of the next value of every configuration of the neighbourhood, `TableSpec`, a rule of the sum
of the neighbourhood, `TotalisticSpec`, a rule of the cell and the sum of its neighbours, `OuterTotalisticSpec`, an
object whose `table()` returns a `TableSpec`, or a function. A function receives one array per neighbour of every
read state, in neighbourhood order, with one value per cell, the uniform array `random` when the rule is stochastic,
and the parameters of the bound group it names, and returns the next values of the states, one array per state. A
state of several components per cell, such as a truncated power series of `PowerSeriesSpec`, reaches a function as
one row of components per cell. A `StateTranslationNeighborhood` reads a group space at the product of every site
and the group element that the first read state holds there.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the group on the cells of a lattice, the nodes of a network, or a Cayley graph.
state:
The name of the state the rule sets, or a tuple of the names of the states it sets together.
rule:
The table, the totalistic rule, or the function.
neighborhood:
The ordered neighbourhood of every cell, such as `VonNeumannNeighborhood(center=True)`.
reads:
The names of the states whose neighbourhoods a function reads, or `None` for the states it sets.
schedule:
`'synchronous'`, `'random-sequential'`, `'random-site'`, `'checkerboard'`, or `'colored'`.
stochastic:
Whether a function receives the uniform random array `random`.
Raises:
ValidationError:
The state or the read states are not names, or the neighbourhood is not a neighbourhood.
Examples:
>>> from skcomplex.core.components.interactions import LocalRuleInteraction
>>> from skcomplex.core.components.spaces import VonNeumannNeighborhood
>>> from skcomplex.core.specs import TotalisticSpec
>>> parity = TotalisticSpec({1: 1, 3: 1})
>>> around = VonNeumannNeighborhood(n_dim=1, center=True)
>>> rule = LocalRuleInteraction('parity', target='cells', state='cell', rule=parity, neighborhood=around)
>>> rule.schedule
'synchronous'
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (LocalRuleLaw,)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('state', 'reads'))],
state: str | tuple[str, ...],
rule: object,
neighborhood: BaseNeighborhood,
reads: Sequence[str] | None = None,
schedule: str = 'synchronous',
stochastic: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.rule = rule
self.neighborhood = neighborhood
self.reads = reads
self.schedule = schedule
self.stochastic = stochastic
def _get_names(self) -> tuple[str, ...]:
"""Return the names of the states the rule sets."""
return (self.state,) if isinstance(self.state, str) else tuple(self.state)
def _check(self) -> None:
"""Raise an error for states that are not names, or a neighbourhood that is not one.
Raises:
ValidationError:
The state or the read states are not names, or the neighbourhood is not a neighbourhood.
"""
super()._check()
names = [*self._get_names(), *(self.reads or ())]
if not names or not all(isinstance(name, str) and name.isidentifier() for name in names):
error_msg = f'The local rule {self.label!r} sets {self.state!r} and reads {self.reads!r}. Give state names.'
raise ValidationError(error_msg)
if not isinstance(self.neighborhood, BaseNeighborhood):
error_msg = f'The neighbourhood of the local rule {self.label!r} is {self.neighborhood!r}. Give one.'
raise ValidationError(error_msg)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the local rule law of the states of the target.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target`.
Returns:
laws:
The local rule law.
"""
targets = tuple(getattr(self._get_member(namespace, 'target'), name) for name in self._get_names())
reads = (
None
if self.reads is None
else tuple(getattr(self._get_member(namespace, 'target'), name) for name in self.reads)
)
return [
LocalRuleLaw(
targets[0] if isinstance(self.state, str) else targets,
self.neighborhood,
self.rule,
reads=reads,
schedule=self.schedule,
stochastic=self.stochastic,
label=self.label,
),
]