"""Provides actors that influence scalar fields.
.. autosummary::
:nosignatures:
~LocalReactionsActor
~ScalarPDEActor
~DiffusionActor
~ReactionDiffusionActor
~CollectionPDEActor
.. codeauthor:: David Zwicker <david.zwicker@ds.mpg.de>
"""
from __future__ import annotations
import inspect
import warnings
from collections.abc import Callable
from typing import Any
import numpy as np
from pde import ScalarField
from pde.backends.numba.utils import jit
from pde.pdes.base import PDEBase
from pde.tools.docstrings import get_text_block
from pde.tools.expressions import ScalarExpression
from ... import Parameter
from ...elements import FieldCollectionElement, FieldElementBase, ScalarFieldElement
from ..base import ActorBase
[docs]
class LocalReactionsActor(ActorBase):
"""Actor simulating a local chemical reactions in a field."""
parameters_default = [
Parameter(
"reaction_flux",
"0",
str,
"An expression for the reaction flux in the mean field. The expression may "
"depend on the concentration and time, which are denoted by the variables "
"`c` and `t` respectively.",
),
]
element_classes = (FieldElementBase,)
def __init__(self, parameters: dict[str, Any] | None = None):
"""
Args:
parameters (dict):
Parameters affecting the actor. Call
:meth:`~LocalReactionsActor.show_parameters` for details.
"""
super().__init__(parameters=parameters)
reaction_flux = self.parameters["reaction_flux"]
self._reaction = ScalarExpression(reaction_flux, signature=["c", "t"])
[docs]
def estimate_dt(self, elements: tuple[FieldElementBase]) -> float: # type: ignore
"""Get the optimal time step for the simulation of the actor.
Args:
elements (tuple of :class:`~emulsim.elements.fields.MeanfieldElement`):
The element affected by the actor
"""
s_max = np.abs(self._reaction(np.linspace(0, 1, 32), t=0)).max()
if s_max == 0:
return float("inf")
else:
# The maximum timestep for s(c) = dc/dt is based on the
# analysis of the differential equation based on a linearized
# expression for the flux as dc(t)/dt = k * c(t),
# which has the solution c(t) = c(0) * exp(k * t),
# where k is a growth or reaction rate.
# Using an explicit Euler stepping, we find that the relative
# error ε during a single time step of length Δt is given by
# ε ≈ 0.5 * (k * Δt)**2 to lowest order in Δt. If we want to limit
# ε ≤ 0.01, we thus have to choose Δt ≤ sqrt(2 * ε) / k = 0.1 / k.
return 0.1 / s_max # type: ignore
[docs]
def make_evolver_numba( # type: ignore
self, elements: tuple[FieldElementBase]
) -> Callable[[tuple[np.ndarray], float, float], None]:
"""Return a function evolve the field from time `t` to `t + dt`
Args:
elements (tuple of :class:`~emulsim.elements.fields.MeanfieldElement`):
The element affected by the actor
Returns:
callable: A function with signature (field_data, t: float,
dt: float), which evolves the field_data.
"""
reation_flux = self._reaction.get_function(backend="numba")
@jit
def evolver(fields_data: tuple[np.ndarray], t: float, dt: float) -> None:
"""Evolve the diffusion equation explicitly."""
(field_data,) = fields_data
field_data += dt * reation_flux(field_data, t)
return evolver # type: ignore
[docs]
def evolve(self, elements: tuple[FieldElementBase], t: float, dt: float): # type: ignore
"""Evolve the field from time `t` to `t + dt`
Args:
elements (tuple of :class:`~emulsim.elements.fields.MeanfieldElement`):
The element affected by the actor
t (float):
The current time point
dt (float):
The time step used to evolve the element
"""
(element,) = elements # extract single element
element.data[...] += dt * self._reaction(element.data, t)
[docs]
class MeanfieldActor(LocalReactionsActor):
def __init__(self, parameters: dict[str, Any] | None = None):
"""
Args:
parameters (dict):
Parameters affecting the actor. Call
:meth:`~LocalReactionsActor.show_parameters` for details.
"""
# MeanfieldActor was deprecated on 2023-05-23
warnings.warn("MeanfieldActor is now LocalReactionsActor", DeprecationWarning)
super().__init__(parameters=parameters)
[docs]
class ScalarPDEActor(ActorBase):
"""Actor evolving a field according to a PDE."""
element_classes = (ScalarFieldElement,)
def __init__(self, pde: PDEBase, parameters: dict[str, Any] | None = None):
"""Initialize the actor and its PDE.
Args:
pde (:class:`~pde.pdes.base.PDEBase`):
The partial differential equation describing the dynamics of the
scalar field.
parameters (dict):
Parameters affecting the actor. Call
:meth:`~ScalarPDEActor.show_parameters` for details.
"""
super().__init__(parameters=parameters)
if inspect.isclass(pde):
self._logger.warning("Got class `%s` instead of an instance", pde)
self.pde: PDEBase = pde()
else:
self.pde = pde
@property
def info(self) -> dict[str, Any]:
"""dict: information about the actor"""
result = super().info
result["pde"] = {"class": self.pde.__class__.__name__}
return result
[docs]
def copy(self) -> ScalarPDEActor:
"""Returns a copy the actor."""
return self.__class__(self.pde, self.parameters.copy())
[docs]
def make_evolver_numba( # type: ignore
self, elements: tuple[ScalarFieldElement]
) -> Callable[[tuple[np.ndarray], float, float], None]:
"""Return a function evolving the field from time `t` to `t + dt`
Args:
elements (tuple of :class:`~emulsim.elements.fields.ScalarFieldElement`):
The element affected by the actor
Returns:
callable: A function with signature (field_data, t: float,
dt: float), which evolves the field_data.
"""
(element,) = elements # extract single element
pde_rhs = self.pde.make_pde_rhs(element._field, backend="numba")
@jit
def evolver(fields_data: tuple[np.ndarray], t: float, dt: float) -> None:
"""Evolve the PDE explicitly."""
(field_data,) = fields_data
field_data += dt * pde_rhs(field_data, t)
return evolver # type: ignore
[docs]
def evolve(self, elements: tuple[ScalarFieldElement], t: float, dt: float): # type: ignore
"""Evolve the field from time `t` to `t + dt`
Args:
elements (tuple of :class:`~emulsim.elements.fields.ScalarFieldElement`):
The element affected by the actor
t (float):
The current time point
dt (float):
The time step used to evolve the element
"""
(element,) = elements # extract single element
rate = self.pde.evolution_rate(element._field, t)
element._field += dt * rate
[docs]
class DiffusionActor(ScalarPDEActor):
"""Actor evolving a field according to a simple diffusion equation."""
parameters_default = [
Parameter(
"diffusivity",
1,
float,
"Diffusivity in the field. This actor only supports constant "
"diffusivities. Diffusivities depending on local concentration are "
"supported by `ReactionDiffusionActor`.",
),
Parameter(
"boundary_conditions",
"auto_periodic_neumann",
object,
"Defines the boundary conditions on the field."
+ get_text_block("ARG_BOUNDARIES"),
),
]
def __init__(self, parameters: dict[str, Any] | None = None):
"""
Args:
parameters (dict):
Parameters affecting the actor. Call
:meth:`~DiffusionActor.show_parameters` for details.
"""
from pde import DiffusionPDE
# parse parameters
parameters = self._parse_parameters(parameters, include_deprecated=True)
# initialize diffusion equation
eq = DiffusionPDE(
diffusivity=parameters["diffusivity"],
bc=parameters["boundary_conditions"],
)
super().__init__(eq, parameters)
[docs]
def estimate_dt(self, elements: tuple[ScalarFieldElement]) -> float: # type: ignore
"""Get the optimal time step for the simulation of the actor.
Args:
elements (tuple of :class:`~emulsim.elements.fields.ScalarFieldElement`):
The element affected by the actor
Returns:
float: the time step
"""
(element,) = elements # extract single element
dx = float(element.grid.discretization.min())
return 0.1 * dx**2 / float(self.pde.diffusivity) # type: ignore
[docs]
class ReactionDiffusionActor(ScalarPDEActor):
"""Actor evolving a field according to a reaction-diffusion equation.
This class relies on the optional `phasesep` package, which needs to be
installed separately.
"""
parameters_default = [
Parameter(
"diffusivity",
"1",
str,
"Diffusivity in the field. This can be an expression depending on the "
"local concentration that is parsed by `sympy`. Alternatively, simple "
"numbers are also supported.",
),
Parameter(
"reaction_flux",
"0",
str,
"An expression for the reaction flux in the field, which can depend on the "
"concentration (denoted by `c` or `phi`), spatial coordinates (denoted by "
"`x[i]`, where `i` is the dimension), and time `t`.",
),
Parameter(
"boundary_conditions",
"auto_periodic_neumann",
object,
"Defines the boundary conditions on the field."
+ get_text_block("ARG_BOUNDARIES"),
),
Parameter(
"expression_constants",
{},
dict,
"A dictionary defining values of constants that can be used in "
"expressions, e.g., the parameter `reaction_flux`.",
),
]
def __init__(self, parameters: dict[str, Any] | None = None):
"""
Args:
parameters (dict):
Parameters affecting the actor. Call
:meth:`~ReactionDiffusionActor.show_parameters` for details
"""
from phasesep.pdes import ReactionDiffusionPDE # type: ignore
# parse parameters
parameters = self._parse_parameters(parameters, include_deprecated=True)
# initialize reaction-diffusion equation
pde_params = {
"diffusivity": parameters["diffusivity"],
"reaction_flux": parameters["reaction_flux"],
"bc": parameters["boundary_conditions"],
"expression_constants": parameters["expression_constants"],
}
super().__init__(ReactionDiffusionPDE(pde_params), parameters)
[docs]
def estimate_dt(self, elements: tuple[ScalarFieldElement]) -> float: # type: ignore
"""Get the optimal time step for the simulation of the actor.
Args:
elements (tuple of :class:`~emulsim.elements.fields.ScalarFieldElement`):
The element affected by the actor
Returns:
float: the time step
"""
(element,) = elements # extract single element
grid = element.grid
# estimate the time step based on the chemical reaction
if hasattr(self.pde, "_reaction"):
# PDE seems to be an instance of ReactionDiffusionPDE
test_field = ScalarField(grid)
s_max = 0
for c in np.linspace(0, 1, 32):
test_field.data = c
rates = self.pde.reaction_rate(test_field) # type: ignore
s_max = max(s_max, np.max(np.abs(rates.data)))
diffusivity = self.pde.diffusivity.value # type: ignore
else:
# PDE seems to be an instance of DiffusionPDE
s_max = 0
diffusivity = self.pde.diffusivity # type: ignore
if s_max == 0:
dt_reaction = float("inf")
else:
dt_reaction = 0.2 / s_max # maximal 2% error during time step
# This estimate is based on an analysis of the differential equation
# dc(t)/dt = k * c(t), which has the solution c(t) = c(0) * exp(k * t), so
# k is a growth or reaction rate. Using an explicit Euler stepping, we find
# that the relative error ε during a single time step of length Δt is given
# by ε ≈ 0.5 * (k * Δt)**2 to lowest order in Δt. If we want to limit
# ε ≤ 0.02, we thus have to choose Δt ≤ sqrt(2 * ε) / k = 0.2 / k. The
# expression is thus a conservative estimate using the maximal reaction
# rate k = max(s).
# estimate the time step required for diffusion
dx = grid.discretization.min()
dt_diffusion = 0.1 * dx**2 / diffusivity
return min(dt_reaction, dt_diffusion) # type: ignore
[docs]
class CollectionPDEActor(ActorBase):
"""Actor evolving a field collection according to a PDE."""
element_classes = (FieldCollectionElement,)
def __init__(self, pde: PDEBase, parameters: dict[str, Any] | None = None):
"""Initialize the actor and its PDE.
Args:
pde (:class:`~pde.pdes.base.PDEBase`):
The partial differential equation describing the dynamics of the
scalar field.
parameters (dict):
Parameters affecting the actor. Call
:meth:`~CollectionPDEActor.show_parameters` for details.
"""
super().__init__(parameters=parameters)
if inspect.isclass(pde):
self._logger.warning("Got class `%s` instead of an instance", pde)
self.pde: PDEBase = pde()
else:
self.pde = pde
@property
def info(self) -> dict[str, Any]:
"""dict: information about the actor"""
result = super().info
result["pde"] = {"class": self.pde.__class__.__name__}
return result
[docs]
def copy(self) -> CollectionPDEActor:
"""Returns a copy the actor."""
return self.__class__(self.pde, self.parameters.copy())
[docs]
def make_evolver_numba( # type: ignore
self, elements: tuple[FieldCollectionElement]
) -> Callable[[tuple[np.ndarray], float, float], None]:
"""Return a function evolving the field from time `t` to `t + dt`
Args:
elements (tuple of :class:`~emulsim.elements.fields.ScalarFieldElement`):
The element affected by the actor
Returns:
callable: A function with signature (field_data, t: float,
dt: float), which evolves the field_data.
"""
(element,) = elements # extract single element
pde_rhs = self.pde.make_pde_rhs(element.field, backend="numba")
@jit
def evolver(fields_data: tuple[np.ndarray], t: float, dt: float) -> None:
"""Evolve the PDE explicitly."""
(field_data,) = fields_data
field_data += dt * pde_rhs(field_data, t)
return evolver # type: ignore
[docs]
def evolve(self, elements: tuple[FieldCollectionElement], t: float, dt: float): # type: ignore
"""Evolve the field from time `t` to `t + dt`
Args:
elements (tuple of :class:`~emulsim.elements.fields.ScalarFieldElement`):
The element affected by the actor
t (float):
The current time point
dt (float):
The time step used to evolve the element
"""
(element,) = elements # extract single element
rate = self.pde.evolution_rate(element.field, t)
element.field.data += dt * rate.data