Solver Module (softpotato.solver)

The softpotato.solver subpackage provides numerical solvers for 1D multi-species chemical diffusion and reaction-diffusion problems across planar Cartesian and curvilinear (spherical, cylindrical) geometries.

Architecture and Spatial Operator Delegation

Solvers follow the Single Responsibility Principle: - Time Stepping & Linear Algebra: Managed exclusively by solver classes (ScipyIVPSolver, CrankNicolson, ImplicitFiniteDifference, ExplicitFiniteDifference). - Spatial Meshes & Discretization: DiffusionProblem natively accepts either structured Mesh1D instances (Uniform1DMesh, GeometricExpandingMesh, ExponentialMesh) from softpotato.grid or raw 1D numpy.ndarray coordinate arrays. - Spatial Operators & Coordinate Metrics: Discretization stencils (\(\nabla^2 c\)), boundary conditions, interfacial fluxes, and Faradaic current calculations on both uniform and non-uniform meshes are delegated directly to Geometry abstractions (PlanarGeometry, SphericalGeometry, CylindricalGeometry).

Solvers and Supported Options

Solver Class & Aliases	Supported Configuration Options	Description & Defaults
ExplicitFiniteDifference ('explicit', 'ftcs')	dt (float, optional)	Explicit Forward-Time Central-Space. If dt=None, automatically defaults to \(0.45 \times \text{CFL limit}\) (\(\Delta t \le \Delta x_{\min}^2 / 2D_{\max}\)).
ImplicitFiniteDifference ('implicit', 'btcs')	dt (float, optional)	Backward-Time Central-Space implicit scheme using banded tridiagonal matrix inversion. If dt=None, defaults to 100 uniform steps: \((t_{\text{end}} - t_{\text{start}}) / 100\).
CrankNicolson ('crank_nicolson', 'crank-nicolson')	dt (float, optional)	Second-order in time and space Crank–Nicolson implicit scheme using banded tridiagonal solve. If dt=None, defaults to 100 uniform steps: \((t_{\text{end}} - t_{\text{start}}) / 100\).
ScipyIVPSolver ('scipy_ivp', 'solve_ivp')	method (str, default 'Radau'), rtol (float, default 1e-6), atol (float, default 1e-8)	Method of Lines ODE integration using scipy.integrate.solve_ivp. Supports stiff ('Radau', 'BDF') and non-stiff ('RK45') integrators.

Examples

1. Planar 1D Diffusion

import numpy as np
from softpotato.solver import DiffusionProblem, DirichletBC, get_solver

# Linear domain x in [0, 1]
x = np.linspace(0.0, 1.0, 50)
problem = DiffusionProblem(
    grid=x,
    diffusivity={"c": 1e-5},
    boundary_conditions={"c": (DirichletBC(0.0), DirichletBC(1.0))},
    initial_conditions={"c": 1.0},
    geometry="planar",  # default
)

solver = get_solver("crank_nicolson", dt=0.01)
result = solver.solve(problem, t_span=(0.0, 1.0))
print(result["c"].shape)  # (N_t, N_x) concentration profile
print(result.fluxes["c"].shape)  # (N_t,) interfacial flux
print(result.currents["c"].shape)  # (N_t,) Faradaic current in Amperes

2. Spherical Microelectrode Step Response with Expanding Mesh

from softpotato.geometry import SphericalGeometry
from softpotato.grid import GeometricExpandingMesh
from softpotato.solver import DiffusionProblem, DirichletBC, get_solver

r0 = 25e-6  # 25 um radius
r_mesh = GeometricExpandingMesh.from_diffusion_layer(
    D=1e-9, t_max=2.0, x_min=r0, n_points=75, gamma=1.06
)
problem = DiffusionProblem(
    grid=r_mesh,
    diffusivity={"c": 1e-9},
    boundary_conditions={"c": (DirichletBC(0.0), DirichletBC(1.0))},
    initial_conditions={"c": 1.0},
    geometry=SphericalGeometry(r0=r0),
    n_electrons=1,
)

solver = get_solver("scipy_ivp")
result = solver.solve(problem, t_span=(0.0, 2.0))
# result.currents['c'] contains exact spherical Cottrell chronoamperometric transient

Module Contents

solver

Numerical solvers for chemical species diffusion and reaction-diffusion problems.

This module provides a unified framework for solving 1D multi-species diffusion equations:

\[ \frac{\partial c_i}{\partial t} = D_i \frac{\partial^2 c_i}{\partial x^2} + R_i(\mathbf{c}, x, t) \]

Available built-in solvers include:

• "scipy_ivp" / "solve_ivp": Method of Lines using scipy.integrate.solve_ivp (Radau/BDF/RK45).
• "crank_nicolson": 2nd-order unconditionally stable Crank-Nicolson finite difference.
• "implicit" / "btcs": Unconditionally stable Backward-Time Central-Space.
• "explicit" / "ftcs": Forward-Time Central-Space with CFL stability enforcement.

Examples:

Solving a 1D diffusion problem and comparing two solvers:

>>> import numpy as np
>>> from softpotato.solver import (
...     DiffusionProblem,
...     DirichletBC,
...     NeumannBC,
...     get_solver,
...     list_solvers,
... )
>>>
>>> list_solvers()
['btcs', 'crank-nicolson', 'crank_nicolson', 'explicit', 'ftcs', 'implicit', 'scipy_ivp', 'solve_ivp']
>>>
>>> grid = np.linspace(0.0, 1.0, 41)
>>> problem = DiffusionProblem(
...     grid=grid,
...     diffusivity={"c": 0.01},
...     boundary_conditions={"c": (DirichletBC(0.0), DirichletBC(0.0))},
...     initial_conditions={"c": lambda x: np.sin(np.pi * x)},
... )
>>>
>>> t_eval = np.linspace(0.0, 0.5, 6)
>>> solver_mol = get_solver("scipy_ivp")
>>> res_mol = solver_mol.solve(problem, t_span=(0.0, 0.5), t_eval=t_eval)
>>> res_mol.success
True
>>>
>>> solver_cn = get_solver("crank_nicolson", dt=0.01)
>>> res_cn = solver_cn.solve(problem, t_span=(0.0, 0.5), t_eval=t_eval)
>>> res_cn.success
True

BaseSolver

Bases: ABC

Abstract Base Class for all diffusion solvers.

Supports automatic registration of subclasses using name="...":

Examples:

Creating and registering a user-made solver:

>>> import numpy as np
>>> from softpotato.solver import BaseSolver, SolverResult, get_solver, list_solvers
>>>
>>> class CustomSolver(BaseSolver, name="custom_demo"):
...     def _run_solver(self, problem, t_span, t_eval, **kwargs):
...         t = np.array([t_span[0], t_span[1]]) if t_eval is None else t_eval
...         c = {sp: np.zeros((len(t), problem.n_points)) for sp in problem.species}
...         return SolverResult(t=t, x=problem.grid, concentrations=c)
>>>
>>> "custom_demo" in list_solvers()
True
>>> my_solver = get_solver("custom_demo")
>>> isinstance(my_solver, CustomSolver)
True

Source code in src/softpotato/solver/base.py

class BaseSolver(ABC, metaclass=_BaseSolverMeta):
    """Abstract Base Class for all diffusion solvers.

    Supports automatic registration of subclasses using `name="..."`:

    Examples
    --------
    Creating and registering a user-made solver:

    >>> import numpy as np
    >>> from softpotato.solver import BaseSolver, SolverResult, get_solver, list_solvers
    >>>
    >>> class CustomSolver(BaseSolver, name="custom_demo"):
    ...     def _run_solver(self, problem, t_span, t_eval, **kwargs):
    ...         t = np.array([t_span[0], t_span[1]]) if t_eval is None else t_eval
    ...         c = {sp: np.zeros((len(t), problem.n_points)) for sp in problem.species}
    ...         return SolverResult(t=t, x=problem.grid, concentrations=c)
    >>>
    >>> "custom_demo" in list_solvers()
    True
    >>> my_solver = get_solver("custom_demo")
    >>> isinstance(my_solver, CustomSolver)
    True
    """

    _registry: ClassVar[dict[str, type[BaseSolver]]] = {}

    def __init_subclass__(
        cls,
        name: str | None = None,
        solver_name: str | None = None,
        **kwargs: Any,
    ) -> None:
        super().__init_subclass__(**kwargs)
        resolved_name = name or solver_name
        if resolved_name:
            cls._registry[resolved_name.lower()] = cls

    def __init__(self, **options: Any) -> None:
        """Initialize solver with optional solver-specific configuration.

        Parameters
        ----------
        **options : Any
            Arbitrary solver configuration options stored in ``self.options``.
            Concrete solver subclasses extract their parameters (such as ``dt``
            or SciPy integration options ``method``, ``rtol``, ``atol``) from
            this dictionary.
        """
        self.options = options

    def solve(
        self,
        problem: DiffusionProblem,
        t_span: tuple[float, float],
        t_eval: np.ndarray | None = None,
        **kwargs: Any,
    ) -> SolverResult:
        """Solve the diffusion problem over the specified time span.

        Parameters
        ----------
        problem : DiffusionProblem
            Problem definition containing grid, species, diffusivity, and boundary conditions.
        t_span : tuple[float, float]
            Time integration interval (t_start, t_end).
        t_eval : np.ndarray, optional
            Specific time points to evaluate and store the solution. If None, default
            evaluation time points will be generated by the solver.
        **kwargs : Any
            Additional run-time options passed to the underlying solver engine,
            temporarily overriding any matching configuration options set in ``__init__``.

        Returns
        -------
        SolverResult
            Standardized container holding time, coordinates, and concentration arrays.
        """
        self.validate(problem, t_span, t_eval)
        result = self._run_solver(problem, t_span, t_eval, **kwargs)
        return self._post_process(problem, result)

    @abstractmethod
    def _run_solver(
        self,
        problem: DiffusionProblem,
        t_span: tuple[float, float],
        t_eval: np.ndarray | None,
        **kwargs: Any,
    ) -> SolverResult:
        """Core numerical calculation loop or integrator dispatch."""
        raise NotImplementedError

    def validate(
        self,
        problem: DiffusionProblem,
        t_span: tuple[float, float],
        t_eval: np.ndarray | None,
    ) -> None:
        """Sanity and physical validity checks before integration."""
        if not isinstance(problem, DiffusionProblem):
            raise TypeError(
                f"problem must be an instance of DiffusionProblem, got {type(problem)}"
            )
        if t_span[0] >= t_span[1]:
            raise ValueError(
                f"t_span[1] must be strictly greater than t_span[0], got {t_span}"
            )
        if t_eval is not None:
            t_eval = np.asarray(t_eval, dtype=float)
            if t_eval.ndim != 1:
                raise ValueError("t_eval must be a 1D array.")
            if len(t_eval) < 2:
                raise ValueError("t_eval must have at least 2 points.")
            if t_eval[0] < t_span[0] or t_eval[-1] > t_span[1]:
                raise ValueError(
                    f"t_eval [{t_eval[0]}, {t_eval[-1]}] must lie within t_span {t_span}."
                )
            if np.any(np.diff(t_eval) <= 0):
                raise ValueError(
                    "t_eval points must be strictly monotonically increasing."
                )

    def _post_process(
        self, problem: DiffusionProblem, result: SolverResult
    ) -> SolverResult:
        """Compute surface fluxes, Faradaic currents, and perform post-run integrity checks."""
        # Calculate surface fluxes and Faradaic currents using problem.geometry.
        # Sign convention (electrode-normal into the solution domain):
        # - Left boundary (r = r_0 / x = 0, inward normal n_left = +r / +x):
        #     J_left(t) = -D * (dc/dn)|_{left}
        # - Right boundary (r = r_max / x = L, inward normal n_right = -r / -x):
        #     J_right(t) = -D * (dc/dn)|_{right} = +D * (dc/dr)|_{right}
        # Faradaic current:
        #     I = n * F * A * J
        geom = (
            get_geometry(problem.geometry)
            if isinstance(problem.geometry, str)
            else problem.geometry
        )
        if result.mesh is None and problem.mesh is not None:
            result.mesh = problem.mesh
        for sp in problem.species:
            if sp not in result.concentrations:
                continue
            c = result.concentrations[sp]
            D = problem.diffusivity[sp]
            if sp not in result.fluxes:
                result.fluxes[sp] = np.asarray(
                    geom.surface_flux(c, problem.grid, diffusivity=D, side="left")
                )
            if sp not in result.right_fluxes:
                result.right_fluxes[sp] = np.asarray(
                    geom.surface_flux(c, problem.grid, diffusivity=D, side="right")
                )

            n_electrons_val = getattr(problem, "n_electrons", 1)
            n_e = (
                n_electrons_val.get(sp, 1)
                if isinstance(n_electrons_val, dict)
                else n_electrons_val
            )

            # Compute Faradaic currents based on problem.electrodes ('left', 'right', or 'both')
            if (
                problem.electrodes in ("left", "both")
                and sp not in result.left_currents
            ):
                result.left_currents[sp] = np.asarray(
                    geom.current(
                        result.fluxes[sp],
                        n_electrons=n_e,
                        side="left",
                        grid=problem.grid,
                    ),
                    dtype=float,
                )
            if (
                problem.electrodes in ("right", "both")
                and sp not in result.right_currents
            ):
                result.right_currents[sp] = np.asarray(
                    geom.current(
                        result.right_fluxes[sp],
                        n_electrons=n_e,
                        side="right",
                        grid=problem.grid,
                    ),
                    dtype=float,
                )

            if sp not in result.currents:
                if problem.electrodes == "left":
                    result.currents[sp] = result.left_currents[sp]
                elif problem.electrodes == "right":
                    result.currents[sp] = result.right_currents[sp]
                elif problem.electrodes == "both":
                    result.currents[sp] = (
                        result.left_currents[sp] + result.right_currents[sp]
                    )
        return result

