Solver Module (`softpotato.solver`)¶

The softpotato.solver subpackage provides numerical solvers for 1D multi-species chemical diffusion and reaction-diffusion problems.

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^2 / 2D_{\max}\)).
`ImplicitFiniteDifference` (`'implicit'`, `'btcs'`)	`dt` (float, optional)	Backward-Time Central-Space implicit scheme. 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. 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.

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 and perform post-run integrity checks."""
        # Calculate surface flux J = -D * dc/dx at left boundary x=0 using 2nd order difference.
        # Sign convention:
        # Fick's first law defines the diffusive flux along the +x outward normal into the solution:
        #   J(0, t) = -D * (dc/dx)|_{x=0}
        # When an electroactive species is consumed at the electrode (c(0, t) < c_bulk), (dc/dx)|_{x=0} > 0,
        # so the interfacial flux J(0, t) is strictly NEGATIVE (J < 0).
        # Under standard IUPAC sign conventions (cathodic reduction < 0, anodic oxidation > 0):
        # - For reduction (Ox + n e^- -> Red):
        #     I_cathodic = n * F * A * J_Ox < 0 (negative)
        # - For oxidation (Red -> Ox + n e^-):
        #     I_anodic = -n * F * A * J_Red > 0 (positive)
        if problem.is_uniform:
            dx = problem.dx
            for sp in problem.species:
                if sp not in result.fluxes:
                    c = result.concentrations[sp]
                    D = problem.diffusivity[sp]
                    # 2nd-order forward difference at x=0: (-3*c0 + 4*c1 - c2) / (2*dx)
                    dcdx_left = (-3.0 * c[:, 0] + 4.0 * c[:, 1] - c[:, 2]) / (2.0 * dx)
                    result.fluxes[sp] = -D * dcdx_left
        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`	1D spatial grid points \(x \in [x_{\min}, x_{\max}]\).	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()`

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
        1D spatial grid points :math:`x \in [x_{\min}, x_{\max}]`.
    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.

    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
    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)

    def __post_init__(self) -> None:
        self.grid = np.asarray(self.grid, dtype=float)
        if self.grid.ndim != 1 or len(self.grid) < 3:
            raise ValueError("grid must be a 1D array with at least 3 points.")
        if np.any(np.diff(self.grid) <= 0):
            raise ValueError("grid points must be strictly monotonically increasing.")

        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 dx(self) -> float:
        """Grid step size for uniform grids."""
        diffs = np.diff(self.grid)
        if not np.allclose(diffs, diffs[0], rtol=1e-5, atol=1e-8):
            raise ValueError("Grid is non-uniform; uniform dx is not defined.")
        return float(diffs[0])

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

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

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

n_points `property` ¶

n_points: int

Number of spatial grid points.

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."""
    ic = self.initial_conditions[sp]
    if callable(ic):
        arr = np.asarray(ic(self.grid), dtype=float)
    elif np.isscalar(ic):
        arr = np.full_like(self.grid, float(cast(Any, ic)), dtype=float)
    else:
        arr = np.asarray(ic, dtype=float)

    if arr.shape != self.grid.shape:
        raise ValueError(
            f"Initial condition for '{sp}' has shape {arr.shape}, expected {self.grid.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()`
`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)},
... )
>>> res.species
['A']
>>> res["A"].shape
(5, 10)
>>> res.fluxes["A"].shape
(5,)

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,).
    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)},
    ... )
    >>> res.species
    ['A']
    >>> res["A"].shape
    (5, 10)
    >>> res.fluxes["A"].shape
    (5,)
    """

    t: np.ndarray
    x: np.ndarray
    concentrations: dict[str, np.ndarray]
    fluxes: dict[str, np.ndarray] = field(default_factory=dict)
    success: bool = True
    message: str = "Success"
    raw_output: Any = 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 species(self) -> list[str]:
        """List of species present in the solution."""
        return list(self.concentrations.keys())

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:
        if not problem.is_uniform:
            raise NotImplementedError("CrankNicolson requires a uniform spatial grid.")

        dx = problem.dx
        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}

        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_crank_nicolson(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=problem.grid,
            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,
    ) -> dict[str, np.ndarray]:
        """Advance all species profiles by one Crank-Nicolson time step."""
        reactions = (
            problem.reactions(c_dict, problem.grid, t) if problem.reactions else None
        )
        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]
            r = D * dt / (2.0 * dx**2)

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

            # Interior rows
            ab[1, 1:-1] = 1.0 + 2.0 * r
            ab[0, 2:] = -r  # A[i, i+1]
            ab[2, :-2] = -r  # A[i, i-1]
            rhs[1:-1] = r * c[:-2] + (1.0 - 2.0 * r) * c[1:-1] + r * c[2:]
            if reactions and sp in reactions:
                rhs[1:-1] += dt * reactions[sp][1:-1]

            # Left boundary
            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 + 2.0 * r
                ab[0, 1] = -2.0 * r
                g_left_curr = bc_left.evaluate(t)
                g_left_next = bc_left.evaluate(t_next)
                rhs[0] = (
                    (1.0 - 2.0 * r) * c[0]
                    + 2.0 * r * c[1]
                    - 2.0 * r * dx * (g_left_curr + g_left_next)
                )
                if reactions and sp in reactions:
                    rhs[0] += dt * reactions[sp][0]

            # Right boundary
            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 + 2.0 * r
                ab[2, -2] = -2.0 * r
                g_right_curr = bc_right.evaluate(t)
                g_right_next = bc_right.evaluate(t_next)
                rhs[-1] = (
                    (1.0 - 2.0 * r) * c[-1]
                    + 2.0 * r * c[-2]
                    + 2.0 * r * dx * (g_right_curr + g_right_next)
                )
                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:

\[ \Delta t \le \frac{\Delta x^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^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:

    .. math::

        \Delta t \le \frac{\Delta x^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^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:
        if not problem.is_uniform:
            raise NotImplementedError(
                "ExplicitFiniteDifference requires a uniform spatial grid."
            )

        dx = problem.dx
        max_D = max(problem.diffusivity.values())
        cfl_limit = (dx**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²/(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}

        # 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=problem.grid,
            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,
    ) -> dict[str, np.ndarray]:
        """Advance all species concentration profiles by one time step dt."""
        reactions = (
            problem.reactions(c_dict, problem.grid, t) if problem.reactions else None
        )
        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]

            d2c = np.empty_like(c)
            # Interior 2nd-derivative
            d2c[1:-1] = (c[2:] - 2.0 * c[1:-1] + c[:-2]) / (dx**2)

            # Left boundary
            if bc_left.btype == "neumann":
                g_left = bc_left.evaluate(t)
                d2c[0] = (2.0 * c[1] - 2.0 * c[0] - 2.0 * dx * g_left) / (dx**2)
            else:
                d2c[0] = 0.0

            # Right boundary
            if bc_right.btype == "neumann":
                g_right = bc_right.evaluate(t)
                d2c[-1] = (2.0 * c[-2] - 2.0 * c[-1] + 2.0 * dx * g_right) / (dx**2)
            else:
                d2c[-1] = 0.0

            # Advance state
            sp_reaction = reactions[sp] if reactions and sp in reactions else 0.0
            next_sp_c = c + dt * (D * d2c + 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:
        if not problem.is_uniform:
            raise NotImplementedError(
                "ImplicitFiniteDifference requires a uniform spatial grid."
            )

        dx = problem.dx
        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}

        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)
                t_current += dt

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

        return SolverResult(
            t=t_eval,
            x=problem.grid,
            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,
    ) -> dict[str, np.ndarray]:
        """Advance all species profiles by one BTCS implicit time step."""
        reactions = (
            problem.reactions(c_dict, problem.grid, t) if problem.reactions else None
        )
        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]
            r = D * dt / (dx**2)

            # Assemble banded matrix ab of shape (3, N)
            # ab[0, 1:] = upper diagonal
            # ab[1, :]  = main diagonal
            # ab[2, :-1] = lower diagonal
            ab = np.zeros((3, N), dtype=float)
            rhs = np.copy(c)
            if reactions and sp in reactions:
                rhs += dt * reactions[sp]

            # Interior rows
            ab[1, 1:-1] = 1.0 + 2.0 * r
            ab[0, 2:] = -r  # A[i, i+1]
            ab[2, :-2] = -r  # A[i, i-1]

            # 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 + 2.0 * r
                ab[0, 1] = -2.0 * r
                g_left = bc_left.evaluate(t_next)
                rhs[0] -= 2.0 * r * dx * 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 + 2.0 * r
                ab[2, -2] = -2.0 * r
                g_right = bc_right.evaluate(t_next)
                rhs[-1] += 2.0 * r * dx * 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:
        grid = problem.grid
        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"]
            ]

        # Precompute spatial discretization weights
        is_uniform = problem.is_uniform
        if is_uniform:
            dx = problem.dx
            inv_dx2 = 1.0 / (dx**2)
        else:
            h = np.diff(grid)  # h[i] = x[i+1] - x[i]

        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]

                d2c = np.empty(N, dtype=float)
                if is_uniform:
                    # Interior points
                    d2c[1:-1] = (
                        c_full[2:] - 2.0 * c_full[1:-1] + c_full[:-2]
                    ) * inv_dx2

                    # Boundary points if Neumann
                    if not info["left_dirichlet"]:
                        g_left = bc_left.evaluate(t)
                        d2c[0] = (
                            2.0 * c_full[1] - 2.0 * c_full[0] - 2.0 * dx * g_left
                        ) * inv_dx2

                    if not info["right_dirichlet"]:
                        g_right = bc_right.evaluate(t)
                        d2c[-1] = (
                            2.0 * c_full[-2] - 2.0 * c_full[-1] + 2.0 * dx * g_right
                        ) * inv_dx2
                else:
                    # General 3-point stencil for non-uniform grid
                    # d²c/dx² ~ 2 / (h_{i-1} + h_i) * [ (c_{i+1}-c_i)/h_i - (c_i-c_{i-1})/h_{i-1} ]
                    hi_m1 = h[:-1]
                    hi = h[1:]
                    d2c[1:-1] = (2.0 / (hi_m1 + hi)) * (
                        (c_full[2:] - c_full[1:-1]) / hi
                        - (c_full[1:-1] - c_full[:-2]) / hi_m1
                    )
                    if not info["left_dirichlet"]:
                        g_left = bc_left.evaluate(t)
                        h0 = h[0]
                        d2c[0] = (2.0 / h0) * ((c_full[1] - c_full[0]) / h0 - g_left)
                    if not info["right_dirichlet"]:
                        g_right = bc_right.evaluate(t)
                        hm1 = h[-1]
                        d2c[-1] = (2.0 / hm1) * (
                            g_right - (c_full[-1] - c_full[-2]) / hm1
                        )

                sp_rate = reactions[sp] if reactions and sp in reactions else 0.0
                sp_dcdt = D * d2c + 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

Solver Module (softpotato.solver)¶

Solvers and Supported Options¶

Module Contents¶

solver ¶

BaseSolver ¶

__init__ ¶

solve ¶

validate ¶

BoundaryCondition dataclass ¶

evaluate ¶

DiffusionProblem dataclass ¶

n_points property ¶

dx property ¶

is_uniform property ¶

get_initial_profile ¶

SolverResult dataclass ¶

species property ¶

__getitem__ ¶

CrankNicolson ¶

__init__ ¶

ExplicitFiniteDifference ¶

__init__ ¶

ImplicitFiniteDifference ¶

__init__ ¶

ScipyIVPSolver ¶

__init__ ¶

DirichletBC ¶

NeumannBC ¶

get_solver ¶

list_solvers ¶

register_solver ¶

Solver Module (`softpotato.solver`)¶

init ¶

BoundaryCondition `dataclass` ¶

DiffusionProblem `dataclass` ¶

n_points `property` ¶

dx `property` ¶

is_uniform `property` ¶

SolverResult `dataclass` ¶

species `property` ¶

getitem ¶

init ¶

init ¶

init ¶

init ¶