FHP lattice gas
The FHP-I lattice gas on a periodic triangular lattice, with a dense disc that spreads into the gas.
- Name:
fhp_lattice_gas - Field:
physics/fluids - Formalism: lattice dynamics
- References: Frisch, Hasslacher, and Pomeau 1986; Dilão 2023
Factory
@register_model(
field='physics/fluids',
tags=('lattice', 'lattice gas', 'stochastic', 'conservation'),
references=('Frisch, Hasslacher, and Pomeau 1986', 'Dilão 2023'),
defaults={'until': 200, 'dt': 1},
)
def fhp_lattice_gas(
shape: Sequence[int] = (64, 64),
*,
density: float = 0.2,
bump: float = 0.6,
radius: float | None = None,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return the FHP-I lattice gas on a periodic triangular lattice, with a dense disc that spreads into the gas.
Every channel of every site is occupied at random with probability `density`, or `bump` inside a disc at the centre
of the lattice, so the excess density spreads as a sound wave and then diffuses.
Args:
shape:
The number of rows and columns of the triangular lattice.
density:
The occupation probability of a channel outside the disc.
bump:
The occupation probability of a channel inside the disc.
radius:
The radius of the disc in lattice spacings, or `None` for a sixth of the smaller side.
random_state:
The seed or generator of the initial occupations.
Returns:
system:
The model, with the channels `gas.n0` to `gas.n5` and the particle numbers `gas.density`.
Examples:
>>> from skcomplex.models.physics.fluids import fhp_lattice_gas
>>> fhp_lattice_gas((8, 8)).elements[0].check().prepare().size_
64
"""
space = LatticeSpace(tuple(shape), geometry='triangular').check().prepare()
positions = space.positions_
centre = positions.mean(axis=0)
size = min(space.shape_) / 6 if radius is None else radius
inside = np.linalg.norm(positions - centre, axis=1) < size
probability = np.where(inside, bump, density)[:, None]
occupied = (np.random.default_rng(random_state).random((space.size_, DIRECTIONS)) < probability).astype(float)
channels = {name: occupied[:, direction] for direction, name in enumerate(CHANNELS)}
gas = LatticeElement(
'gas',
element=FHPSiteElement,
space=space,
values={'density': occupied.sum(axis=1), **channels},
visual=DENSITY_VISUAL,
)
return System(
elements=[gas],
interactions=[
LocalRuleInteraction(
'collision',
target='gas',
state=CHANNELS,
rule=_collide,
neighborhood=SITE,
stochastic=True,
),
LocalRuleInteraction('propagation', target='gas', state=CHANNELS, rule=_propagate, neighborhood=TRIANGULAR),
LocalRuleInteraction(
'count',
target='gas',
state='density',
rule=_count,
neighborhood=SITE,
reads=CHANNELS,
),
],
)
Components
class FHPSiteElement(Element):
"""Site of the FHP-I lattice gas, with one boolean channel per direction of the triangular lattice.
The local rules of the gas collide the particles at every site, head-on pairs turning by a sixth of a turn to either
side with equal probability and symmetric triples turning by a sixth of a turn, and then move every particle to the
neighbouring site along its direction. Collisions conserve the number of particles and their momentum at every
site. 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 at 0 degrees.
n1:
The initial occupation of the channel at 60 degrees.
n2:
The initial occupation of the channel at 120 degrees.
n3:
The initial occupation of the channel at 180 degrees.
n4:
The initial occupation of the channel at 240 degrees.
n5:
The initial occupation of the channel at 300 degrees.
density:
The initial number of particles at the site.
Examples:
>>> from skcomplex.models.physics.fluids import FHPSiteElement
>>> FHPSiteElement('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,
n4: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0,
n5: 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.n4 = n4
self.n5 = n5
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.
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.interactions import LocalRuleInteraction
>>> from skcomplex.core.spaces import VonNeumannNeighborhood
>>> from skcomplex.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.check().schedule
'synchronous'
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
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,
),
]