__init__

__init__(**options: Any) -> None

Initialize solver with optional solver-specific configuration.

Parameters:

Name	Type	Description	Default
**options	Any	Arbitrary solver configuration options stored in self.options. Concrete solver subclasses extract their parameters (such as dt or SciPy integration options method, rtol, atol) from this dictionary.	{}

Source code in src/softpotato/solver/base.py

def __init__(self, **options: Any) -> None:
    """Initialize solver with optional solver-specific configuration.

    Parameters
    ----------
    **options : Any
        Arbitrary solver configuration options stored in ``self.options``.
        Concrete solver subclasses extract their parameters (such as ``dt``
        or SciPy integration options ``method``, ``rtol``, ``atol``) from
        this dictionary.
    """
    self.options = options

solve

solve(problem: DiffusionProblem, t_span: tuple[float, float], t_eval: ndarray | None = None, **kwargs: Any) -> SolverResult

Solve the diffusion problem over the specified time span.

Parameters:

Name	Type	Description	Default
problem	DiffusionProblem	Problem definition containing grid, species, diffusivity, and boundary conditions.	required
t_span	tuple[float, float]	Time integration interval (t_start, t_end).	required
t_eval	ndarray	Specific time points to evaluate and store the solution. If None, default evaluation time points will be generated by the solver.	None
**kwargs	Any	Additional run-time options passed to the underlying solver engine, temporarily overriding any matching configuration options set in __init__.	{}

Returns:

Type	Description
SolverResult	Standardized container holding time, coordinates, and concentration arrays.

Source code in src/softpotato/solver/base.py

def solve(
    self,
    problem: DiffusionProblem,
    t_span: tuple[float, float],
    t_eval: np.ndarray | None = None,
    **kwargs: Any,
) -> SolverResult:
    """Solve the diffusion problem over the specified time span.

    Parameters
    ----------
    problem : DiffusionProblem
        Problem definition containing grid, species, diffusivity, and boundary conditions.
    t_span : tuple[float, float]
        Time integration interval (t_start, t_end).
    t_eval : np.ndarray, optional
        Specific time points to evaluate and store the solution. If None, default
        evaluation time points will be generated by the solver.
    **kwargs : Any
        Additional run-time options passed to the underlying solver engine,
        temporarily overriding any matching configuration options set in ``__init__``.

    Returns
    -------
    SolverResult
        Standardized container holding time, coordinates, and concentration arrays.
    """
    self.validate(problem, t_span, t_eval)
    result = self._run_solver(problem, t_span, t_eval, **kwargs)
    return self._post_process(problem, result)

validate

validate(problem: DiffusionProblem, t_span: tuple[float, float], t_eval: ndarray | None) -> None

Sanity and physical validity checks before integration.

Source code in src/softpotato/solver/base.py

def validate(
    self,
    problem: DiffusionProblem,
    t_span: tuple[float, float],
    t_eval: np.ndarray | None,
) -> None:
    """Sanity and physical validity checks before integration."""
    if not isinstance(problem, DiffusionProblem):
        raise TypeError(
            f"problem must be an instance of DiffusionProblem, got {type(problem)}"
        )
    if t_span[0] >= t_span[1]:
        raise ValueError(
            f"t_span[1] must be strictly greater than t_span[0], got {t_span}"
        )
    if t_eval is not None:
        t_eval = np.asarray(t_eval, dtype=float)
        if t_eval.ndim != 1:
            raise ValueError("t_eval must be a 1D array.")
        if len(t_eval) < 2:
            raise ValueError("t_eval must have at least 2 points.")
        if t_eval[0] < t_span[0] or t_eval[-1] > t_span[1]:
            raise ValueError(
                f"t_eval [{t_eval[0]}, {t_eval[-1]}] must lie within t_span {t_span}."
            )
        if np.any(np.diff(t_eval) <= 0):
            raise ValueError(
                "t_eval points must be strictly monotonically increasing."
            )

BoundaryCondition dataclass

Specification of a 1D boundary condition.

Parameters:

Name	Type	Description	Default
btype	('dirichlet', 'neumann')	Type of boundary condition: 'dirichlet': fixed concentration, \(c(x_{\text{bound}}, t) = \text{value}(t)\). 'neumann': fixed concentration gradient, \(\left.\frac{\partial c}{\partial x}\right|_{x_{\text{bound}}} = \text{value}(t)\).	'dirichlet'
value	float or Callable[[float], float]	Fixed value or a time-dependent callable \(f(t)\) returning a float.	required

Examples:

Constant concentration boundary condition:

>>> from softpotato.solver import BoundaryCondition
>>> bc = BoundaryCondition("dirichlet", 1.0)
>>> bc.evaluate(t=5.0)
1.0

Time-dependent concentration gradient (Neumann flux):

>>> bc_dynamic = BoundaryCondition("neumann", lambda t: -0.1 * t)
>>> bc_dynamic.evaluate(t=2.0)
-0.2

Source code in src/softpotato/solver/base.py

@dataclass
class BoundaryCondition:
    r"""Specification of a 1D boundary condition.

    Parameters
    ----------
    btype : {'dirichlet', 'neumann'}
        Type of boundary condition:

        - ``'dirichlet'``: fixed concentration, :math:`c(x_{\text{bound}}, t) = \text{value}(t)`.
        - ``'neumann'``: fixed concentration gradient, :math:`\left.\frac{\partial c}{\partial x}\right|_{x_{\text{bound}}} = \text{value}(t)`.
    value : float or Callable[[float], float]
        Fixed value or a time-dependent callable :math:`f(t)` returning a float.

    Examples
    --------
    Constant concentration boundary condition:

    >>> from softpotato.solver import BoundaryCondition
    >>> bc = BoundaryCondition("dirichlet", 1.0)
    >>> bc.evaluate(t=5.0)
    1.0

    Time-dependent concentration gradient (Neumann flux):

    >>> bc_dynamic = BoundaryCondition("neumann", lambda t: -0.1 * t)
    >>> bc_dynamic.evaluate(t=2.0)
    -0.2
    """

    btype: Literal["dirichlet", "neumann"]
    value: float | Callable[[float], float]

    def evaluate(self, t: float) -> float:
        """Evaluate boundary value at time t."""
        if callable(self.value):
            return float(self.value(t))
        return float(self.value)

evaluate

evaluate(t: float) -> float

Evaluate boundary value at time t.

Source code in src/softpotato/solver/base.py

def evaluate(self, t: float) -> float:
    """Evaluate boundary value at time t."""
    if callable(self.value):
        return float(self.value(t))
    return float(self.value)

DiffusionProblem dataclass

Formulation of a 1D multi-species chemical diffusion problem.

Governing equation for each species \(i\):

\[ \frac{\partial c_i}{\partial t} = D_i \frac{\partial^2 c_i}{\partial x^2} + R_i(\mathbf{c}, x, t) \]

Parameters:

Name	Type	Description	Default
grid	ndarray or Mesh1D	1D spatial grid points \(x \in [x_{\min}, x_{\max}]\) or a Mesh1D mesh object (e.g., Uniform1DMesh, GeometricExpandingMesh, ExponentialMesh).	required
diffusivity	dict[str, float]	Mapping from chemical species name to its diffusion coefficient \(D_i > 0\).	required
boundary_conditions	dict[str, tuple[BoundaryCondition, BoundaryCondition]]	Mapping from species name to a tuple of (left_bc, right_bc).	required
initial_conditions	dict[str, float | ndarray | Callable[[ndarray], ndarray]]	Mapping from species name to initial concentration profile (scalar, array, or callable).	required
reactions	Callable[[dict[str, ndarray], ndarray, float], dict[str, ndarray]]	Reaction kinetics term returning {species: reaction_rate_array}. Defaults to None.	None
species	list[str]	List of chemical species names. Inferred from diffusivity if not explicitly provided.	list()
geometry	Geometry | str	Coordinate geometry specifying differential operators (\(\nabla^2\)), normal surface flux evaluations, and active electrode area. Can be a registered alias (e.g., 'planar', 'spherical', 'cylindrical') or a Geometry instance. Defaults to 'planar'.	'planar'
n_electrons	int | dict[str, int]	Number of electrons transferred per reaction (\(n > 0\)) for Faradaic current calculations. Can be an integer or a species-to-integer mapping.	1
electrodes	('left', 'right', 'both')	Domain boundary or boundaries configured as active electrodes for Faradaic current evaluation. Defaults to 'left'. If 'both', currents from both sides are evaluated on SolverResult.left_currents and SolverResult.right_currents, and summed on SolverResult.currents.	'left'
electrode_side	('left', 'right', 'both')	Alias for electrodes.	'left'

