Fixed-energy sandpile

The fixed-energy sandpile, a Bak-Tang-Wiesenfeld pile with no drive and no loss on a periodic lattice.

• Name: fixed_energy_sandpile
• Field: physics/criticality
• Formalism: lattice dynamics
• References: Dickman, Vespignani, and Zapperi 1998; Vespignani et al. 2000; Gros 2015

Select this system

from skcomplex.systems import create_simulation, create_system, describe

system = create_system('fixed_energy_sandpile')
entry = describe('fixed_energy_sandpile')
trajectory = create_simulation('fixed_energy_sandpile').prepare().run().trajectory_

Factory

@register_system(
    field='physics/criticality',
    tags=('lattice', 'absorbing states', 'phase transition', 'conservation', 'sandpile'),
    references=('Dickman, Vespignani, and Zapperi 1998', 'Vespignani et al. 2000', 'Gros 2015'),
    defaults={'until': 500, 'dt': 1},
    observables=lambda system: [ObservableSpec('activity', _active, group='pile')],
)
def fixed_energy_sandpile(
    shape: Sequence[int] = (32, 32),
    *,
    density: float = 2.3,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the fixed-energy sandpile, a Bak-Tang-Wiesenfeld pile with no drive and no loss on a periodic lattice.

    Every step, every site with at least four grains gives one grain to each of its four neighbours, all at once.
    The lattice is periodic, so the density of grains `density` is conserved, and it controls an absorbing phase
    transition. Below the critical density, about 2.125 on the square lattice, the topplings stop in an absorbing
    configuration, and above it the pile stays active forever. The density of active sites, the observable `activity`,
    is the order parameter. The grains start scattered uniformly at random over the sites.

    Args:
        shape:
            The number of rows and columns.

        density:
            The mean number of grains per site.

        random_state:
            The seed or generator of the initial grains.

    Returns:
        system:
            The system, with the heights `pile.height` and the observable `activity`, the density of active sites.

    Examples:
        >>> from skcomplex.systems.physics.criticality import fixed_energy_sandpile
        >>> fixed_energy_sandpile((8, 8), density=2.0).elements[0].prepare().size_
        64
    """
    size = int(np.prod(shape))
    rng = np.random.default_rng(random_state)
    grains = rng.integers(size, size=round(density * size))
    heights = np.bincount(grains, minlength=size).reshape(tuple(shape)).astype(float)
    space = LatticeSpace(tuple(shape), boundary='periodic')
    return System(
        elements=[LatticeElement('pile', element=SandpileCellElement, space=space, values={'height': heights})],
        interactions=[ToppleInteraction('topple', target='pile', threshold=THRESHOLD)],
    )

Components

Elements:

• SandpileCellElement, from Element

Interactions:

• ToppleInteraction, from LocalRuleInteraction

class SandpileCellElement(Element):
    """Site of a sandpile, which holds a number of grains.

    An `AvalancheInteraction` drops the grains and topples the sites that hold too many.

    Args:
        label:
            The label of the site or the lattice group.

        visual:
            How the site is drawn, or `None` for `default_visual_`.

        height:
            The initial number of grains.

    Attributes:
        default_visual_:
            The visual of the site when it sets none.

    Examples:
        >>> from skcomplex.systems.physics.common import SandpileCellElement
        >>> SandpileCellElement('site').height
        0
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        height: Annotated[float, StateSpec(domain=NN)] = 0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.height = height

class ToppleInteraction(LocalRuleInteraction):
    """Parallel toppling of a sandpile without drive, where every site at or above the threshold topples once a step.

    Every step, every site that holds at least `threshold` grains loses `threshold` grains, and every site gains one
    grain from each toppling neighbour, all at once. On a periodic lattice no grain is lost, so the number of grains is
    conserved, and the activity either dies out in an absorbing configuration or lasts forever.

    Args:
        label:
            The label of the interaction, unique among the parts of its system.

        target:
            The path of the lattice group of sites.

        threshold:
            The number of grains at which a site topples, the number of its neighbours.

    Examples:
        >>> from skcomplex.systems.physics.criticality import ToppleInteraction
        >>> ToppleInteraction('topple', target='pile').threshold
        4
    """

    def __init__(self, label: str, *, target: Annotated[str, RoleSpec()], threshold: int = 4) -> None:
        super().__init__(
            label,
            target=target,
            state='height',
            rule=self._topple,
            neighborhood=VonNeumannNeighborhood(center=True),
        )
        self.threshold = threshold

    def _topple(self, height: Array, *around: Array) -> Array:
        """Return the heights after every site at or above the threshold topples once."""
        gained = np.sum([neighbor >= self.threshold for neighbor in around], axis=0)
        return height - self.threshold * (height >= self.threshold) + gained