Examples:

Setting up a 1D problem for two species ('A' and 'B'):

>>> import numpy as np
>>> from softpotato.solver import DiffusionProblem, DirichletBC, NeumannBC
>>>
>>> grid = np.linspace(0.0, 1.0, 50)
>>> problem = DiffusionProblem(
...     grid=grid,
...     diffusivity={"A": 1e-4, "B": 2e-4},
...     boundary_conditions={
...         "A": (DirichletBC(0.0), DirichletBC(1.0)),
...         "B": (NeumannBC(0.0), DirichletBC(0.0)),
...     },
...     initial_conditions={"A": 1.0, "B": 0.0},
... )
>>> problem.species
['A', 'B']
>>> problem.n_points
50
>>> bool(np.isclose(problem.dx, 1.0 / 49.0))
True

Source code in src/softpotato/solver/base.py

@dataclass
class DiffusionProblem:
    r"""Formulation of a 1D multi-species chemical diffusion problem.

    Governing equation for each species :math:`i`:

    .. math::

        \frac{\partial c_i}{\partial t} = D_i \frac{\partial^2 c_i}{\partial x^2} + R_i(\mathbf{c}, x, t)

    Parameters
    ----------
    grid : np.ndarray or Mesh1D
        1D spatial grid points :math:`x \in [x_{\min}, x_{\max}]` or a
        :class:`~softpotato.grid.Mesh1D` mesh object (e.g.,
        :class:`~softpotato.grid.Uniform1DMesh`,
        :class:`~softpotato.grid.GeometricExpandingMesh`,
        :class:`~softpotato.grid.ExponentialMesh`).
    diffusivity : dict[str, float]
        Mapping from chemical species name to its diffusion coefficient :math:`D_i > 0`.
    boundary_conditions : dict[str, tuple[BoundaryCondition, BoundaryCondition]]
        Mapping from species name to a tuple of ``(left_bc, right_bc)``.
    initial_conditions : dict[str, float | np.ndarray | Callable[[np.ndarray], np.ndarray]]
        Mapping from species name to initial concentration profile (scalar, array, or callable).
    reactions : Callable[[dict[str, np.ndarray], np.ndarray, float], dict[str, np.ndarray]], optional
        Reaction kinetics term returning ``{species: reaction_rate_array}``. Defaults to None.
    species : list[str], optional
        List of chemical species names. Inferred from diffusivity if not explicitly provided.
    geometry : Geometry | str, default 'planar'
        Coordinate geometry specifying differential operators (:math:`\nabla^2`),
        normal surface flux evaluations, and active electrode area. Can be a
        registered alias (e.g., ``'planar'``, ``'spherical'``, ``'cylindrical'``)
        or a :class:`~softpotato.geometry.Geometry` instance. Defaults to ``'planar'``.
    n_electrons : int | dict[str, int], default 1
        Number of electrons transferred per reaction (:math:`n > 0`) for Faradaic
        current calculations. Can be an integer or a species-to-integer mapping.
    electrodes : {'left', 'right', 'both'}, default 'left'
        Domain boundary or boundaries configured as active electrodes for Faradaic
        current evaluation. Defaults to ``'left'``. If ``'both'``, currents from both
        sides are evaluated on ``SolverResult.left_currents`` and ``SolverResult.right_currents``,
        and summed on ``SolverResult.currents``.
    electrode_side : {'left', 'right', 'both'}, optional
        Alias for ``electrodes``.

    Examples
    --------
    Setting up a 1D problem for two species ('A' and 'B'):

    >>> import numpy as np
    >>> from softpotato.solver import DiffusionProblem, DirichletBC, NeumannBC
    >>>
    >>> grid = np.linspace(0.0, 1.0, 50)
    >>> problem = DiffusionProblem(
    ...     grid=grid,
    ...     diffusivity={"A": 1e-4, "B": 2e-4},
    ...     boundary_conditions={
    ...         "A": (DirichletBC(0.0), DirichletBC(1.0)),
    ...         "B": (NeumannBC(0.0), DirichletBC(0.0)),
    ...     },
    ...     initial_conditions={"A": 1.0, "B": 0.0},
    ... )
    >>> problem.species
    ['A', 'B']
    >>> problem.n_points
    50
    >>> bool(np.isclose(problem.dx, 1.0 / 49.0))
    True
    """

    grid: np.ndarray | Mesh1D
    diffusivity: dict[str, float]
    boundary_conditions: dict[str, tuple[BoundaryCondition, BoundaryCondition]]
    initial_conditions: dict[
        str, float | np.ndarray | Callable[[np.ndarray], np.ndarray]
    ]
    reactions: (
        Callable[[dict[str, np.ndarray], np.ndarray, float], dict[str, np.ndarray]]
        | None
    ) = None
    species: list[str] = field(default_factory=list)
    geometry: Geometry | str = "planar"
    n_electrons: int | dict[str, int] = 1
    electrodes: Literal["left", "right", "both"] = "left"
    electrode_side: Literal["left", "right", "both"] | None = None
    mesh: Mesh1D | None = field(init=False, default=None, repr=False)

    def __post_init__(self) -> None:
        if self.electrode_side is not None:
            self.electrodes = self.electrode_side
        if self.electrodes not in ("left", "right", "both"):
            raise ValueError(
                f"electrodes must be 'left', 'right', or 'both', got {self.electrodes!r}."
            )
        self.electrode_side = self.electrodes

        if isinstance(self.geometry, str):
            self.geometry = get_geometry(self.geometry)
        elif not isinstance(self.geometry, Geometry):
            raise TypeError(
                f"geometry must be a str or Geometry instance, got {type(self.geometry)}."
            )
        if isinstance(self.grid, Mesh1D):
            self.mesh = self.grid
            self.grid = self.geometry.validate_grid(self.grid.nodes)
        else:
            self.mesh = None
            self.grid = self.geometry.validate_grid(self.grid)

        if not self.species:
            self.species = list(self.diffusivity.keys())

        for sp in self.species:
            if sp not in self.diffusivity:
                raise KeyError(f"Species '{sp}' missing in diffusivity dict.")
            if self.diffusivity[sp] <= 0:
                raise ValueError(f"Diffusivity for species '{sp}' must be positive.")
            if sp not in self.boundary_conditions:
                raise KeyError(f"Species '{sp}' missing in boundary_conditions dict.")
            if sp not in self.initial_conditions:
                raise KeyError(f"Species '{sp}' missing in initial_conditions dict.")

    @property
    def n_points(self) -> int:
        """Number of spatial grid points."""
        return len(self.grid)

    @property
    def bounds(self) -> tuple[float, float]:
        """Spatial domain boundaries (x_min, x_max)."""
        return float(self.grid[0]), float(self.grid[-1])

    @property
    def dx_array(self) -> np.ndarray:
        """1D array of consecutive spatial interval spacings."""
        if self.mesh is not None:
            return self.mesh.dx
        return np.diff(self.grid)

    @property
    def dx_min(self) -> float:
        """Minimum spatial interval spacing across the grid."""
        return float(np.min(self.dx_array))

    @property
    def dx_max(self) -> float:
        """Maximum spatial interval spacing across the grid."""
        return float(np.max(self.dx_array))

    @property
    def dx(self) -> float:
        """Grid step size for uniform grids."""
        if not self.is_uniform:
            raise ValueError("Grid is non-uniform; uniform dx is not defined.")
        return float(self.dx_array[0])

    @property
    def is_uniform(self) -> bool:
        """Check whether the spatial grid is uniformly spaced."""
        if self.mesh is not None:
            return self.mesh.is_uniform
        diffs = np.diff(self.grid)
        return bool(np.allclose(diffs, diffs[0], rtol=1e-5, atol=0.0))

    def get_initial_profile(self, sp: str) -> np.ndarray:
        """Evaluate initial concentration profile for a given species."""
        grid_arr = np.asarray(self.grid, dtype=float)
        ic = self.initial_conditions[sp]
        if callable(ic):
            arr = np.asarray(ic(grid_arr), dtype=float)
        elif np.isscalar(ic):
            arr = np.full_like(grid_arr, float(cast(Any, ic)), dtype=float)
        else:
            arr = np.asarray(ic, dtype=float)

        if arr.shape != grid_arr.shape:
            raise ValueError(
                f"Initial condition for '{sp}' has shape {arr.shape}, expected {grid_arr.shape}"
            )
        return arr

n_points property

n_points: int

Number of spatial grid points.

bounds property

bounds: tuple[float, float]

Spatial domain boundaries (x_min, x_max).

dx_array property

dx_array: ndarray

1D array of consecutive spatial interval spacings.

dx_min property

dx_min: float

Minimum spatial interval spacing across the grid.

dx_max property

dx_max: float

Maximum spatial interval spacing across the grid.

dx property

dx: float

Grid step size for uniform grids.

is_uniform property

is_uniform: bool

Check whether the spatial grid is uniformly spaced.

get_initial_profile

get_initial_profile(sp: str) -> np.ndarray

Evaluate initial concentration profile for a given species.

Source code in src/softpotato/solver/base.py

def get_initial_profile(self, sp: str) -> np.ndarray:
    """Evaluate initial concentration profile for a given species."""
    grid_arr = np.asarray(self.grid, dtype=float)
    ic = self.initial_conditions[sp]
    if callable(ic):
        arr = np.asarray(ic(grid_arr), dtype=float)
    elif np.isscalar(ic):
        arr = np.full_like(grid_arr, float(cast(Any, ic)), dtype=float)
    else:
        arr = np.asarray(ic, dtype=float)

    if arr.shape != grid_arr.shape:
        raise ValueError(
            f"Initial condition for '{sp}' has shape {arr.shape}, expected {grid_arr.shape}"
        )
    return arr

SolverResult dataclass

Standardized output container for diffusion simulations.

Parameters:

Name	Type	Description	Default
t	ndarray	1D array of time steps of shape (N_t,).	required
x	ndarray	1D array of spatial coordinates of shape (N_x,).	required
concentrations	dict[str, ndarray]	Dictionary mapping species name to 2D concentration array of shape (N_t, N_x).	required
fluxes	dict[str, ndarray]	Dictionary mapping species name to surface flux at left boundary x=0 of shape (N_t,).	dict()
right_fluxes	dict[str, ndarray]	Dictionary mapping species name to surface flux at right boundary x=L of shape (N_t,) under the electrode-normal convention (pointing into the solution domain).	dict()
currents	dict[str, ndarray]	Dictionary mapping species name to Faradaic current at the active electrode(s) of shape (N_t,) in Amperes. If electrodes='both', contains the net cell current (\(I_{ ext{net}} = I_{ ext{left}} + I_{ ext{right}}\)).	dict()
left_currents	dict[str, ndarray]	Dictionary mapping species name to Faradaic current at the left boundary of shape (N_t,) in Amperes.	dict()
right_currents	dict[str, ndarray]	Dictionary mapping species name to Faradaic current at the right boundary of shape (N_t,) in Amperes.	dict()
success	bool	Flag indicating if the numerical solution completed successfully.	True
message	str	Descriptive status message.	"Success"
raw_output	Any	Underlying solver output (e.g., OdeResult for SciPy).	None

Examples:

Accessing simulation results:

>>> import numpy as np
>>> from softpotato.solver import SolverResult
>>>
>>> res = SolverResult(
...     t=np.linspace(0, 1, 5),
...     x=np.linspace(0, 1, 10),
...     concentrations={"A": np.ones((5, 10))},
...     fluxes={"A": np.zeros(5)},
...     right_fluxes={"A": np.zeros(5)},
... )
>>> res.species
['A']
>>> res["A"].shape
(5, 10)
>>> res.fluxes["A"].shape
(5,)
>>> res.left_fluxes["A"].shape
(5,)
>>> res.right_fluxes["A"].shape
(5,)
>>> res.left_currents
{}
>>> res.right_currents
{}

Source code in src/softpotato/solver/base.py

@dataclass
class SolverResult:
    """Standardized output container for diffusion simulations.

    Parameters
    ----------
    t : np.ndarray
        1D array of time steps of shape (N_t,).
    x : np.ndarray
        1D array of spatial coordinates of shape (N_x,).
    concentrations : dict[str, np.ndarray]
        Dictionary mapping species name to 2D concentration array of shape (N_t, N_x).
    fluxes : dict[str, np.ndarray], optional
        Dictionary mapping species name to surface flux at left boundary x=0 of shape (N_t,).
    right_fluxes : dict[str, np.ndarray], optional
        Dictionary mapping species name to surface flux at right boundary x=L of shape (N_t,)
        under the electrode-normal convention (pointing into the solution domain).
    currents : dict[str, np.ndarray], optional
        Dictionary mapping species name to Faradaic current at the active electrode(s)
        of shape (N_t,) in Amperes. If ``electrodes='both'``, contains the net cell current
        (:math:`I_{\text{net}} = I_{\text{left}} + I_{\text{right}}`).
    left_currents : dict[str, np.ndarray], optional
        Dictionary mapping species name to Faradaic current at the left boundary of shape (N_t,)
        in Amperes.
    right_currents : dict[str, np.ndarray], optional
        Dictionary mapping species name to Faradaic current at the right boundary of shape (N_t,)
        in Amperes.
    success : bool, default True
        Flag indicating if the numerical solution completed successfully.
    message : str, default "Success"
        Descriptive status message.
    raw_output : Any, optional
        Underlying solver output (e.g., OdeResult for SciPy).

    Examples
    --------
    Accessing simulation results:

    >>> import numpy as np
    >>> from softpotato.solver import SolverResult
    >>>
    >>> res = SolverResult(
    ...     t=np.linspace(0, 1, 5),
    ...     x=np.linspace(0, 1, 10),
    ...     concentrations={"A": np.ones((5, 10))},
    ...     fluxes={"A": np.zeros(5)},
    ...     right_fluxes={"A": np.zeros(5)},
    ... )
    >>> res.species
    ['A']
    >>> res["A"].shape
    (5, 10)
    >>> res.fluxes["A"].shape
    (5,)
    >>> res.left_fluxes["A"].shape
    (5,)
    >>> res.right_fluxes["A"].shape
    (5,)
    >>> res.left_currents
    {}
    >>> res.right_currents
    {}
    """

    t: np.ndarray
    x: np.ndarray
    concentrations: dict[str, np.ndarray]
    fluxes: dict[str, np.ndarray] = field(default_factory=dict)
    right_fluxes: dict[str, np.ndarray] = field(default_factory=dict)
    currents: dict[str, np.ndarray] = field(default_factory=dict)
    left_currents: dict[str, np.ndarray] = field(default_factory=dict)
    right_currents: dict[str, np.ndarray] = field(default_factory=dict)
    success: bool = True
    message: str = "Success"
    raw_output: Any = None
    mesh: Mesh1D | None = None

    def __getitem__(self, species_name: str) -> np.ndarray:
        """Shorthand to access concentrations for a species: result['A']."""
        return self.concentrations[species_name]

    @property
    def left_fluxes(self) -> dict[str, np.ndarray]:
        """Surface fluxes at the left boundary (alias for ``fluxes``)."""
        return self.fluxes

    @property
    def species(self) -> list[str]:
        """List of species present in the solution."""
        return list(self.concentrations.keys())

left_fluxes property

left_fluxes: dict[str, ndarray]

Surface fluxes at the left boundary (alias for fluxes).

species property

species: list[str]

List of species present in the solution.

__getitem__

__getitem__(species_name: str) -> np.ndarray

Shorthand to access concentrations for a species: result['A'].

Source code in src/softpotato/solver/base.py

def __getitem__(self, species_name: str) -> np.ndarray:
    """Shorthand to access concentrations for a species: result['A']."""
    return self.concentrations[species_name]

CrankNicolson

Bases: BaseSolver

Crank-Nicolson (\(\theta = 0.5\)) finite difference diffusion solver.

Second-order accurate in both time and space: \(\mathcal{O}(\Delta t^2 + \Delta x^2)\). Unconditionally stable. Solves the coupled linear system at each time step:

\[ \frac{c_i^{n+1} - c_i^n}{\Delta t} = \frac{D}{2} \left[ \frac{c_{i+1}^{n+1} - 2c_i^{n+1} + c_{i-1}^{n+1}}{\Delta x^2} + \frac{c_{i+1}^n - 2c_i^n + c_{i-1}^n}{\Delta x^2} \right] + R_i^{n+1/2} \]

Employs scipy.linalg.solve_banded to invert the tridiagonal system in \(\mathcal{O}(N)\) time at each step.

Parameters:

Name	Type	Description	Default
dt	float	Fixed time step size \(\Delta t\). If None (default), automatically selects 100 uniform steps over the integration interval: \(\Delta t = (t_{\text{end}} - t_{\text{start}}) / 100\).	None
**options	Any	Additional solver options stored in self.options.	{}

Examples:

>>> import numpy as np
>>> from softpotato.solver import CrankNicolson, DiffusionProblem, DirichletBC
>>>
>>> grid = np.linspace(0.0, 1.0, 41)
>>> problem = DiffusionProblem(
...     grid=grid,
...     diffusivity={"c": 0.02},
...     boundary_conditions={"c": (DirichletBC(0.0), DirichletBC(0.0))},
...     initial_conditions={"c": lambda x: np.sin(np.pi * x)},
... )
>>> # 2nd-order accurate time marching
>>> solver = CrankNicolson(dt=0.01)
>>> result = solver.solve(problem, t_span=(0.0, 0.2))
>>> result.success
True
>>> result["c"].shape[1]
41

Source code in src/softpotato/solver/crank_nicolson.py

class CrankNicolson(BaseSolver, name="crank_nicolson"):
    r"""Crank-Nicolson (:math:`\theta = 0.5`) finite difference diffusion solver.

    Second-order accurate in both time and space: :math:`\mathcal{O}(\Delta t^2 + \Delta x^2)`.
    Unconditionally stable. Solves the coupled linear system at each time step:

    .. math::

        \frac{c_i^{n+1} - c_i^n}{\Delta t} = \frac{D}{2} \left[ \frac{c_{i+1}^{n+1} - 2c_i^{n+1} + c_{i-1}^{n+1}}{\Delta x^2} + \frac{c_{i+1}^n - 2c_i^n + c_{i-1}^n}{\Delta x^2} \right] + R_i^{n+1/2}

    Employs ``scipy.linalg.solve_banded`` to invert the tridiagonal system in :math:`\mathcal{O}(N)` time at each step.

    Parameters
    ----------
    dt : float, optional
        Fixed time step size :math:`\Delta t`.
        If None (default), automatically selects 100 uniform steps over the integration
        interval: :math:`\Delta t = (t_{\text{end}} - t_{\text{start}}) / 100`.
    **options : Any
        Additional solver options stored in ``self.options``.

    Examples
    --------
    >>> import numpy as np
    >>> from softpotato.solver import CrankNicolson, DiffusionProblem, DirichletBC
    >>>
    >>> grid = np.linspace(0.0, 1.0, 41)
    >>> problem = DiffusionProblem(
    ...     grid=grid,
    ...     diffusivity={"c": 0.02},
    ...     boundary_conditions={"c": (DirichletBC(0.0), DirichletBC(0.0))},
    ...     initial_conditions={"c": lambda x: np.sin(np.pi * x)},
    ... )
    >>> # 2nd-order accurate time marching
    >>> solver = CrankNicolson(dt=0.01)
    >>> result = solver.solve(problem, t_span=(0.0, 0.2))
    >>> result.success
    True
    >>> result["c"].shape[1]
    41
    """

    def __init__(self, dt: float | None = None, **options: Any) -> None:
        """Initialize Crank-Nicolson solver."""
        super().__init__(dt=dt, **options)

    def _run_solver(
        self,
        problem: DiffusionProblem,
        t_span: tuple[float, float],
        t_eval: np.ndarray | None,
        **kwargs: Any,
    ) -> SolverResult:
        geom = problem.geometry
        assert isinstance(geom, Geometry)
        dx = problem.dx if problem.is_uniform else None
        t_start, t_end = t_span

        user_dt = (
            kwargs.get("dt")
            if kwargs.get("dt") is not None
            else self.options.get("dt", None)
        )
        if user_dt is not None:
            base_dt = float(user_dt)
        else:
            base_dt = (t_end - t_start) / 100.0

        if t_eval is None:
            n_steps = max(2, int(np.ceil((t_end - t_start) / base_dt)) + 1)
            t_eval = np.linspace(t_start, t_end, n_steps)
        else:
            t_eval = np.asarray(t_eval, dtype=float)

        n_times = len(t_eval)
        N = problem.n_points
        curr_c = {sp: problem.get_initial_profile(sp) for sp in problem.species}
        laplacian_ops = {
            sp: geom.laplacian_matrix(
                problem.grid,
                bc_left=problem.boundary_conditions[sp][0],
                bc_right=problem.boundary_conditions[sp][1],
            )
            for sp in problem.species
        }

        has_initial_discontinuity = any(
            (
                problem.boundary_conditions[sp][0].btype == "dirichlet"
                and not np.isclose(
                    curr_c[sp][0],
                    problem.boundary_conditions[sp][0].evaluate(t_start),
                    rtol=1e-6,
                    atol=1e-9,
                )
            )
            or (
                problem.boundary_conditions[sp][1].btype == "dirichlet"
                and not np.isclose(
                    curr_c[sp][-1],
                    problem.boundary_conditions[sp][1].evaluate(t_start),
                    rtol=1e-6,
                    atol=1e-9,
                )
            )
            for sp in problem.species
        )
        rannacher_steps_remaining = 2 if has_initial_discontinuity else 0

        # If t_eval starts after t_start, advance from t_start to t_eval[0]
        t_current = t_start
        if t_eval[0] > t_start:
            interval = t_eval[0] - t_start
            substeps = max(1, int(np.ceil(interval / base_dt)))
            dt = interval / substeps
            for _ in range(substeps):
                if rannacher_steps_remaining > 0:
                    for _ in range(2):
                        curr_c = self._step_crank_nicolson(
                            problem,
                            curr_c,
                            t_current,
                            dt / 2.0,
                            dx,
                            theta=1.0,
                            laplacian_ops=laplacian_ops,
                        )
                        t_current += dt / 2.0
                    rannacher_steps_remaining -= 1
                else:
                    curr_c = self._step_crank_nicolson(
                        problem,
                        curr_c,
                        t_current,
                        dt,
                        dx,
                        laplacian_ops=laplacian_ops,
                    )
                    t_current += dt

        records = {sp: np.zeros((n_times, N), dtype=float) for sp in problem.species}
        for sp in problem.species:
            records[sp][0, :] = curr_c[sp]

        # Time marching with substepping
        for k in range(n_times - 1):
            t_k = t_eval[k]
            t_next = t_eval[k + 1]
            interval = t_next - t_k
            substeps = max(1, int(np.ceil(interval / base_dt)))
            dt = interval / substeps

            t_current = t_k
            for _ in range(substeps):
                if rannacher_steps_remaining > 0:
                    for _ in range(2):
                        curr_c = self._step_crank_nicolson(
                            problem,
                            curr_c,
                            t_current,
                            dt / 2.0,
                            dx,
                            theta=1.0,
                            laplacian_ops=laplacian_ops,
                        )
                        t_current += dt / 2.0
                    rannacher_steps_remaining -= 1
                else:
                    curr_c = self._step_crank_nicolson(
                        problem,
                        curr_c,
                        t_current,
                        dt,
                        dx,
                        laplacian_ops=laplacian_ops,
                    )
                    t_current += dt

            for sp in problem.species:
                records[sp][k + 1, :] = curr_c[sp]

        return SolverResult(
            t=t_eval,
            x=np.asarray(problem.grid, dtype=float),
            concentrations=records,
            success=True,
            message="Crank-Nicolson solution completed successfully.",
        )

    def _step_crank_nicolson(
        self,
        problem: DiffusionProblem,
        c_dict: dict[str, np.ndarray],
        t: float,
        dt: float,
        dx: float | None = None,
        *,
        theta: float = 0.5,
        laplacian_ops: dict[str, tuple[np.ndarray, np.ndarray]] | None = None,
    ) -> dict[str, np.ndarray]:
        """Advance all species profiles by one Crank-Nicolson time step."""
        grid_arr = np.asarray(problem.grid, dtype=float)
        reactions = (
            problem.reactions(c_dict, grid_arr, t) if problem.reactions else None
        )
        geom = problem.geometry
        assert isinstance(geom, Geometry)
        new_c: dict[str, np.ndarray] = {}
        N = problem.n_points
        t_next = t + dt

        for sp in problem.species:
            c = c_dict[sp]
            D = problem.diffusivity[sp]
            bc_left, bc_right = problem.boundary_conditions[sp]

            # Assemble or reuse banded Laplacian L and boundary scaling from geometry
            if laplacian_ops is not None and sp in laplacian_ops:
                l_banded, bc_scale = laplacian_ops[sp]
            else:
                l_banded, bc_scale = geom.laplacian_matrix(
                    problem.grid, bc_left=bc_left, bc_right=bc_right
                )

            ab = np.zeros((3, N), dtype=float)
            rhs = np.zeros(N, dtype=float)

            # Interior rows (1 <= i <= N - 2)
            ab[1, 1:-1] = 1.0 - theta * dt * D * l_banded[1, 1:-1]
            ab[0, 2:] = -theta * dt * D * l_banded[0, 2:]
            ab[2, :-2] = -theta * dt * D * l_banded[2, :-2]

            l_dot_c = (
                l_banded[2, :-2] * c[:-2]
                + l_banded[1, 1:-1] * c[1:-1]
                + l_banded[0, 2:] * c[2:]
            )
            rhs[1:-1] = c[1:-1] + (1.0 - theta) * dt * D * l_dot_c
            if reactions and sp in reactions:
                rhs[1:-1] += dt * reactions[sp][1:-1]

            # Left boundary (i = 0)
            if bc_left.btype == "dirichlet":
                ab[1, 0] = 1.0
                ab[0, 1] = 0.0
                rhs[0] = bc_left.evaluate(t_next)
            else:  # Neumann
                ab[1, 0] = 1.0 - theta * dt * D * l_banded[1, 0]
                ab[0, 1] = -theta * dt * D * l_banded[0, 1]
                g_left_curr = bc_left.evaluate(t)
                g_left_next = bc_left.evaluate(t_next)
                g_left_eff = theta * g_left_next + (1.0 - theta) * g_left_curr
                l_dot_c_0 = l_banded[1, 0] * c[0] + l_banded[0, 1] * c[1]
                rhs[0] = (
                    c[0]
                    + (1.0 - theta) * dt * D * l_dot_c_0
                    + dt * D * bc_scale[0] * g_left_eff
                )
                if reactions and sp in reactions:
                    rhs[0] += dt * reactions[sp][0]

            # Right boundary (i = N - 1)
            if bc_right.btype == "dirichlet":
                ab[1, -1] = 1.0
                ab[2, -2] = 0.0
                rhs[-1] = bc_right.evaluate(t_next)
            else:  # Neumann
                ab[1, -1] = 1.0 - theta * dt * D * l_banded[1, -1]
                ab[2, -2] = -theta * dt * D * l_banded[2, -2]
                g_right_curr = bc_right.evaluate(t)
                g_right_next = bc_right.evaluate(t_next)
                g_right_eff = theta * g_right_next + (1.0 - theta) * g_right_curr
                l_dot_c_end = l_banded[1, -1] * c[-1] + l_banded[2, -2] * c[-2]
                rhs[-1] = (
                    c[-1]
                    + (1.0 - theta) * dt * D * l_dot_c_end
                    + dt * D * bc_scale[-1] * g_right_eff
                )
                if reactions and sp in reactions:
                    rhs[-1] += dt * reactions[sp][-1]

            new_c[sp] = solve_banded((1, 1), ab, rhs)

        return new_c

__init__

__init__(dt: float | None = None, **options: Any) -> None

Initialize Crank-Nicolson solver.

Source code in src/softpotato/solver/crank_nicolson.py

def __init__(self, dt: float | None = None, **options: Any) -> None:
    """Initialize Crank-Nicolson solver."""
    super().__init__(dt=dt, **options)

ExplicitFiniteDifference

Bases: BaseSolver

Explicit Forward-Time Central-Space (FTCS) finite difference diffusion solver.

Numerically solves 1D multi-species diffusion equations:

\[ c_i^{n+1} = c_i^n + \Delta t \left[ D \frac{c_{i+1}^n - 2c_i^n + c_{i-1}^n}{\Delta x^2} + R_i^n \right] \]

Enforces the Courant-Friedrichs-Lewy (CFL) stability criterion using the minimum spatial interval \(\Delta x_{\min} = \min_j \Delta x_j\):

\[ \Delta t \le \frac{(\Delta x_{\min})^2}{2 \max_i(D_i)} \]

Parameters:

Name	Type	Description	Default
dt	float	Fixed time step size \(\Delta t\). Must satisfy the Courant-Friedrichs-Lewy (CFL) stability condition: \[ \Delta t \le \frac{(\Delta x_{\min})^2}{2 \max_i(D_i)} \] If None (default), automatically selects a conservative stable step size equal to \(0.45 \times \text{CFL limit}\).	None
**options	Any	Additional solver options stored in self.options.	{}

Examples:

>>> import numpy as np
>>> from softpotato.solver import ExplicitFiniteDifference, DiffusionProblem, DirichletBC
>>>
>>> grid = np.linspace(0.0, 1.0, 21)
>>> problem = DiffusionProblem(
...     grid=grid,
...     diffusivity={"A": 0.01},
...     boundary_conditions={"A": (DirichletBC(0.0), DirichletBC(1.0))},
...     initial_conditions={"A": 0.0},
... )
>>> # Use automatic CFL-safe time stepping
>>> solver = ExplicitFiniteDifference()
>>> result = solver.solve(problem, t_span=(0.0, 0.1))
>>> result.success
True
>>> result["A"].shape[1]
21

Source code in src/softpotato/solver/explicit.py

class ExplicitFiniteDifference(BaseSolver, name="explicit"):
    r"""Explicit Forward-Time Central-Space (FTCS) finite difference diffusion solver.

    Numerically solves 1D multi-species diffusion equations:

    .. math::

        c_i^{n+1} = c_i^n + \Delta t \left[ D \frac{c_{i+1}^n - 2c_i^n + c_{i-1}^n}{\Delta x^2} + R_i^n \right]

    Enforces the Courant-Friedrichs-Lewy (CFL) stability criterion using the minimum
    spatial interval :math:`\Delta x_{\min} = \min_j \Delta x_j`:

    .. math::

        \Delta t \le \frac{(\Delta x_{\min})^2}{2 \max_i(D_i)}

    Parameters
    ----------
    dt : float, optional
        Fixed time step size :math:`\Delta t`. Must satisfy the Courant-Friedrichs-Lewy (CFL)
        stability condition:

        .. math::

            \Delta t \le \frac{(\Delta x_{\min})^2}{2 \max_i(D_i)}

        If None (default), automatically selects a conservative stable step size
        equal to :math:`0.45 \times \text{CFL limit}`.
    **options : Any
        Additional solver options stored in ``self.options``.

    Examples
    --------
    >>> import numpy as np
    >>> from softpotato.solver import ExplicitFiniteDifference, DiffusionProblem, DirichletBC
    >>>
    >>> grid = np.linspace(0.0, 1.0, 21)
    >>> problem = DiffusionProblem(
    ...     grid=grid,
    ...     diffusivity={"A": 0.01},
    ...     boundary_conditions={"A": (DirichletBC(0.0), DirichletBC(1.0))},
    ...     initial_conditions={"A": 0.0},
    ... )
    >>> # Use automatic CFL-safe time stepping
    >>> solver = ExplicitFiniteDifference()
    >>> result = solver.solve(problem, t_span=(0.0, 0.1))
    >>> result.success
    True
    >>> result["A"].shape[1]
    21
    """

    def __init__(self, dt: float | None = None, **options: Any) -> None:
        """Initialize Explicit Finite Difference (FTCS) solver."""
        super().__init__(dt=dt, **options)

    def _run_solver(
        self,
        problem: DiffusionProblem,
        t_span: tuple[float, float],
        t_eval: np.ndarray | None,
        **kwargs: Any,
    ) -> SolverResult:
        geom = problem.geometry
        assert isinstance(geom, Geometry)
        dx_min = problem.dx_min
        dx = problem.dx if problem.is_uniform else dx_min
        max_D = max(problem.diffusivity.values())
        cfl_limit = (dx_min**2) / (2.0 * max_D)

        user_dt = (
            kwargs.get("dt")
            if kwargs.get("dt") is not None
            else self.options.get("dt", None)
        )
        if user_dt is not None:
            dt_req = float(user_dt)
            if dt_req > cfl_limit:
                raise ValueError(
                    f"Specified time step dt={dt_req:.3e} exceeds the CFL stability limit "
                    f"dx_min²/(2*D_max)={cfl_limit:.3e}."
                )
            base_dt = dt_req
        else:
            # Default to a safe sub-critical time step (45% of stability limit)
            base_dt = 0.45 * cfl_limit

        t_start, t_end = t_span
        if t_eval is None:
            # Generate evaluation time grid
            n_steps = max(2, int(np.ceil((t_end - t_start) / base_dt)) + 1)
            t_eval = np.linspace(t_start, t_end, n_steps)
        else:
            t_eval = np.asarray(t_eval, dtype=float)

        # Initialize current state
        curr_c = {sp: problem.get_initial_profile(sp) for sp in problem.species}

        # If t_eval starts after t_start, advance from t_start to t_eval[0]
        t_current = t_start
        if t_eval[0] > t_start:
            interval = t_eval[0] - t_start
            substeps = max(1, int(np.ceil(interval / base_dt)))
            dt = interval / substeps
            for _ in range(substeps):
                curr_c = self._step_explicit(problem, curr_c, t_current, dt, dx)
                t_current += dt

        # Solution records: shape (n_times, n_x) per species
        n_times = len(t_eval)
        records = {
            sp: np.zeros((n_times, problem.n_points), dtype=float)
            for sp in problem.species
        }
        for sp in problem.species:
            records[sp][0, :] = curr_c[sp]

        # Time marching across t_eval intervals with adaptive substepping
        for k in range(n_times - 1):
            t_k = t_eval[k]
            t_next = t_eval[k + 1]
            interval = t_next - t_k
            substeps = max(1, int(np.ceil(interval / base_dt)))
            dt = interval / substeps

            t_current = t_k
            for _ in range(substeps):
                curr_c = self._step_explicit(problem, curr_c, t_current, dt, dx)
                t_current += dt

            for sp in problem.species:
                records[sp][k + 1, :] = curr_c[sp]

        return SolverResult(
            t=t_eval,
            x=np.asarray(problem.grid, dtype=float),
            concentrations=records,
            success=True,
            message="Explicit finite difference solution completed successfully.",
        )

    def _step_explicit(
        self,
        problem: DiffusionProblem,
        c_dict: dict[str, np.ndarray],
        t: float,
        dt: float,
        dx: float | None = None,
    ) -> dict[str, np.ndarray]:
        """Advance all species concentration profiles by one time step dt."""
        grid_arr = np.asarray(problem.grid, dtype=float)
        reactions = (
            problem.reactions(c_dict, grid_arr, t) if problem.reactions else None
        )
        geom = problem.geometry
        assert isinstance(geom, Geometry)
        new_c: dict[str, np.ndarray] = {}
        t_next = t + dt

        for sp in problem.species:
            c = c_dict[sp]
            D = problem.diffusivity[sp]
            bc_left, bc_right = problem.boundary_conditions[sp]

            lap = geom.laplacian(
                c, problem.grid, bc_left=bc_left, bc_right=bc_right, t=t
            )

            # Advance state: c^{n+1} = c^n + dt * (D * lap + R)
            sp_reaction = reactions[sp] if reactions and sp in reactions else 0.0
            next_sp_c = c + dt * (D * lap + sp_reaction)

            # Enforce Dirichlet boundaries at t_next
            if bc_left.btype == "dirichlet":
                next_sp_c[0] = bc_left.evaluate(t_next)
            if bc_right.btype == "dirichlet":
                next_sp_c[-1] = bc_right.evaluate(t_next)

            new_c[sp] = next_sp_c

        return new_c

__init__

__init__(dt: float | None = None, **options: Any) -> None

Initialize Explicit Finite Difference (FTCS) solver.

Source code in src/softpotato/solver/explicit.py

def __init__(self, dt: float | None = None, **options: Any) -> None:
    """Initialize Explicit Finite Difference (FTCS) solver."""
    super().__init__(dt=dt, **options)

ImplicitFiniteDifference

Bases: BaseSolver

Implicit Backward-Time Central-Space (BTCS) finite difference diffusion solver.

Numerically solves 1D multi-species diffusion equations:

\[ \frac{c_i^{n+1} - c_i^n}{\Delta t} = D \frac{c_{i+1}^{n+1} - 2c_i^{n+1} + c_{i-1}^{n+1}}{\Delta x^2} + R_i^{n+1} \]

Unconditionally stable for arbitrary time steps \(\Delta t\). Solves the resulting tridiagonal system using scipy.linalg.solve_banded in \(\mathcal{O}(N)\) time per step.

Parameters:

Name	Type	Description	Default
dt	float	Fixed time step size \(\Delta t\). If None (default), automatically selects 100 uniform steps over the integration interval: \(\Delta t = (t_{\text{end}} - t_{\text{start}}) / 100\).	None
**options	Any	Additional solver options stored in self.options.	{}

Examples:

>>> import numpy as np
>>> from softpotato.solver import ImplicitFiniteDifference, DiffusionProblem, DirichletBC
>>>
>>> grid = np.linspace(0.0, 1.0, 31)
>>> problem = DiffusionProblem(
...     grid=grid,
...     diffusivity={"A": 0.05},
...     boundary_conditions={"A": (DirichletBC(0.0), DirichletBC(1.0))},
...     initial_conditions={"A": 0.0},
... )
>>> # Implicit schemes remain stable even with large time steps
>>> solver = ImplicitFiniteDifference(dt=0.05)
>>> result = solver.solve(problem, t_span=(0.0, 0.5))
>>> result.success
True
>>> result["A"].shape[1]
31

Source code in src/softpotato/solver/implicit.py

class ImplicitFiniteDifference(BaseSolver, name="implicit"):
    r"""Implicit Backward-Time Central-Space (BTCS) finite difference diffusion solver.

    Numerically solves 1D multi-species diffusion equations:

    .. math::

        \frac{c_i^{n+1} - c_i^n}{\Delta t} = D \frac{c_{i+1}^{n+1} - 2c_i^{n+1} + c_{i-1}^{n+1}}{\Delta x^2} + R_i^{n+1}

    Unconditionally stable for arbitrary time steps :math:`\Delta t`. Solves the resulting
    tridiagonal system using ``scipy.linalg.solve_banded`` in :math:`\mathcal{O}(N)` time per step.

    Parameters
    ----------
    dt : float, optional
        Fixed time step size :math:`\Delta t`.
        If None (default), automatically selects 100 uniform steps over the integration
        interval: :math:`\Delta t = (t_{\text{end}} - t_{\text{start}}) / 100`.
    **options : Any
        Additional solver options stored in ``self.options``.

    Examples
    --------
    >>> import numpy as np
    >>> from softpotato.solver import ImplicitFiniteDifference, DiffusionProblem, DirichletBC
    >>>
    >>> grid = np.linspace(0.0, 1.0, 31)
    >>> problem = DiffusionProblem(
    ...     grid=grid,
    ...     diffusivity={"A": 0.05},
    ...     boundary_conditions={"A": (DirichletBC(0.0), DirichletBC(1.0))},
    ...     initial_conditions={"A": 0.0},
    ... )
    >>> # Implicit schemes remain stable even with large time steps
    >>> solver = ImplicitFiniteDifference(dt=0.05)
    >>> result = solver.solve(problem, t_span=(0.0, 0.5))
    >>> result.success
    True
    >>> result["A"].shape[1]
    31
    """

    def __init__(self, dt: float | None = None, **options: Any) -> None:
        """Initialize Implicit Finite Difference (BTCS) solver."""
        super().__init__(dt=dt, **options)

    def _run_solver(
        self,
        problem: DiffusionProblem,
        t_span: tuple[float, float],
        t_eval: np.ndarray | None,
        **kwargs: Any,
    ) -> SolverResult:
        geom = problem.geometry
        assert isinstance(geom, Geometry)
        dx = problem.dx if problem.is_uniform else None
        t_start, t_end = t_span

        user_dt = (
            kwargs.get("dt")
            if kwargs.get("dt") is not None
            else self.options.get("dt", None)
        )
        if user_dt is not None:
            base_dt = float(user_dt)
        else:
            # Default to dividing total duration into 100 steps
            base_dt = (t_end - t_start) / 100.0

        if t_eval is None:
            n_steps = max(2, int(np.ceil((t_end - t_start) / base_dt)) + 1)
            t_eval = np.linspace(t_start, t_end, n_steps)
        else:
            t_eval = np.asarray(t_eval, dtype=float)

        n_times = len(t_eval)
        N = problem.n_points
        curr_c = {sp: problem.get_initial_profile(sp) for sp in problem.species}
        laplacian_ops = {
            sp: geom.laplacian_matrix(
                problem.grid,
                bc_left=problem.boundary_conditions[sp][0],
                bc_right=problem.boundary_conditions[sp][1],
            )
            for sp in problem.species
        }

        # If t_eval starts after t_start, advance from t_start to t_eval[0]
        t_current = t_start
        if t_eval[0] > t_start:
            interval = t_eval[0] - t_start
            substeps = max(1, int(np.ceil(interval / base_dt)))
            dt = interval / substeps
            for _ in range(substeps):
                curr_c = self._step_implicit(
                    problem, curr_c, t_current, dt, dx, laplacian_ops=laplacian_ops
                )
                t_current += dt

        records = {sp: np.zeros((n_times, N), dtype=float) for sp in problem.species}
        for sp in problem.species:
            records[sp][0, :] = curr_c[sp]

        # Time marching with substepping
        for k in range(n_times - 1):
            t_k = t_eval[k]
            t_next = t_eval[k + 1]
            interval = t_next - t_k
            substeps = max(1, int(np.ceil(interval / base_dt)))
            dt = interval / substeps

            t_current = t_k
            for _ in range(substeps):
                curr_c = self._step_implicit(
                    problem, curr_c, t_current, dt, dx, laplacian_ops=laplacian_ops
                )
                t_current += dt

            for sp in problem.species:
                records[sp][k + 1, :] = curr_c[sp]

        return SolverResult(
            t=t_eval,
            x=np.asarray(problem.grid, dtype=float),
            concentrations=records,
            success=True,
            message="Implicit finite difference solution completed successfully.",
        )

    def _step_implicit(
        self,
        problem: DiffusionProblem,
        c_dict: dict[str, np.ndarray],
        t: float,
        dt: float,
        dx: float | None = None,
        *,
        laplacian_ops: dict[str, tuple[np.ndarray, np.ndarray]] | None = None,
    ) -> dict[str, np.ndarray]:
        """Advance all species profiles by one BTCS implicit time step."""
        grid_arr = np.asarray(problem.grid, dtype=float)
        reactions = (
            problem.reactions(c_dict, grid_arr, t) if problem.reactions else None
        )
        geom = problem.geometry
        assert isinstance(geom, Geometry)
        new_c: dict[str, np.ndarray] = {}
        N = problem.n_points
        t_next = t + dt
        theta = 1.0

        for sp in problem.species:
            c = c_dict[sp]
            D = problem.diffusivity[sp]
            bc_left, bc_right = problem.boundary_conditions[sp]

            # Assemble or reuse banded Laplacian L and boundary scaling from geometry
            if laplacian_ops is not None and sp in laplacian_ops:
                l_banded, bc_scale = laplacian_ops[sp]
            else:
                l_banded, bc_scale = geom.laplacian_matrix(
                    problem.grid, bc_left=bc_left, bc_right=bc_right
                )

            ab = np.zeros((3, N), dtype=float)
            rhs = np.copy(c)
            if reactions and sp in reactions:
                rhs += dt * reactions[sp]

            # Interior rows (1 <= i <= N - 2)
            ab[1, 1:-1] = 1.0 - theta * dt * D * l_banded[1, 1:-1]
            ab[0, 2:] = -theta * dt * D * l_banded[0, 2:]
            ab[2, :-2] = -theta * dt * D * l_banded[2, :-2]

            # Left boundary (row 0)
            if bc_left.btype == "dirichlet":
                ab[1, 0] = 1.0
                ab[0, 1] = 0.0
                rhs[0] = bc_left.evaluate(t_next)
            else:  # Neumann
                ab[1, 0] = 1.0 - theta * dt * D * l_banded[1, 0]
                ab[0, 1] = -theta * dt * D * l_banded[0, 1]
                g_left = bc_left.evaluate(t_next)
                rhs[0] += dt * D * bc_scale[0] * g_left

            # Right boundary (row N-1)
            if bc_right.btype == "dirichlet":
                ab[1, -1] = 1.0
                ab[2, -2] = 0.0
                rhs[-1] = bc_right.evaluate(t_next)
            else:  # Neumann
                ab[1, -1] = 1.0 - theta * dt * D * l_banded[1, -1]
                ab[2, -2] = -theta * dt * D * l_banded[2, -2]
                g_right = bc_right.evaluate(t_next)
                rhs[-1] += dt * D * bc_scale[-1] * g_right

            # Solve tridiagonal system in O(N)
            new_c[sp] = solve_banded((1, 1), ab, rhs)

        return new_c

__init__

__init__(dt: float | None = None, **options: Any) -> None

Initialize Implicit Finite Difference (BTCS) solver.

Source code in src/softpotato/solver/implicit.py

def __init__(self, dt: float | None = None, **options: Any) -> None:
    """Initialize Implicit Finite Difference (BTCS) solver."""
    super().__init__(dt=dt, **options)

ScipyIVPSolver

Bases: BaseSolver

Diffusion solver using Method of Lines (MOL) and scipy.integrate.solve_ivp.

Discretizes the spatial diffusion operator into a coupled system of ODEs:

\[ \frac{dc_i}{dt} = D_i \frac{\partial^2 c_i}{\partial x^2} + R_i(\mathbf{c}, x, t) \]

Integrates the resulting stiff or non-stiff system using SciPy's ODE suite. Defaults to the stiff L-stable 'Radau' integrator, suitable for diffusion-reaction systems.

Parameters:

Name	Type	Description	Default
method	('Radau', 'BDF', 'RK45', 'RK23', 'DOP853', 'LSODA')	ODE integration algorithm passed to scipy.integrate.solve_ivp. Recommended options: 'Radau' (default): 5th-order implicit Runge-Kutta; highly recommended for stiff diffusion and fast reaction-diffusion equations. 'BDF': Variable-order multi-step Backward Differentiation Formula; well-suited for large stiff systems. 'RK45': Explicit Runge-Kutta order 4(5); suitable for non-stiff pure diffusion.	'Radau'
rtol	float	Relative error tolerance for adaptive time stepping.	1e-6
atol	float	Absolute error tolerance for adaptive time stepping.	1e-8
**options	Any	Additional solver options passed to scipy.integrate.solve_ivp.	{}

Examples:

Simulating reaction-diffusion with first-order decay:

>>> import numpy as np
>>> from softpotato.solver import ScipyIVPSolver, DiffusionProblem, NeumannBC
>>>
>>> grid = np.linspace(0.0, 1.0, 31)
>>> problem = DiffusionProblem(
...     grid=grid,
...     diffusivity={"A": 0.01},
...     boundary_conditions={"A": (NeumannBC(0.0), NeumannBC(0.0))},
...     initial_conditions={"A": 1.0},
...     reactions=lambda c_dict, x, t: {"A": -0.5 * c_dict["A"]},
... )
>>> solver = ScipyIVPSolver(method="Radau", rtol=1e-5)
>>> result = solver.solve(problem, t_span=(0.0, 1.0))
>>> result.success
True
>>> result["A"].shape[1]
31

Source code in src/softpotato/solver/scipy_ivp.py

class ScipyIVPSolver(BaseSolver, name="scipy_ivp"):
    r"""Diffusion solver using Method of Lines (MOL) and ``scipy.integrate.solve_ivp``.

    Discretizes the spatial diffusion operator into a coupled system of ODEs:

    .. math::

        \frac{dc_i}{dt} = D_i \frac{\partial^2 c_i}{\partial x^2} + R_i(\mathbf{c}, x, t)

    Integrates the resulting stiff or non-stiff system using SciPy's ODE suite.
    Defaults to the stiff L-stable ``'Radau'`` integrator, suitable for diffusion-reaction systems.

    Parameters
    ----------
    method : {'Radau', 'BDF', 'RK45', 'RK23', 'DOP853', 'LSODA'}, default 'Radau'
        ODE integration algorithm passed to ``scipy.integrate.solve_ivp``.
        Recommended options:

        - ``'Radau'`` (default): 5th-order implicit Runge-Kutta; highly recommended for stiff
          diffusion and fast reaction-diffusion equations.
        - ``'BDF'``: Variable-order multi-step Backward Differentiation Formula; well-suited
          for large stiff systems.
        - ``'RK45'``: Explicit Runge-Kutta order 4(5); suitable for non-stiff pure diffusion.
    rtol : float, default 1e-6
        Relative error tolerance for adaptive time stepping.
    atol : float, default 1e-8
        Absolute error tolerance for adaptive time stepping.
    **options : Any
        Additional solver options passed to ``scipy.integrate.solve_ivp``.

    Examples
    --------
    Simulating reaction-diffusion with first-order decay:

    >>> import numpy as np
    >>> from softpotato.solver import ScipyIVPSolver, DiffusionProblem, NeumannBC
    >>>
    >>> grid = np.linspace(0.0, 1.0, 31)
    >>> problem = DiffusionProblem(
    ...     grid=grid,
    ...     diffusivity={"A": 0.01},
    ...     boundary_conditions={"A": (NeumannBC(0.0), NeumannBC(0.0))},
    ...     initial_conditions={"A": 1.0},
    ...     reactions=lambda c_dict, x, t: {"A": -0.5 * c_dict["A"]},
    ... )
    >>> solver = ScipyIVPSolver(method="Radau", rtol=1e-5)
    >>> result = solver.solve(problem, t_span=(0.0, 1.0))
    >>> result.success
    True
    >>> result["A"].shape[1]
    31
    """

    def __init__(
        self,
        method: str = "Radau",
        rtol: float = 1e-6,
        atol: float = 1e-8,
        **options: Any,
    ) -> None:
        """Initialize SciPy Method of Lines IVP solver."""
        super().__init__(method=method, rtol=rtol, atol=atol, **options)

    def _run_solver(
        self,
        problem: DiffusionProblem,
        t_span: tuple[float, float],
        t_eval: np.ndarray | None,
        **kwargs: Any,
    ) -> SolverResult:
        geom = problem.geometry
        assert isinstance(geom, Geometry)
        grid = np.asarray(problem.grid, dtype=float)
        N = problem.n_points
        t_start, t_end = t_span

        # Determine unknown indices per species based on boundary conditions
        species_info: dict[str, dict[str, Any]] = {}
        total_unknowns = 0

        for sp in problem.species:
            bc_left, bc_right = problem.boundary_conditions[sp]
            idx_start = 1 if bc_left.btype == "dirichlet" else 0
            idx_end = N - 1 if bc_right.btype == "dirichlet" else N
            size = idx_end - idx_start

            species_info[sp] = {
                "idx_start": idx_start,
                "idx_end": idx_end,
                "size": size,
                "offset": total_unknowns,
                "left_dirichlet": bc_left.btype == "dirichlet",
                "right_dirichlet": bc_right.btype == "dirichlet",
            }
            total_unknowns += size

        # Assemble initial state vector y0
        y0 = np.empty(total_unknowns, dtype=float)
        for sp in problem.species:
            info = species_info[sp]
            ic_full = problem.get_initial_profile(sp)
            y0[info["offset"] : info["offset"] + info["size"]] = ic_full[
                info["idx_start"] : info["idx_end"]
            ]

        def rhs(t: float, y: np.ndarray) -> np.ndarray:
            """Compute ODE system time derivatives."""
            # 1. Unpack y into full concentration profiles
            c_dict: dict[str, np.ndarray] = {}
            for sp in problem.species:
                info = species_info[sp]
                c_full = np.empty(N, dtype=float)
                bc_left, bc_right = problem.boundary_conditions[sp]

                c_full[info["idx_start"] : info["idx_end"]] = y[
                    info["offset"] : info["offset"] + info["size"]
                ]
                if info["left_dirichlet"]:
                    c_full[0] = bc_left.evaluate(t)
                if info["right_dirichlet"]:
                    c_full[-1] = bc_right.evaluate(t)

                c_dict[sp] = c_full

            # 2. Evaluate chemical reaction rates if specified
            reactions = (
                problem.reactions(c_dict, grid, t) if problem.reactions else None
            )

            # 3. Compute spatial derivatives and fill dy
            dy = np.empty(total_unknowns, dtype=float)
            for sp in problem.species:
                info = species_info[sp]
                c_full = c_dict[sp]
                D = problem.diffusivity[sp]
                bc_left, bc_right = problem.boundary_conditions[sp]

                lap = geom.laplacian(
                    c_full, grid, bc_left=bc_left, bc_right=bc_right, t=t
                )
                sp_rate = reactions[sp] if reactions and sp in reactions else 0.0
                sp_dcdt = D * lap + sp_rate

                dy[info["offset"] : info["offset"] + info["size"]] = sp_dcdt[
                    info["idx_start"] : info["idx_end"]
                ]

            return dy

        # Solver options
        method = (
            kwargs.get("method")
            if kwargs.get("method") is not None
            else self.options.get("method", "Radau")
        )
        rtol = (
            kwargs.get("rtol")
            if kwargs.get("rtol") is not None
            else self.options.get("rtol", 1e-6)
        )
        atol = (
            kwargs.get("atol")
            if kwargs.get("atol") is not None
            else self.options.get("atol", 1e-8)
        )

        if t_eval is None:
            # Default to 101 points across t_span
            t_eval = np.linspace(t_start, t_end, 101)
        else:
            t_eval = np.asarray(t_eval, dtype=float)

        # Call scipy.integrate.solve_ivp
        ode_res = solve_ivp(
            fun=rhs,
            t_span=t_span,
            y0=y0,
            method=method,
            t_eval=t_eval,
            rtol=rtol,
            atol=atol,
        )

        if not ode_res.success:
            return SolverResult(
                t=t_eval,
                x=grid,
                concentrations={},
                success=False,
                message=f"solve_ivp failed: {ode_res.message}",
                raw_output=ode_res,
            )

        # Unpack solution into (n_times, N) per species
        n_times = len(ode_res.t)
        records: dict[str, np.ndarray] = {}

        for sp in problem.species:
            info = species_info[sp]
            bc_left, bc_right = problem.boundary_conditions[sp]
            sol_sp = np.empty((n_times, N), dtype=float)

            # Insert solved unknowns
            y_sp = ode_res.y[info["offset"] : info["offset"] + info["size"], :]
            sol_sp[:, info["idx_start"] : info["idx_end"]] = y_sp.T

            # Fill Dirichlet boundaries at evaluated times
            if info["left_dirichlet"]:
                sol_sp[:, 0] = [bc_left.evaluate(ti) for ti in ode_res.t]
            if info["right_dirichlet"]:
                sol_sp[:, -1] = [bc_right.evaluate(ti) for ti in ode_res.t]

            records[sp] = sol_sp

        return SolverResult(
            t=ode_res.t,
            x=grid,
            concentrations=records,
            success=True,
            message="SciPy solve_ivp completed successfully.",
            raw_output=ode_res,
        )

__init__

__init__(method: str = 'Radau', rtol: float = 1e-06, atol: float = 1e-08, **options: Any) -> None

Initialize SciPy Method of Lines IVP solver.

Source code in src/softpotato/solver/scipy_ivp.py

def __init__(
    self,
    method: str = "Radau",
    rtol: float = 1e-6,
    atol: float = 1e-8,
    **options: Any,
) -> None:
    """Initialize SciPy Method of Lines IVP solver."""
    super().__init__(method=method, rtol=rtol, atol=atol, **options)

DirichletBC

DirichletBC(value: float | Callable[[float], float]) -> BoundaryCondition

Convenience constructor for a Dirichlet boundary condition \(c = \text{value}\).

Parameters:

Name	Type	Description	Default
value	float or Callable[[float], float]	Fixed concentration or time-dependent function \(f(t)\).	required

Returns:

Type	Description
BoundaryCondition	A boundary condition instance configured with btype='dirichlet'.

Examples:

>>> from softpotato.solver import DirichletBC
>>> bc = DirichletBC(0.0)
>>> bc.btype
'dirichlet'
>>> bc.evaluate(10.0)
0.0

Source code in src/softpotato/solver/base.py

def DirichletBC(value: float | Callable[[float], float]) -> BoundaryCondition:
    r"""Convenience constructor for a Dirichlet boundary condition :math:`c = \text{value}`.

    Parameters
    ----------
    value : float or Callable[[float], float]
        Fixed concentration or time-dependent function :math:`f(t)`.

    Returns
    -------
    BoundaryCondition
        A boundary condition instance configured with ``btype='dirichlet'``.

    Examples
    --------
    >>> from softpotato.solver import DirichletBC
    >>> bc = DirichletBC(0.0)
    >>> bc.btype
    'dirichlet'
    >>> bc.evaluate(10.0)
    0.0
    """
    return BoundaryCondition(btype="dirichlet", value=value)

NeumannBC

NeumannBC(value: float | Callable[[float], float]) -> BoundaryCondition

Convenience constructor for a Neumann boundary condition \(\frac{\partial c}{\partial x} = \text{value}\).

Parameters:

Name	Type	Description	Default
value	float or Callable[[float], float]	Fixed gradient \(\partial c / \partial x\) or time-dependent function \(g(t)\).	required

Returns:

Type	Description
BoundaryCondition	A boundary condition instance configured with btype='neumann'.

Examples:

>>> from softpotato.solver import NeumannBC
>>> bc = NeumannBC(0.0)
>>> bc.btype
'neumann'
>>> bc.evaluate(0.0)
0.0

Source code in src/softpotato/solver/base.py

def NeumannBC(value: float | Callable[[float], float]) -> BoundaryCondition:
    r"""Convenience constructor for a Neumann boundary condition :math:`\frac{\partial c}{\partial x} = \text{value}`.

    Parameters
    ----------
    value : float or Callable[[float], float]
        Fixed gradient :math:`\partial c / \partial x` or time-dependent function :math:`g(t)`.

    Returns
    -------
    BoundaryCondition
        A boundary condition instance configured with ``btype='neumann'``.

    Examples
    --------
    >>> from softpotato.solver import NeumannBC
    >>> bc = NeumannBC(0.0)
    >>> bc.btype
    'neumann'
    >>> bc.evaluate(0.0)
    0.0
    """
    return BoundaryCondition(btype="neumann", value=value)

get_solver

get_solver(name: str, **kwargs: Any) -> BaseSolver

Instantiate a solver by its registered name.

Parameters:

Name	Type	Description	Default
name	str	Registered solver identifier (case-insensitive), e.g. 'scipy_ivp', 'crank_nicolson', 'implicit', 'explicit'.	required
**kwargs	Any	Configuration options forwarded to the solver's constructor: For 'explicit' / 'ftcs': dt (float, optional): Fixed time step size satisfying the CFL stability limit. For 'implicit' / 'btcs': dt (float, optional): Time step size (defaults to 100 intervals). For 'crank_nicolson': dt (float, optional): Time step size (defaults to 100 intervals). For 'scipy_ivp' / 'solve_ivp': method (str, default 'Radau'): Integration algorithm ('Radau', 'BDF', 'RK45', etc.). rtol (float, default 1e-6): Relative error tolerance. atol (float, default 1e-8): Absolute error tolerance.	{}

Returns:

Type	Description
BaseSolver	An instance of the requested solver.

Examples:

>>> from softpotato.solver import get_solver
>>> solver = get_solver("scipy_ivp", method="Radau")
>>> solver.options["method"]
'Radau'
>>> solver_cn = get_solver("crank_nicolson", dt=0.01)
>>> solver_cn.options["dt"]
0.01

Source code in src/softpotato/solver/base.py

def get_solver(name: str, **kwargs: Any) -> BaseSolver:
    """Instantiate a solver by its registered name.

    Parameters
    ----------
    name : str
        Registered solver identifier (case-insensitive), e.g. 'scipy_ivp', 'crank_nicolson',
        'implicit', 'explicit'.
    **kwargs : Any
        Configuration options forwarded to the solver's constructor:

        - For ``'explicit'`` / ``'ftcs'``:
          - ``dt`` (float, optional): Fixed time step size satisfying the CFL stability limit.
        - For ``'implicit'`` / ``'btcs'``:
          - ``dt`` (float, optional): Time step size (defaults to 100 intervals).
        - For ``'crank_nicolson'``:
          - ``dt`` (float, optional): Time step size (defaults to 100 intervals).
        - For ``'scipy_ivp'`` / ``'solve_ivp'``:
          - ``method`` (str, default 'Radau'): Integration algorithm ('Radau', 'BDF', 'RK45', etc.).
          - ``rtol`` (float, default 1e-6): Relative error tolerance.
          - ``atol`` (float, default 1e-8): Absolute error tolerance.

    Returns
    -------
    BaseSolver
        An instance of the requested solver.

    Examples
    --------
    >>> from softpotato.solver import get_solver
    >>> solver = get_solver("scipy_ivp", method="Radau")
    >>> solver.options["method"]
    'Radau'
    >>> solver_cn = get_solver("crank_nicolson", dt=0.01)
    >>> solver_cn.options["dt"]
    0.01
    """
    normalized = name.lower()
    if normalized not in BaseSolver._registry:
        available = ", ".join(sorted(BaseSolver._registry.keys()))
        raise KeyError(f"Unknown solver '{name}'. Available solvers: {available}")
    return BaseSolver._registry[normalized](**kwargs)

list_solvers

list_solvers() -> list[str]

Return a list of all currently registered solver names.

Returns:

Type	Description
list[str]	Alphabetically sorted list of registered solver names.

Examples:

>>> from softpotato.solver import list_solvers
>>> solvers = list_solvers()
>>> "scipy_ivp" in solvers
True
>>> "crank_nicolson" in solvers
True

Source code in src/softpotato/solver/base.py

def list_solvers() -> list[str]:
    """Return a list of all currently registered solver names.

    Returns
    -------
    list[str]
        Alphabetically sorted list of registered solver names.

    Examples
    --------
    >>> from softpotato.solver import list_solvers
    >>> solvers = list_solvers()
    >>> "scipy_ivp" in solvers
    True
    >>> "crank_nicolson" in solvers
    True
    """
    return sorted(BaseSolver._registry.keys())

register_solver

register_solver(name: str) -> Callable[[type[BaseSolver]], type[BaseSolver]]

Decorator to register a custom solver class under a given name.

Parameters:

Name	Type	Description	Default
name	str	Name identifier to associate with the decorated solver class.	required

Returns:

Type	Description
Callable[[type[BaseSolver]], type[BaseSolver]]	Decorator registering the class into BaseSolver._registry.

Examples:

>>> import numpy as np
>>> from softpotato.solver import BaseSolver, SolverResult, register_solver, get_solver
>>>
>>> @register_solver("spectral_decorator_demo")
... class SpectralSolver(BaseSolver):
...     def _run_solver(self, problem, t_span, t_eval, **kwargs):
...         return SolverResult(t=np.array([0, 1]), x=problem.grid, concentrations={})
>>>
>>> solver = get_solver("spectral_decorator_demo")
>>> isinstance(solver, SpectralSolver)
True

Source code in src/softpotato/solver/base.py

def register_solver(name: str) -> Callable[[type[BaseSolver]], type[BaseSolver]]:
    """Decorator to register a custom solver class under a given name.

    Parameters
    ----------
    name : str
        Name identifier to associate with the decorated solver class.

    Returns
    -------
    Callable[[type[BaseSolver]], type[BaseSolver]]
        Decorator registering the class into BaseSolver._registry.

    Examples
    --------
    >>> import numpy as np
    >>> from softpotato.solver import BaseSolver, SolverResult, register_solver, get_solver
    >>>
    >>> @register_solver("spectral_decorator_demo")
    ... class SpectralSolver(BaseSolver):
    ...     def _run_solver(self, problem, t_span, t_eval, **kwargs):
    ...         return SolverResult(t=np.array([0, 1]), x=problem.grid, concentrations={})
    >>>
    >>> solver = get_solver("spectral_decorator_demo")
    >>> isinstance(solver, SpectralSolver)
    True
    """

    def decorator(cls: type[BaseSolver]) -> type[BaseSolver]:
        if not issubclass(cls, BaseSolver):
            raise TypeError(f"Class {cls.__name__} must inherit from BaseSolver.")
        BaseSolver._registry[name.lower()] = cls
        return cls

    return decorator