Geometry Module (softpotato.geometry)

The softpotato.geometry subpackage provides spatial coordinate abstractions, finite-difference differential operators (\(\nabla^2\), \(\nabla\)), normal interfacial flux evaluation, and active electrode surface integration across planar, spherical, and cylindrical domains.

Coordinate Geometries and Differential Operators

Geometry Class & Aliases	Symmetry Factor (\(p\))	Laplacian Operator (\(\nabla^2 c\))	Active Surface Area (\(A\))
PlanarGeometry ('planar', 'cartesian', '1d', 'linear')	\(p = 0\)	\(\dfrac{\partial^2 c}{\partial x^2}\)	User-specified \(A\) (default \(1.0\))
CylindricalGeometry ('cylindrical', 'cylinder', 'wire')	\(p = 1\)	\(\dfrac{\partial^2 c}{\partial r^2} + \dfrac{1}{r}\dfrac{\partial c}{\partial r}\)	\(2 \pi r_0 L\)
SphericalGeometry ('spherical', 'sphere', 'hemispherical', 'hemisphere')	\(p = 2\)	\(\dfrac{\partial^2 c}{\partial r^2} + \dfrac{2}{r}\dfrac{\partial c}{\partial r}\)	\(4 \pi r_0^2\) (full sphere) or \(2 \pi r_0^2\) (hemisphere)

Module Contents

geometry

Spatial coordinate geometries and differential operators for electrochemical simulation.

This subpackage abstracts coordinate system metrics, spatial finite-difference operators (\(\nabla^2\), \(\nabla\)), normal interfacial flux evaluation, and active electrode surface integration:

\[ \nabla^2 c = \frac{\partial^2 c}{\partial r^2} + \frac{p}{r} \frac{\partial c}{\partial r} \]

Available built-in geometries include:

• PlanarGeometry ("planar", \(p = 0\)): 1D Cartesian linear diffusion perpendicular to a planar macroelectrode.
• SphericalGeometry ("spherical", "hemispherical", \(p = 2\)): 1D radial diffusion to a full spherical or hemispherical electrode.
• CylindricalGeometry ("cylindrical", \(p = 1\)): 1D axisymmetric radial diffusion to a cylindrical wire or fiber microelectrode.

Examples:

Evaluating spatial operators across planar, spherical, and cylindrical geometries:

>>> import numpy as np
>>> from softpotato.geometry import (
...     PlanarGeometry,
...     SphericalGeometry,
...     CylindricalGeometry,
...     get_geometry,
... )
>>>
>>> planar = get_geometry("planar", area=1.0)
>>> x = np.linspace(0.0, 1.0, 51)
>>> bool(np.allclose(planar.laplacian(x**2, x), 2.0))
True
>>>
>>> sphere = SphericalGeometry(r0=1e-4)
>>> r = np.linspace(1.0, 5.0, 201)
>>> bool(np.max(np.abs(sphere.laplacian(1.0 / r, r)[1:-1])) < 1e-4)
True
>>>
>>> cylinder = CylindricalGeometry(r0=5e-6, length=1e-2)
>>> bool(np.max(np.abs(cylinder.laplacian(np.log(r), r)[1:-1])) < 1e-4)
True

Geometry

Bases: ABC

Abstract base class representing spatial domain coordinate geometry.

Encapsulates coordinate metrics, finite-difference differential operators (\(\nabla^2\), \(\nabla\)), normal interfacial flux evaluation, and active electrode surface area integration for 1D and curvilinear coordinate systems:

\[ \nabla^2 c = \frac{1}{r^p} \frac{\partial}{\partial r} \left( r^p \frac{\partial c}{\partial r} \right) = \frac{\partial^2 c}{\partial r^2} + \frac{p}{r} \frac{\partial c}{\partial r} \]

where \(p = 0\) corresponds to planar Cartesian coordinates (\(x\)), \(p = 1\) to cylindrical radial coordinates (\(r\)), and \(p = 2\) to spherical radial coordinates (\(r\)).

Attributes:

Name	Type	Description
name	str	Canonical identifier of the coordinate geometry (e.g., 'planar', 'spherical', 'cylindrical').
dimension	int	Number of active spatial dimensions discretized by the geometry.
coordinate_name	str	Symbol of the primary spatial coordinate (e.g., 'x' or 'r').
p	int	Radial symmetry factor (\(p = 0\) planar, \(p = 1\) cylindrical, \(p = 2\) spherical).

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

class Geometry(ABC, metaclass=_GeometryMeta):
    r"""Abstract base class representing spatial domain coordinate geometry.

    Encapsulates coordinate metrics, finite-difference differential operators
    (:math:`\nabla^2`, :math:`\nabla`), normal interfacial flux evaluation, and
    active electrode surface area integration for 1D and curvilinear coordinate
    systems:

    .. math::

        \nabla^2 c = \frac{1}{r^p} \frac{\partial}{\partial r}
                     \left( r^p \frac{\partial c}{\partial r} \right)
                   = \frac{\partial^2 c}{\partial r^2} + \frac{p}{r} \frac{\partial c}{\partial r}

    where :math:`p = 0` corresponds to planar Cartesian coordinates (:math:`x`),
    :math:`p = 1` to cylindrical radial coordinates (:math:`r`), and :math:`p = 2`
    to spherical radial coordinates (:math:`r`).

    Attributes
    ----------
    name : str
        Canonical identifier of the coordinate geometry (e.g., ``'planar'``,
        ``'spherical'``, ``'cylindrical'``).
    dimension : int
        Number of active spatial dimensions discretized by the geometry.
    coordinate_name : str
        Symbol of the primary spatial coordinate (e.g., ``'x'`` or ``'r'``).
    p : int
        Radial symmetry factor (:math:`p = 0` planar, :math:`p = 1` cylindrical,
        :math:`p = 2` spherical).
    """

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

    name: str = "geometry"
    dimension: int = 1
    coordinate_name: str = "x"
    p: int = 0

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

    @property
    @abstractmethod
    def area(self) -> float:
        r"""Active electrode surface area :math:`A` in :math:`\text{m}^2` or :math:`\text{cm}^2`."""
        raise NotImplementedError

    @property
    def surface_area(self) -> float:
        r"""Active electrode surface area :math:`A` (alias for :attr:`area`)."""
        return self.area

    @property
    def symmetry_factor(self) -> int:
        r"""Radial coordinate symmetry exponent :math:`p \in \{0, 1, 2\}`."""
        return self.p

    def validate_grid(self, grid: np.ndarray | Mesh1D) -> np.ndarray:
        """Validate that a spatial grid satisfies the coordinate domain constraints.

        Parameters
        ----------
        grid : np.ndarray or Mesh1D
            1D array or mesh of spatial coordinates.

        Returns
        -------
        np.ndarray
            Validated 1D floating-point coordinate grid.

        Raises
        ------
        ValueError
            If ``grid`` is not 1D, has fewer than 3 points, is not strictly
            monotonically increasing, or violates radial domain bounds (:math:`r_0 > 0`).
        """
        arr = np.asarray(grid.nodes if isinstance(grid, Mesh1D) else grid, dtype=float)
        if arr.ndim != 1 or len(arr) < 3:
            raise ValueError("grid must be a 1D array with at least 3 points.")
        if np.any(np.diff(arr) <= 0.0):
            raise ValueError("grid points must be strictly monotonically increasing.")
        if self.p > 0 and arr[0] <= 0.0:
            raise ValueError(
                f"{self.__class__.__name__} requires strictly positive radial coordinates "
                f"(r0 = grid[0] > 0), got grid[0]={arr[0]}."
            )
        return arr

    @abstractmethod
    def laplacian(
        self,
        c: np.ndarray,
        grid: np.ndarray | Mesh1D,
        bc_left: Any = None,
        bc_right: Any = None,
        t: float = 0.0,
    ) -> np.ndarray:
        r"""Evaluate the spatial Laplacian :math:`\nabla^2 c` at interior and boundary points.

        Parameters
        ----------
        c : np.ndarray
            1D concentration profile array of shape ``(N,)``.
        grid : np.ndarray or Mesh1D
            1D spatial coordinate array or mesh of shape ``(N,)``.
        bc_left : BoundaryCondition, optional
            Boundary condition at the inner/left boundary (:math:`i = 0`).
        bc_right : BoundaryCondition, optional
            Boundary condition at the outer/right boundary (:math:`i = N - 1`).
        t : float, default 0.0
            Evaluation time in seconds for time-dependent boundary conditions.

        Returns
        -------
        np.ndarray
            1D array of :math:`\nabla^2 c` values of shape ``(N,)``.
        """
        raise NotImplementedError

    @abstractmethod
    def laplacian_matrix(
        self,
        grid: np.ndarray | Mesh1D,
        bc_left: Any = None,
        bc_right: Any = None,
    ) -> tuple[np.ndarray, np.ndarray]:
        r"""Return banded matrix representation of :math:`\nabla^2` for implicit solvers.

        Constructs the tridiagonal finite-difference operator :math:`\mathbf{L}` in
        ``scipy.linalg.solve_banded`` ``(1, 1)`` format of shape ``(3, N)``, together
        with a boundary gradient scaling vector :math:`\mathbf{b}_{\text{scale}}` of
        shape ``(N,)`` such that for Neumann boundary gradients :math:`g_{\text{left}}`,
        :math:`g_{\text{right}}`:

        .. math::

            \nabla^2 \mathbf{c} = \mathbf{L}\mathbf{c}
                + \mathbf{b}_{\text{scale}} \odot [g_{\text{left}}, 0, \dots, 0, g_{\text{right}}]^T

        Parameters
        ----------
        grid : np.ndarray or Mesh1D
            1D spatial coordinate array or mesh of shape ``(N,)``.
        bc_left : BoundaryCondition, optional
            Boundary condition at the left boundary (:math:`i = 0`).
        bc_right : BoundaryCondition, optional
            Boundary condition at the right boundary (:math:`i = N - 1`).

        Returns
        -------
        L_banded : np.ndarray
            Banded tridiagonal matrix of shape ``(3, N)`` where row 0 holds the
            superdiagonal, row 1 holds the main diagonal, and row 2 holds the subdiagonal.
        bc_scale : np.ndarray
            1D boundary gradient multiplier array of shape ``(N,)``.
        """
        raise NotImplementedError

    @abstractmethod
    def surface_flux(
        self,
        c: np.ndarray,
        grid: np.ndarray | Mesh1D,
        diffusivity: float,
        side: str = "left",
    ) -> float | np.ndarray:
        r"""Calculate normal interfacial diffusion flux :math:`J = -D (\partial c / \partial n)` at a boundary.

        Follows the electrode-normal convention where the unit normal :math:`\mathbf{n}`
        points from the boundary into the solution domain (:math:`+\hat{x}` or :math:`+\hat{r}`
        at the left boundary, and :math:`-\hat{x}` or :math:`-\hat{r}` at the right boundary).

        Parameters
        ----------
        c : np.ndarray
            Concentration profile of shape ``(N_x,)`` or time series of shape ``(N_t, N_x)``.
        grid : np.ndarray or Mesh1D
            1D spatial coordinate array or mesh of shape ``(N_x,)``.
        diffusivity : float
            Species diffusion coefficient :math:`D > 0`.
        side : {'left', 'right'}, default 'left'
            Domain boundary at which to evaluate the normal flux.

        Returns
        -------
        float or np.ndarray
            Interfacial flux :math:`J` (scalar ``float`` if ``c`` is 1D, or 1D
            ``np.ndarray`` of shape ``(N_t,)`` if ``c`` is 2D).
        """
        raise NotImplementedError

    @abstractmethod
    def current(
        self,
        flux: float | np.ndarray,
        n_electrons: int = 1,
        side: Literal["left", "right"] = "left",
        grid: np.ndarray | Mesh1D | None = None,
    ) -> float | np.ndarray:
        r"""Integrate normal flux over the active electrode surface to compute Faradaic current.

        Computes the total Faradaic current :math:`I = n F A J` according to IUPAC
        sign conventions (cathodic reduction :math:`I < 0` when reactant flux into
        solution is negative, :math:`J < 0`).

        Parameters
        ----------
        flux : float or np.ndarray
            Normal interfacial diffusion flux :math:`J` at the electrode surface.
        n_electrons : int, default 1
            Number of electrons transferred per reaction (:math:`n > 0`).
        side : {'left', 'right'}, default 'left'
            Domain boundary at which to compute the Faradaic current.
        grid : np.ndarray or Mesh1D, optional
            Spatial grid coordinates used to determine surface area if coordinate-dependent.

        Returns
        -------
        float or np.ndarray
            Faradaic current :math:`I` in Amperes (:math:`\text{A}`).
        """
        raise NotImplementedError

    def _eval_1d_laplacian(
        self,
        c: np.ndarray,
        grid: np.ndarray | Mesh1D,
        bc_left: Any = None,
        bc_right: Any = None,
        t: float = 0.0,
    ) -> np.ndarray:
        """Evaluate 1D Cartesian or curvilinear radial Laplacian with 2nd-order stencils."""
        grid_arr = self.validate_grid(grid)
        c_arr = np.asarray(c, dtype=float)
        if c_arr.shape != grid_arr.shape:
            raise ValueError(
                f"Concentration array shape {c_arr.shape} must match grid shape {grid_arr.shape}."
            )

        n_pts = len(grid_arr)
        diffs = np.diff(grid_arr)
        is_uniform = (
            grid.is_uniform
            if isinstance(grid, Mesh1D)
            else bool(np.allclose(diffs, diffs[0], rtol=1e-5, atol=0.0))
        )
        p = self.p

        lap = np.empty(n_pts, dtype=float)

        # 1. Interior points (i = 1, ..., N - 2)
        if is_uniform:
            dr = float(diffs[0])
            inv_dr2 = 1.0 / (dr**2)
            d2c = (c_arr[2:] - 2.0 * c_arr[1:-1] + c_arr[:-2]) * inv_dr2
            if p == 0:
                lap[1:-1] = d2c
            else:
                dc = (c_arr[2:] - c_arr[:-2]) / (2.0 * dr)
                lap[1:-1] = d2c + (p / grid_arr[1:-1]) * dc
        else:
            hm1 = diffs[:-1]
            hi = diffs[1:]
            d2c = (2.0 / (hm1 + hi)) * (
                (c_arr[2:] - c_arr[1:-1]) / hi - (c_arr[1:-1] - c_arr[:-2]) / hm1
            )
            if p == 0:
                lap[1:-1] = d2c
            else:
                dc = (
                    (hm1**2) * (c_arr[2:] - c_arr[1:-1])
                    + (hi**2) * (c_arr[1:-1] - c_arr[:-2])
                ) / (hm1 * hi * (hm1 + hi))
                lap[1:-1] = d2c + (p / grid_arr[1:-1]) * dc

        # 2. Left boundary (i = 0)
        if bc_left is not None:
            btype_left = getattr(bc_left, "btype", None)
            if btype_left == "dirichlet":
                lap[0] = 0.0
            elif btype_left == "neumann":
                g_left = float(bc_left.evaluate(t))
                h0 = float(diffs[0])
                if is_uniform:
                    d2c_0 = (2.0 * c_arr[1] - 2.0 * c_arr[0] - 2.0 * h0 * g_left) / (
                        h0**2
                    )
                else:
                    d2c_0 = (2.0 / h0) * ((c_arr[1] - c_arr[0]) / h0 - g_left)
                lap[0] = d2c_0 if p == 0 else d2c_0 + (p / grid_arr[0]) * g_left
            else:
                raise ValueError(
                    f"Unsupported left boundary condition type: {btype_left!r}."
                )
        else:
            # Pure operator evaluation at i = 0 using one-sided 2nd-order stencils
            h0, h1 = float(diffs[0]), float(diffs[1])
            dc_0 = (
                -(2.0 * h0 + h1) / (h0 * (h0 + h1)) * c_arr[0]
                + (h0 + h1) / (h0 * h1) * c_arr[1]
                - h0 / (h1 * (h0 + h1)) * c_arr[2]
            )
            if n_pts >= 4:
                d1 = grid_arr[1] - grid_arr[0]
                d2 = grid_arr[2] - grid_arr[0]
                d3 = grid_arr[3] - grid_arr[0]
                d2c_0 = -2.0 * (
                    -(d1 + d2 + d3) / (d1 * d2 * d3) * c_arr[0]
                    + (d2 + d3) / (d1 * (d1 - d2) * (d1 - d3)) * c_arr[1]
                    + (d1 + d3) / (d2 * (d2 - d1) * (d2 - d3)) * c_arr[2]
                    + (d1 + d2) / (d3 * (d3 - d1) * (d3 - d2)) * c_arr[3]
                )
            else:
                d2c_0 = (2.0 / (h0 + h1)) * (
                    (c_arr[2] - c_arr[1]) / h1 - (c_arr[1] - c_arr[0]) / h0
                )
            lap[0] = d2c_0 if p == 0 else d2c_0 + (p / grid_arr[0]) * dc_0

        # 3. Right boundary (i = N - 1)
        if bc_right is not None:
            btype_right = getattr(bc_right, "btype", None)
            if btype_right == "dirichlet":
                lap[-1] = 0.0
            elif btype_right == "neumann":
                g_right = float(bc_right.evaluate(t))
                h_end = float(diffs[-1])
                if is_uniform:
                    d2c_end = (
                        2.0 * c_arr[-2] - 2.0 * c_arr[-1] + 2.0 * h_end * g_right
                    ) / (h_end**2)
                else:
                    d2c_end = (2.0 / h_end) * (
                        g_right - (c_arr[-1] - c_arr[-2]) / h_end
                    )
                lap[-1] = d2c_end if p == 0 else d2c_end + (p / grid_arr[-1]) * g_right
            else:
                raise ValueError(
                    f"Unsupported right boundary condition type: {btype_right!r}."
                )
        else:
            # Pure operator evaluation at i = N - 1 using one-sided 2nd-order stencils
            h_end, h_end2 = float(diffs[-1]), float(diffs[-2])
            dc_end = (
                (2.0 * h_end + h_end2) / (h_end * (h_end + h_end2)) * c_arr[-1]
                - (h_end + h_end2) / (h_end * h_end2) * c_arr[-2]
                + h_end / (h_end2 * (h_end + h_end2)) * c_arr[-3]
            )
            if n_pts >= 4:
                d1 = grid_arr[-2] - grid_arr[-1]
                d2 = grid_arr[-3] - grid_arr[-1]
                d3 = grid_arr[-4] - grid_arr[-1]
                d2c_end = -2.0 * (
                    -(d1 + d2 + d3) / (d1 * d2 * d3) * c_arr[-1]
                    + (d2 + d3) / (d1 * (d1 - d2) * (d1 - d3)) * c_arr[-2]
                    + (d1 + d3) / (d2 * (d2 - d1) * (d2 - d3)) * c_arr[-3]
                    + (d1 + d2) / (d3 * (d3 - d1) * (d3 - d2)) * c_arr[-4]
                )
            else:
                d2c_end = (2.0 / (h_end2 + h_end)) * (
                    (c_arr[-1] - c_arr[-2]) / h_end - (c_arr[-2] - c_arr[-3]) / h_end2
                )
            lap[-1] = d2c_end if p == 0 else d2c_end + (p / grid_arr[-1]) * dc_end

        return lap

    def _build_1d_laplacian_matrix(
        self,
        grid: np.ndarray | Mesh1D,
        bc_left: Any = None,
        bc_right: Any = None,
    ) -> tuple[np.ndarray, np.ndarray]:
        """Assemble banded tridiagonal Laplacian matrix and boundary source scaling."""
        grid_arr = self.validate_grid(grid)
        n_pts = len(grid_arr)
        diffs = np.diff(grid_arr)
        is_uniform = (
            grid.is_uniform
            if isinstance(grid, Mesh1D)
            else bool(np.allclose(diffs, diffs[0], rtol=1e-5, atol=0.0))
        )
        p = self.p

        l_banded = np.zeros((3, n_pts), dtype=float)
        bc_scale = np.zeros(n_pts, dtype=float)

        # Interior rows (i = 1, ..., N - 2)
        if is_uniform:
            dr = float(diffs[0])
            inv_dr2 = 1.0 / (dr**2)
            l_banded[1, 1:-1] = -2.0 * inv_dr2
            if p == 0:
                l_banded[0, 2:] = inv_dr2
                l_banded[2, :-2] = inv_dr2
            else:
                radial_term = p / (2.0 * grid_arr[1:-1] * dr)
                l_banded[0, 2:] = inv_dr2 + radial_term
                l_banded[2, :-2] = inv_dr2 - radial_term
        else:
            hm1 = diffs[:-1]
            hi = diffs[1:]
            denom_sum = hm1 + hi
            w2_lower = 2.0 / (hm1 * denom_sum)
            w2_diag = -2.0 / (hm1 * hi)
            w2_upper = 2.0 / (hi * denom_sum)

            if p == 0:
                l_banded[2, :-2] = w2_lower
                l_banded[1, 1:-1] = w2_diag
                l_banded[0, 2:] = w2_upper
            else:
                p_over_r = p / grid_arr[1:-1]
                w1_lower = -hi / (hm1 * denom_sum)
                w1_diag = (hi - hm1) / (hm1 * hi)
                w1_upper = hm1 / (hi * denom_sum)

                l_banded[2, :-2] = w2_lower + p_over_r * w1_lower
                l_banded[1, 1:-1] = w2_diag + p_over_r * w1_diag
                l_banded[0, 2:] = w2_upper + p_over_r * w1_upper

        # Left boundary (i = 0)
        btype_left = (
            getattr(bc_left, "btype", "neumann") if bc_left is not None else "neumann"
        )
        if btype_left == "dirichlet":
            l_banded[1, 0] = 0.0
            l_banded[0, 1] = 0.0
            bc_scale[0] = 0.0
        elif btype_left == "neumann":
            h0 = float(diffs[0])
            l_banded[1, 0] = -2.0 / (h0**2)
            l_banded[0, 1] = 2.0 / (h0**2)
            bc_scale[0] = -(2.0 / h0 - (p / grid_arr[0] if p > 0 else 0.0))
        else:
            raise ValueError(
                f"Unsupported left boundary condition type: {btype_left!r}."
            )

        # Right boundary (i = N - 1)
        btype_right = (
            getattr(bc_right, "btype", "neumann") if bc_right is not None else "neumann"
        )
        if btype_right == "dirichlet":
            l_banded[1, -1] = 0.0
            l_banded[2, -2] = 0.0
            bc_scale[-1] = 0.0
        elif btype_right == "neumann":
            hm1_end = float(diffs[-1])
            l_banded[1, -1] = -2.0 / (hm1_end**2)
            l_banded[2, -2] = 2.0 / (hm1_end**2)
            bc_scale[-1] = 2.0 / hm1_end + (p / grid_arr[-1] if p > 0 else 0.0)
        else:
            raise ValueError(
                f"Unsupported right boundary condition type: {btype_right!r}."
            )

        return l_banded, bc_scale

    def _eval_1d_surface_flux(
        self,
        c: np.ndarray,
        grid: np.ndarray | Mesh1D,
        diffusivity: float,
        side: str = "left",
    ) -> float | np.ndarray:
        """Compute 2nd-order one-sided interfacial diffusion flux at left or right boundary."""
        if diffusivity <= 0.0:
            raise ValueError(
                f"Diffusivity must be positive (diffusivity > 0), got {diffusivity}."
            )
        if side not in ("left", "right"):
            raise ValueError(f"side must be 'left' or 'right', got {side!r}.")

        grid_arr = self.validate_grid(grid)
        c_arr = np.asarray(c, dtype=float)
        if c_arr.ndim not in (1, 2) or c_arr.shape[-1] != len(grid_arr):
            raise ValueError(
                f"Concentration array shape {c_arr.shape} is incompatible with grid of length {len(grid_arr)}."
            )

        diffs = np.diff(grid_arr)
        is_uniform = (
            grid.is_uniform
            if isinstance(grid, Mesh1D)
            else bool(np.allclose(diffs, diffs[0], rtol=1e-5, atol=0.0))
        )

        if side == "left":
            if is_uniform:
                dx = float(diffs[0])
                grad = (-3.0 * c_arr[..., 0] + 4.0 * c_arr[..., 1] - c_arr[..., 2]) / (
                    2.0 * dx
                )
            else:
                h0, h1 = float(diffs[0]), float(diffs[1])
                grad = (
                    -(2.0 * h0 + h1) / (h0 * (h0 + h1)) * c_arr[..., 0]
                    + (h0 + h1) / (h0 * h1) * c_arr[..., 1]
                    - h0 / (h1 * (h0 + h1)) * c_arr[..., 2]
                )
            flux = -diffusivity * grad
        else:
            if is_uniform:
                dx = float(diffs[0])
                grad = (
                    3.0 * c_arr[..., -1] - 4.0 * c_arr[..., -2] + c_arr[..., -3]
                ) / (2.0 * dx)
            else:
                hm1, hm2 = float(diffs[-1]), float(diffs[-2])
                grad = (
                    (2.0 * hm1 + hm2) / (hm1 * (hm1 + hm2)) * c_arr[..., -1]
                    - (hm1 + hm2) / (hm1 * hm2) * c_arr[..., -2]
                    + hm1 / (hm2 * (hm1 + hm2)) * c_arr[..., -3]
                )
            flux = diffusivity * grad

        if c_arr.ndim == 1:
            return float(flux)
        return np.asarray(flux, dtype=float)

    def _eval_1d_current(
        self,
        flux: float | np.ndarray,
        n_electrons: int = 1,
        area: float | None = None,
    ) -> float | np.ndarray:
        """Compute Faradaic current I = n * F * A * J."""
        if n_electrons <= 0:
            raise ValueError(
                f"Number of electrons n_electrons must be positive (n_electrons > 0), got {n_electrons}."
            )
        active_area = self.area if area is None else float(area)
        scale = float(n_electrons) * FARADAY * active_area
        if np.isscalar(flux):
            return float(scale * float(flux))  # type: ignore[arg-type]
        flux_arr = np.asarray(flux, dtype=float)
        if flux_arr.ndim == 0:
            return float(scale * float(flux_arr.item()))
        return scale * flux_arr

area abstractmethod property

area: float

Active electrode surface area \(A\) in \(\text{m}^2\) or \(\text{cm}^2\).

surface_area property

surface_area: float

Active electrode surface area \(A\) (alias for area).

symmetry_factor property

symmetry_factor: int

Radial coordinate symmetry exponent \(p \in \{0, 1, 2\}\).

validate_grid

validate_grid(grid: ndarray | Mesh1D) -> np.ndarray

Validate that a spatial grid satisfies the coordinate domain constraints.

Parameters:

Name	Type	Description	Default
grid	ndarray or Mesh1D	1D array or mesh of spatial coordinates.	required

Returns:

Type	Description
ndarray	Validated 1D floating-point coordinate grid.

Raises:

Type	Description
ValueError	If grid is not 1D, has fewer than 3 points, is not strictly monotonically increasing, or violates radial domain bounds (\(r_0 > 0\)).

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

def validate_grid(self, grid: np.ndarray | Mesh1D) -> np.ndarray:
    """Validate that a spatial grid satisfies the coordinate domain constraints.

    Parameters
    ----------
    grid : np.ndarray or Mesh1D
        1D array or mesh of spatial coordinates.

    Returns
    -------
    np.ndarray
        Validated 1D floating-point coordinate grid.

    Raises
    ------
    ValueError
        If ``grid`` is not 1D, has fewer than 3 points, is not strictly
        monotonically increasing, or violates radial domain bounds (:math:`r_0 > 0`).
    """
    arr = np.asarray(grid.nodes if isinstance(grid, Mesh1D) else grid, dtype=float)
    if arr.ndim != 1 or len(arr) < 3:
        raise ValueError("grid must be a 1D array with at least 3 points.")
    if np.any(np.diff(arr) <= 0.0):
        raise ValueError("grid points must be strictly monotonically increasing.")
    if self.p > 0 and arr[0] <= 0.0:
        raise ValueError(
            f"{self.__class__.__name__} requires strictly positive radial coordinates "
            f"(r0 = grid[0] > 0), got grid[0]={arr[0]}."
        )
    return arr

laplacian abstractmethod

laplacian(c: ndarray, grid: ndarray | Mesh1D, bc_left: Any = None, bc_right: Any = None, t: float = 0.0) -> np.ndarray

Evaluate the spatial Laplacian \(\nabla^2 c\) at interior and boundary points.

Parameters:

Name	Type	Description	Default
c	ndarray	1D concentration profile array of shape (N,).	required
grid	ndarray or Mesh1D	1D spatial coordinate array or mesh of shape (N,).	required
bc_left	BoundaryCondition	Boundary condition at the inner/left boundary (\(i = 0\)).	None
bc_right	BoundaryCondition	Boundary condition at the outer/right boundary (\(i = N - 1\)).	None
t	float	Evaluation time in seconds for time-dependent boundary conditions.	0.0

Returns:

Type	Description
ndarray	1D array of \(\nabla^2 c\) values of shape (N,).

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

@abstractmethod
def laplacian(
    self,
    c: np.ndarray,
    grid: np.ndarray | Mesh1D,
    bc_left: Any = None,
    bc_right: Any = None,
    t: float = 0.0,
) -> np.ndarray:
    r"""Evaluate the spatial Laplacian :math:`\nabla^2 c` at interior and boundary points.

    Parameters
    ----------
    c : np.ndarray
        1D concentration profile array of shape ``(N,)``.
    grid : np.ndarray or Mesh1D
        1D spatial coordinate array or mesh of shape ``(N,)``.
    bc_left : BoundaryCondition, optional
        Boundary condition at the inner/left boundary (:math:`i = 0`).
    bc_right : BoundaryCondition, optional
        Boundary condition at the outer/right boundary (:math:`i = N - 1`).
    t : float, default 0.0
        Evaluation time in seconds for time-dependent boundary conditions.

    Returns
    -------
    np.ndarray
        1D array of :math:`\nabla^2 c` values of shape ``(N,)``.
    """
    raise NotImplementedError

laplacian_matrix abstractmethod

laplacian_matrix(grid: ndarray | Mesh1D, bc_left: Any = None, bc_right: Any = None) -> tuple[np.ndarray, np.ndarray]

Return banded matrix representation of \(\nabla^2\) for implicit solvers.

Constructs the tridiagonal finite-difference operator \(\mathbf{L}\) in scipy.linalg.solve_banded (1, 1) format of shape (3, N), together with a boundary gradient scaling vector \(\mathbf{b}_{\text{scale}}\) of shape (N,) such that for Neumann boundary gradients \(g_{\text{left}}\), \(g_{\text{right}}\):

\[ \nabla^2 \mathbf{c} = \mathbf{L}\mathbf{c} + \mathbf{b}_{\text{scale}} \odot [g_{\text{left}}, 0, \dots, 0, g_{\text{right}}]^T \]

Parameters:

Name	Type	Description	Default
grid	ndarray or Mesh1D	1D spatial coordinate array or mesh of shape (N,).	required
bc_left	BoundaryCondition	Boundary condition at the left boundary (\(i = 0\)).	None
bc_right	BoundaryCondition	Boundary condition at the right boundary (\(i = N - 1\)).	None

Returns:

Name	Type	Description
L_banded	ndarray	Banded tridiagonal matrix of shape (3, N) where row 0 holds the superdiagonal, row 1 holds the main diagonal, and row 2 holds the subdiagonal.
bc_scale	ndarray	1D boundary gradient multiplier array of shape (N,).

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

@abstractmethod
def laplacian_matrix(
    self,
    grid: np.ndarray | Mesh1D,
    bc_left: Any = None,
    bc_right: Any = None,
) -> tuple[np.ndarray, np.ndarray]:
    r"""Return banded matrix representation of :math:`\nabla^2` for implicit solvers.

    Constructs the tridiagonal finite-difference operator :math:`\mathbf{L}` in
    ``scipy.linalg.solve_banded`` ``(1, 1)`` format of shape ``(3, N)``, together
    with a boundary gradient scaling vector :math:`\mathbf{b}_{\text{scale}}` of
    shape ``(N,)`` such that for Neumann boundary gradients :math:`g_{\text{left}}`,
    :math:`g_{\text{right}}`:

    .. math::

        \nabla^2 \mathbf{c} = \mathbf{L}\mathbf{c}
            + \mathbf{b}_{\text{scale}} \odot [g_{\text{left}}, 0, \dots, 0, g_{\text{right}}]^T

    Parameters
    ----------
    grid : np.ndarray or Mesh1D
        1D spatial coordinate array or mesh of shape ``(N,)``.
    bc_left : BoundaryCondition, optional
        Boundary condition at the left boundary (:math:`i = 0`).
    bc_right : BoundaryCondition, optional
        Boundary condition at the right boundary (:math:`i = N - 1`).

    Returns
    -------
    L_banded : np.ndarray
        Banded tridiagonal matrix of shape ``(3, N)`` where row 0 holds the
        superdiagonal, row 1 holds the main diagonal, and row 2 holds the subdiagonal.
    bc_scale : np.ndarray
        1D boundary gradient multiplier array of shape ``(N,)``.
    """
    raise NotImplementedError

surface_flux abstractmethod

surface_flux(c: ndarray, grid: ndarray | Mesh1D, diffusivity: float, side: str = 'left') -> float | np.ndarray

Calculate normal interfacial diffusion flux \(J = -D (\partial c / \partial n)\) at a boundary.

Follows the electrode-normal convention where the unit normal \(\mathbf{n}\) points from the boundary into the solution domain (\(+\hat{x}\) or \(+\hat{r}\) at the left boundary, and \(-\hat{x}\) or \(-\hat{r}\) at the right boundary).

Parameters:

Name	Type	Description	Default
c	ndarray	Concentration profile of shape (N_x,) or time series of shape (N_t, N_x).	required
grid	ndarray or Mesh1D	1D spatial coordinate array or mesh of shape (N_x,).	required
diffusivity	float	Species diffusion coefficient \(D > 0\).	required
side	('left', 'right')	Domain boundary at which to evaluate the normal flux.	'left'

Returns:

Type	Description
float or ndarray	Interfacial flux \(J\) (scalar float if c is 1D, or 1D np.ndarray of shape (N_t,) if c is 2D).

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

@abstractmethod
def surface_flux(
    self,
    c: np.ndarray,
    grid: np.ndarray | Mesh1D,
    diffusivity: float,
    side: str = "left",
) -> float | np.ndarray:
    r"""Calculate normal interfacial diffusion flux :math:`J = -D (\partial c / \partial n)` at a boundary.

    Follows the electrode-normal convention where the unit normal :math:`\mathbf{n}`
    points from the boundary into the solution domain (:math:`+\hat{x}` or :math:`+\hat{r}`
    at the left boundary, and :math:`-\hat{x}` or :math:`-\hat{r}` at the right boundary).

    Parameters
    ----------
    c : np.ndarray
        Concentration profile of shape ``(N_x,)`` or time series of shape ``(N_t, N_x)``.
    grid : np.ndarray or Mesh1D
        1D spatial coordinate array or mesh of shape ``(N_x,)``.
    diffusivity : float
        Species diffusion coefficient :math:`D > 0`.
    side : {'left', 'right'}, default 'left'
        Domain boundary at which to evaluate the normal flux.

    Returns
    -------
    float or np.ndarray
        Interfacial flux :math:`J` (scalar ``float`` if ``c`` is 1D, or 1D
        ``np.ndarray`` of shape ``(N_t,)`` if ``c`` is 2D).
    """
    raise NotImplementedError

current abstractmethod

current(flux: float | ndarray, n_electrons: int = 1, side: Literal['left', 'right'] = 'left', grid: ndarray | Mesh1D | None = None) -> float | np.ndarray

Integrate normal flux over the active electrode surface to compute Faradaic current.

Computes the total Faradaic current \(I = n F A J\) according to IUPAC sign conventions (cathodic reduction \(I < 0\) when reactant flux into solution is negative, \(J < 0\)).

Parameters:

Name	Type	Description	Default
flux	float or ndarray	Normal interfacial diffusion flux \(J\) at the electrode surface.	required
n_electrons	int	Number of electrons transferred per reaction (\(n > 0\)).	1
side	('left', 'right')	Domain boundary at which to compute the Faradaic current.	'left'
grid	ndarray or Mesh1D	Spatial grid coordinates used to determine surface area if coordinate-dependent.	None

Returns:

Type	Description
float or ndarray	Faradaic current \(I\) in Amperes (\(\text{A}\)).

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

@abstractmethod
def current(
    self,
    flux: float | np.ndarray,
    n_electrons: int = 1,
    side: Literal["left", "right"] = "left",
    grid: np.ndarray | Mesh1D | None = None,
) -> float | np.ndarray:
    r"""Integrate normal flux over the active electrode surface to compute Faradaic current.

    Computes the total Faradaic current :math:`I = n F A J` according to IUPAC
    sign conventions (cathodic reduction :math:`I < 0` when reactant flux into
    solution is negative, :math:`J < 0`).

    Parameters
    ----------
    flux : float or np.ndarray
        Normal interfacial diffusion flux :math:`J` at the electrode surface.
    n_electrons : int, default 1
        Number of electrons transferred per reaction (:math:`n > 0`).
    side : {'left', 'right'}, default 'left'
        Domain boundary at which to compute the Faradaic current.
    grid : np.ndarray or Mesh1D, optional
        Spatial grid coordinates used to determine surface area if coordinate-dependent.

    Returns
    -------
    float or np.ndarray
        Faradaic current :math:`I` in Amperes (:math:`\text{A}`).
    """
    raise NotImplementedError

UnknownGeometryError

Bases: KeyError, ValueError

Exception raised when an unrecognized geometry identifier is requested.

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

class UnknownGeometryError(KeyError, ValueError):
    """Exception raised when an unrecognized geometry identifier is requested."""

    def __str__(self) -> str:
        return str(self.args[0]) if self.args else ""

CylindricalGeometry

Bases: Geometry

1D cylindrical radial coordinate geometry (\(p = 1\)).

Models axisymmetric radial diffusion to or from a cylindrical wire or carbon fiber microelectrode of radius \(r_0 > 0\) and length \(L > 0\) across radial domain \(r \in [r_0, r_{\max}]\):

\[ \nabla^2 c = \frac{\partial^2 c}{\partial r^2} + \frac{1}{r} \frac{\partial c}{\partial r} = \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial c}{\partial r} \right) \]

Neglecting end effects (\(L \gg r_0\)), the active lateral electrode surface area at \(r = r_0\) is:

\[ A = 2 \pi r_0 L, \quad I = n F A J \]

Parameters:

Name	Type	Description	Default
r0	float	Cylinder radius \(r_0 > 0\) in meters or cm. Defaults to 1.0 until bound to a radial coordinate grid.	None
length	float	Active cylinder length \(L > 0\) in meters or cm.	1.0
radius	float	Keyword alias for r0.	None
L	float	Keyword alias for length.	None
area	float	Explicit surface area override \(A > 0\). If None, computed analytically as \(A = 2 \pi r_0 L\).	None

Attributes:

Name	Type	Description
name	str	Geometry identifier ('cylindrical').
dimension	int	Spatial dimension (1).
coordinate_name	str	Primary coordinate symbol ('r').
p	int	Symmetry factor (1).
r0	float	Cylinder radius \(r_0 > 0\).
length	float	Cylinder length \(L > 0\).

Raises:

Type	Description
ValueError	If r0 <= 0, length <= 0, or area <= 0.

Examples:

>>> import numpy as np
>>> from softpotato.geometry import CylindricalGeometry
>>> geom = CylindricalGeometry(r0=5e-6, length=1e-2)
>>> r = np.linspace(1.0, 5.0, 201)
>>> c = np.log(r)  # Harmonic profile in cylindrical coordinates: ∇²(ln r) = 0
>>> lap = geom.laplacian(c, r)
>>> bool(np.max(np.abs(lap[1:-1])) < 1e-4)
True

Source code in src/softpotato/geometry/cylindrical.py

class CylindricalGeometry(Geometry, name="cylindrical"):
    r"""1D cylindrical radial coordinate geometry (:math:`p = 1`).

    Models axisymmetric radial diffusion to or from a cylindrical wire or carbon
    fiber microelectrode of radius :math:`r_0 > 0` and length :math:`L > 0` across
    radial domain :math:`r \in [r_0, r_{\max}]`:

    .. math::

        \nabla^2 c = \frac{\partial^2 c}{\partial r^2} + \frac{1}{r} \frac{\partial c}{\partial r}
                   = \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial c}{\partial r} \right)

    Neglecting end effects (:math:`L \gg r_0`), the active lateral electrode surface area
    at :math:`r = r_0` is:

    .. math::

        A = 2 \pi r_0 L, \quad I = n F A J

    Parameters
    ----------
    r0 : float, optional
        Cylinder radius :math:`r_0 > 0` in meters or cm. Defaults to ``1.0`` until
        bound to a radial coordinate grid.
    length : float, default 1.0
        Active cylinder length :math:`L > 0` in meters or cm.
    radius : float, optional
        Keyword alias for ``r0``.
    L : float, optional
        Keyword alias for ``length``.
    area : float, optional
        Explicit surface area override :math:`A > 0`. If None, computed analytically
        as :math:`A = 2 \pi r_0 L`.

    Attributes
    ----------
    name : str
        Geometry identifier (``'cylindrical'``).
    dimension : int
        Spatial dimension (``1``).
    coordinate_name : str
        Primary coordinate symbol (``'r'``).
    p : int
        Symmetry factor (``1``).
    r0 : float
        Cylinder radius :math:`r_0 > 0`.
    length : float
        Cylinder length :math:`L > 0`.

    Raises
    ------
    ValueError
        If ``r0 <= 0``, ``length <= 0``, or ``area <= 0``.

    Examples
    --------
    >>> import numpy as np
    >>> from softpotato.geometry import CylindricalGeometry
    >>> geom = CylindricalGeometry(r0=5e-6, length=1e-2)
    >>> r = np.linspace(1.0, 5.0, 201)
    >>> c = np.log(r)  # Harmonic profile in cylindrical coordinates: ∇²(ln r) = 0
    >>> lap = geom.laplacian(c, r)
    >>> bool(np.max(np.abs(lap[1:-1])) < 1e-4)
    True
    """

    name: str = "cylindrical"
    dimension: int = 1
    coordinate_name: str = "r"
    p: int = 1

    def __init__(
        self,
        r0: float | None = None,
        length: float | None = None,
        *,
        radius: float | None = None,
        L: float | None = None,
        area: float | None = None,
    ) -> None:
        if r0 is not None and radius is not None and r0 != radius:
            raise ValueError(
                "Cannot specify both 'r0' and 'radius' with different values."
            )
        if length is not None and L is not None and length != L:
            raise ValueError(
                "Cannot specify both 'length' and 'L' with different values."
            )

        resolved_r0 = radius if radius is not None else r0
        if resolved_r0 is not None:
            if resolved_r0 <= 0.0:
                raise ValueError(
                    f"Cylinder radius r0 must be positive (r0 > 0), got {resolved_r0}."
                )
            self.r0 = float(resolved_r0)
            self._r0_explicit = True
        else:
            self.r0 = 1.0
            self._r0_explicit = False

        resolved_length = (
            1.0 if (length is None and L is None) else (L if L is not None else length)
        )
        assert resolved_length is not None
        if resolved_length <= 0.0:
            raise ValueError(
                f"Cylinder length must be positive (length > 0), got {resolved_length}."
            )
        self.length = float(resolved_length)

        if area is not None:
            if area <= 0.0:
                raise ValueError(
                    f"Electrode area must be positive (area > 0), got {area}."
                )
            self._explicit_area: float | None = float(area)
        else:
            self._explicit_area = None

    @property
    def radius(self) -> float:
        r"""Cylinder radius :math:`r_0` (alias for :attr:`r0`)."""
        return self.r0

    @property
    def L(self) -> float:
        r"""Cylinder length :math:`L` (alias for :attr:`length`)."""
        return self.length

    @property
    def area(self) -> float:
        r"""Active lateral cylinder surface area :math:`A = 2\pi r_0 L`."""
        if self._explicit_area is not None:
            return self._explicit_area
        return 2.0 * math.pi * self.r0 * self.length

    def validate_grid(self, grid: np.ndarray | Mesh1D) -> np.ndarray:
        """Validate radial grid bounds (:math:`r_0 > 0`) and infer ``r0`` if not explicitly set."""
        arr = super().validate_grid(grid)
        if not self._r0_explicit:
            self.r0 = float(arr[0])
        return arr

    def laplacian(
        self,
        c: np.ndarray,
        grid: np.ndarray | Mesh1D,
        bc_left: Any = None,
        bc_right: Any = None,
        t: float = 0.0,
    ) -> np.ndarray:
        r"""Evaluate the cylindrical Laplacian :math:`\nabla^2 c = \frac{\partial^2 c}{\partial r^2} + \frac{1}{r}\frac{\partial c}{\partial r}`.

        Parameters
        ----------
        c : np.ndarray
            1D concentration array of shape ``(N,)``.
        grid : np.ndarray or Mesh1D
            1D radial coordinate grid :math:`r \in [r_0, r_{\max}]` of shape ``(N,)`` with :math:`r_0 > 0`.
        bc_left : BoundaryCondition, optional
            Inner boundary condition at :math:`r = r_0`.
        bc_right : BoundaryCondition, optional
            Outer boundary condition at :math:`r = r_{\max}`.
        t : float, default 0.0
            Time in seconds for evaluating time-dependent boundary conditions.

        Returns
        -------
        np.ndarray
            1D array of :math:`\nabla^2 c` values of shape ``(N,)``.
        """
        return self._eval_1d_laplacian(c, grid, bc_left=bc_left, bc_right=bc_right, t=t)

    def laplacian_matrix(
        self,
        grid: np.ndarray | Mesh1D,
        bc_left: Any = None,
        bc_right: Any = None,
    ) -> tuple[np.ndarray, np.ndarray]:
        r"""Assemble banded matrix representation of the cylindrical Laplacian :math:`\nabla^2`.

        Parameters
        ----------
        grid : np.ndarray or Mesh1D
            1D radial grid :math:`r` of shape ``(N,)`` with :math:`r_0 > 0`.
        bc_left : BoundaryCondition, optional
            Inner boundary condition at :math:`r = r_0`.
        bc_right : BoundaryCondition, optional
            Outer boundary condition at :math:`r = r_{\max}`.

        Returns
        -------
        L_banded : np.ndarray
            Banded tridiagonal operator matrix of shape ``(3, N)``.
        bc_scale : np.ndarray
            Boundary gradient scaling vector of shape ``(N,)``.
        """
        return self._build_1d_laplacian_matrix(grid, bc_left=bc_left, bc_right=bc_right)

    def surface_flux(
        self,
        c: np.ndarray,
        grid: np.ndarray | Mesh1D,
        diffusivity: float,
        side: str = "left",
    ) -> float | np.ndarray:
        r"""Calculate normal interfacial diffusion flux :math:`J = -D (\partial c / \partial n)`.

        Parameters
        ----------
        c : np.ndarray
            Concentration array of shape ``(N_r,)`` or ``(N_t, N_r)``.
        grid : np.ndarray or Mesh1D
            1D radial grid :math:`r` of shape ``(N_r,)`` with :math:`r_0 > 0`.
        diffusivity : float
            Species diffusion coefficient :math:`D > 0`.
        side : {'left', 'right'}, default 'left'
            Boundary side at which to evaluate the normal flux.

        Returns
        -------
        float or np.ndarray
            Normal interfacial diffusion flux :math:`J`.
        """
        return self._eval_1d_surface_flux(c, grid, diffusivity=diffusivity, side=side)

    def current(
        self,
        flux: float | np.ndarray,
        n_electrons: int = 1,
        side: Literal["left", "right"] = "left",
        grid: np.ndarray | Mesh1D | None = None,
    ) -> float | np.ndarray:
        r"""Integrate normal flux over the cylindrical electrode surface :math:`I = n F A J`.

        Parameters
        ----------
        flux : float or np.ndarray
            Interfacial diffusion flux :math:`J`.
        n_electrons : int, default 1
            Number of electrons transferred per reaction (:math:`n > 0`).
        side : {'left', 'right'}, default 'left'
            Domain boundary at which to compute the Faradaic current.
        grid : np.ndarray or Mesh1D, optional
            Spatial grid coordinates used to determine outer boundary area when ``side='right'``.

        Returns
        -------
        float or np.ndarray
            Faradaic current :math:`I` in Amperes (:math:`\text{A}`).
        """
        if side == "right" and grid is not None and self._explicit_area is None:
            r_right = float(np.asarray(grid)[-1])
            area = 2.0 * math.pi * r_right * self.length
            return self._eval_1d_current(flux, n_electrons=n_electrons, area=area)
        return self._eval_1d_current(flux, n_electrons=n_electrons)

radius property

radius: float

Cylinder radius \(r_0\) (alias for r0).

L property

L: float

Cylinder length \(L\) (alias for length).

area property

area: float

Active lateral cylinder surface area \(A = 2\pi r_0 L\).

validate_grid

validate_grid(grid: ndarray | Mesh1D) -> np.ndarray

Validate radial grid bounds (\(r_0 > 0\)) and infer r0 if not explicitly set.

Source code in src/softpotato/geometry/cylindrical.py

def validate_grid(self, grid: np.ndarray | Mesh1D) -> np.ndarray:
    """Validate radial grid bounds (:math:`r_0 > 0`) and infer ``r0`` if not explicitly set."""
    arr = super().validate_grid(grid)
    if not self._r0_explicit:
        self.r0 = float(arr[0])
    return arr

laplacian

laplacian(c: ndarray, grid: ndarray | Mesh1D, bc_left: Any = None, bc_right: Any = None, t: float = 0.0) -> np.ndarray

Evaluate the cylindrical Laplacian \(\nabla^2 c = \frac{\partial^2 c}{\partial r^2} + \frac{1}{r}\frac{\partial c}{\partial r}\).

Parameters:

Name	Type	Description	Default
c	ndarray	1D concentration array of shape (N,).	required
grid	ndarray or Mesh1D	1D radial coordinate grid \(r \in [r_0, r_{\max}]\) of shape (N,) with \(r_0 > 0\).	required
bc_left	BoundaryCondition	Inner boundary condition at \(r = r_0\).	None
bc_right	BoundaryCondition	Outer boundary condition at \(r = r_{\max}\).	None
t	float	Time in seconds for evaluating time-dependent boundary conditions.	0.0

Returns:

Type	Description
ndarray	1D array of \(\nabla^2 c\) values of shape (N,).

Source code in src/softpotato/geometry/cylindrical.py

def laplacian(
    self,
    c: np.ndarray,
    grid: np.ndarray | Mesh1D,
    bc_left: Any = None,
    bc_right: Any = None,
    t: float = 0.0,
) -> np.ndarray:
    r"""Evaluate the cylindrical Laplacian :math:`\nabla^2 c = \frac{\partial^2 c}{\partial r^2} + \frac{1}{r}\frac{\partial c}{\partial r}`.

    Parameters
    ----------
    c : np.ndarray
        1D concentration array of shape ``(N,)``.
    grid : np.ndarray or Mesh1D
        1D radial coordinate grid :math:`r \in [r_0, r_{\max}]` of shape ``(N,)`` with :math:`r_0 > 0`.
    bc_left : BoundaryCondition, optional
        Inner boundary condition at :math:`r = r_0`.
    bc_right : BoundaryCondition, optional
        Outer boundary condition at :math:`r = r_{\max}`.
    t : float, default 0.0
        Time in seconds for evaluating time-dependent boundary conditions.

    Returns
    -------
    np.ndarray
        1D array of :math:`\nabla^2 c` values of shape ``(N,)``.
    """
    return self._eval_1d_laplacian(c, grid, bc_left=bc_left, bc_right=bc_right, t=t)

laplacian_matrix

laplacian_matrix(grid: ndarray | Mesh1D, bc_left: Any = None, bc_right: Any = None) -> tuple[np.ndarray, np.ndarray]

Assemble banded matrix representation of the cylindrical Laplacian \(\nabla^2\).

Parameters:

Name	Type	Description	Default
grid	ndarray or Mesh1D	1D radial grid \(r\) of shape (N,) with \(r_0 > 0\).	required
bc_left	BoundaryCondition	Inner boundary condition at \(r = r_0\).	None
bc_right	BoundaryCondition	Outer boundary condition at \(r = r_{\max}\).	None

Returns:

Name	Type	Description
L_banded	ndarray	Banded tridiagonal operator matrix of shape (3, N).
bc_scale	ndarray	Boundary gradient scaling vector of shape (N,).

Source code in src/softpotato/geometry/cylindrical.py

def laplacian_matrix(
    self,
    grid: np.ndarray | Mesh1D,
    bc_left: Any = None,
    bc_right: Any = None,
) -> tuple[np.ndarray, np.ndarray]:
    r"""Assemble banded matrix representation of the cylindrical Laplacian :math:`\nabla^2`.

    Parameters
    ----------
    grid : np.ndarray or Mesh1D
        1D radial grid :math:`r` of shape ``(N,)`` with :math:`r_0 > 0`.
    bc_left : BoundaryCondition, optional
        Inner boundary condition at :math:`r = r_0`.
    bc_right : BoundaryCondition, optional
        Outer boundary condition at :math:`r = r_{\max}`.

    Returns
    -------
    L_banded : np.ndarray
        Banded tridiagonal operator matrix of shape ``(3, N)``.
    bc_scale : np.ndarray
        Boundary gradient scaling vector of shape ``(N,)``.
    """
    return self._build_1d_laplacian_matrix(grid, bc_left=bc_left, bc_right=bc_right)

surface_flux

surface_flux(c: ndarray, grid: ndarray | Mesh1D, diffusivity: float, side: str = 'left') -> float | np.ndarray

Calculate normal interfacial diffusion flux \(J = -D (\partial c / \partial n)\).

Parameters:

Name	Type	Description	Default
c	ndarray	Concentration array of shape (N_r,) or (N_t, N_r).	required
grid	ndarray or Mesh1D	1D radial grid \(r\) of shape (N_r,) with \(r_0 > 0\).	required
diffusivity	float	Species diffusion coefficient \(D > 0\).	required
side	('left', 'right')	Boundary side at which to evaluate the normal flux.	'left'

Returns:

Type	Description
float or ndarray	Normal interfacial diffusion flux \(J\).

Source code in src/softpotato/geometry/cylindrical.py

def surface_flux(
    self,
    c: np.ndarray,
    grid: np.ndarray | Mesh1D,
    diffusivity: float,
    side: str = "left",
) -> float | np.ndarray:
    r"""Calculate normal interfacial diffusion flux :math:`J = -D (\partial c / \partial n)`.

    Parameters
    ----------
    c : np.ndarray
        Concentration array of shape ``(N_r,)`` or ``(N_t, N_r)``.
    grid : np.ndarray or Mesh1D
        1D radial grid :math:`r` of shape ``(N_r,)`` with :math:`r_0 > 0`.
    diffusivity : float
        Species diffusion coefficient :math:`D > 0`.
    side : {'left', 'right'}, default 'left'
        Boundary side at which to evaluate the normal flux.

    Returns
    -------
    float or np.ndarray
        Normal interfacial diffusion flux :math:`J`.
    """
    return self._eval_1d_surface_flux(c, grid, diffusivity=diffusivity, side=side)

current

current(flux: float | ndarray, n_electrons: int = 1, side: Literal['left', 'right'] = 'left', grid: ndarray | Mesh1D | None = None) -> float | np.ndarray

Integrate normal flux over the cylindrical electrode surface \(I = n F A J\).

Parameters:

Name	Type	Description	Default
flux	float or ndarray	Interfacial diffusion flux \(J\).	required
n_electrons	int	Number of electrons transferred per reaction (\(n > 0\)).	1
side	('left', 'right')	Domain boundary at which to compute the Faradaic current.	'left'
grid	ndarray or Mesh1D	Spatial grid coordinates used to determine outer boundary area when side='right'.	None

Returns:

Type	Description
float or ndarray	Faradaic current \(I\) in Amperes (\(\text{A}\)).

Source code in src/softpotato/geometry/cylindrical.py

def current(
    self,
    flux: float | np.ndarray,
    n_electrons: int = 1,
    side: Literal["left", "right"] = "left",
    grid: np.ndarray | Mesh1D | None = None,
) -> float | np.ndarray:
    r"""Integrate normal flux over the cylindrical electrode surface :math:`I = n F A J`.

    Parameters
    ----------
    flux : float or np.ndarray
        Interfacial diffusion flux :math:`J`.
    n_electrons : int, default 1
        Number of electrons transferred per reaction (:math:`n > 0`).
    side : {'left', 'right'}, default 'left'
        Domain boundary at which to compute the Faradaic current.
    grid : np.ndarray or Mesh1D, optional
        Spatial grid coordinates used to determine outer boundary area when ``side='right'``.

    Returns
    -------
    float or np.ndarray
        Faradaic current :math:`I` in Amperes (:math:`\text{A}`).
    """
    if side == "right" and grid is not None and self._explicit_area is None:
        r_right = float(np.asarray(grid)[-1])
        area = 2.0 * math.pi * r_right * self.length
        return self._eval_1d_current(flux, n_electrons=n_electrons, area=area)
    return self._eval_1d_current(flux, n_electrons=n_electrons)

PlanarGeometry

Bases: Geometry

1D planar Cartesian coordinate geometry (\(p = 0\)).

Models semi-infinite or finite linear diffusion perpendicular to a flat macroelectrode along spatial coordinate \(x \in [x_0, x_{\max}]\):

\[ \nabla^2 c = \frac{\partial^2 c}{\partial x^2} \]

The active electrode surface area is constant (\(A\)), yielding Faradaic current:

\[ I = n F A J, \quad J_{\text{left}} = -D \left.\frac{\partial c}{\partial x}\right|_{x_0} \]

Parameters:

Name	Type	Description	Default
area	float	Active electrode surface area \(A\) in \(\text{m}^2\) or \(\text{cm}^2\) (\(A > 0\)).	1.0

Attributes:

Name	Type	Description
name	str	Geometry identifier ('planar').
dimension	int	Spatial dimension (1).
coordinate_name	str	Primary coordinate symbol ('x').
p	int	Symmetry factor (0).

Raises:

Type	Description
ValueError	If area <= 0.

Examples:

>>> import numpy as np
>>> from softpotato.geometry import PlanarGeometry
>>> geom = PlanarGeometry(area=0.5)
>>> x = np.linspace(0.0, 1.0, 21)
>>> c = x**2  # d²(x²)/dx² = 2.0 everywhere
>>> lap = geom.laplacian(c, x)
>>> bool(np.allclose(lap, 2.0))
True
>>> flux = geom.surface_flux(x, x, diffusivity=1e-5, side="left")
>>> round(flux, 8)
-1e-05

Source code in src/softpotato/geometry/planar.py

class PlanarGeometry(Geometry, name="planar"):
    r"""1D planar Cartesian coordinate geometry (:math:`p = 0`).

    Models semi-infinite or finite linear diffusion perpendicular to a flat macroelectrode
    along spatial coordinate :math:`x \in [x_0, x_{\max}]`:

    .. math::

        \nabla^2 c = \frac{\partial^2 c}{\partial x^2}

    The active electrode surface area is constant (:math:`A`), yielding Faradaic current:

    .. math::

        I = n F A J, \quad J_{\text{left}} = -D \left.\frac{\partial c}{\partial x}\right|_{x_0}

    Parameters
    ----------
    area : float, default 1.0
        Active electrode surface area :math:`A` in :math:`\text{m}^2` or :math:`\text{cm}^2`
        (:math:`A > 0`).

    Attributes
    ----------
    name : str
        Geometry identifier (``'planar'``).
    dimension : int
        Spatial dimension (``1``).
    coordinate_name : str
        Primary coordinate symbol (``'x'``).
    p : int
        Symmetry factor (``0``).

    Raises
    ------
    ValueError
        If ``area <= 0``.

    Examples
    --------
    >>> import numpy as np
    >>> from softpotato.geometry import PlanarGeometry
    >>> geom = PlanarGeometry(area=0.5)
    >>> x = np.linspace(0.0, 1.0, 21)
    >>> c = x**2  # d²(x²)/dx² = 2.0 everywhere
    >>> lap = geom.laplacian(c, x)
    >>> bool(np.allclose(lap, 2.0))
    True
    >>> flux = geom.surface_flux(x, x, diffusivity=1e-5, side="left")
    >>> round(flux, 8)
    -1e-05
    """

    name: str = "planar"
    dimension: int = 1
    coordinate_name: str = "x"
    p: int = 0

    def __init__(self, area: float = 1.0) -> None:
        if area <= 0.0:
            raise ValueError(f"Electrode area must be positive (area > 0), got {area}.")
        self._area = float(area)

    @property
    def area(self) -> float:
        r"""Active electrode surface area :math:`A`."""
        return self._area

    def laplacian(
        self,
        c: np.ndarray,
        grid: np.ndarray | Mesh1D,
        bc_left: Any = None,
        bc_right: Any = None,
        t: float = 0.0,
    ) -> np.ndarray:
        r"""Evaluate the 1D Cartesian Laplacian :math:`\nabla^2 c = \partial^2 c / \partial x^2`.

        Parameters
        ----------
        c : np.ndarray
            1D concentration array of shape ``(N,)``.
        grid : np.ndarray or Mesh1D
            1D spatial grid :math:`x` of shape ``(N,)``.
        bc_left : BoundaryCondition, optional
            Left boundary condition at :math:`x_0`.
        bc_right : BoundaryCondition, optional
            Right boundary condition at :math:`x_{N-1}`.
        t : float, default 0.0
            Time in seconds for evaluating time-dependent boundary conditions.

        Returns
        -------
        np.ndarray
            1D array of :math:`\partial^2 c / \partial x^2` values of shape ``(N,)``.
        """
        return self._eval_1d_laplacian(c, grid, bc_left=bc_left, bc_right=bc_right, t=t)

    def laplacian_matrix(
        self,
        grid: np.ndarray | Mesh1D,
        bc_left: Any = None,
        bc_right: Any = None,
    ) -> tuple[np.ndarray, np.ndarray]:
        r"""Assemble banded matrix representation of :math:`\partial^2 / \partial x^2`.

        Parameters
        ----------
        grid : np.ndarray or Mesh1D
            1D spatial grid :math:`x` of shape ``(N,)``.
        bc_left : BoundaryCondition, optional
            Left boundary condition at :math:`x_0`.
        bc_right : BoundaryCondition, optional
            Right boundary condition at :math:`x_{N-1}`.

        Returns
        -------
        L_banded : np.ndarray
            Banded tridiagonal operator matrix of shape ``(3, N)``.
        bc_scale : np.ndarray
            Boundary gradient scaling vector of shape ``(N,)``.
        """
        return self._build_1d_laplacian_matrix(grid, bc_left=bc_left, bc_right=bc_right)

    def surface_flux(
        self,
        c: np.ndarray,
        grid: np.ndarray | Mesh1D,
        diffusivity: float,
        side: str = "left",
    ) -> float | np.ndarray:
        r"""Calculate normal interfacial diffusion flux :math:`J = -D (\partial c / \partial n)`.

        Parameters
        ----------
        c : np.ndarray
            Concentration array of shape ``(N_x,)`` or ``(N_t, N_x)``.
        grid : np.ndarray or Mesh1D
            1D spatial grid :math:`x` of shape ``(N_x,)``.
        diffusivity : float
            Species diffusion coefficient :math:`D > 0`.
        side : {'left', 'right'}, default 'left'
            Boundary side at which to evaluate the normal flux.

        Returns
        -------
        float or np.ndarray
            Normal interfacial diffusion flux :math:`J`.
        """
        return self._eval_1d_surface_flux(c, grid, diffusivity=diffusivity, side=side)

    def current(
        self,
        flux: float | np.ndarray,
        n_electrons: int = 1,
        side: Literal["left", "right"] = "left",
        grid: np.ndarray | Mesh1D | None = None,
    ) -> float | np.ndarray:
        r"""Compute Faradaic current :math:`I = n F A J`.

        Parameters
        ----------
        flux : float or np.ndarray
            Interfacial diffusion flux :math:`J`.
        n_electrons : int, default 1
            Number of electrons transferred per reaction (:math:`n > 0`).
        side : {'left', 'right'}, default 'left'
            Domain boundary at which to compute the Faradaic current.
        grid : np.ndarray, optional
            Spatial grid coordinates (unused for planar geometry).

        Returns
        -------
        float or np.ndarray
            Faradaic current :math:`I` in Amperes (:math:`\text{A}`).
        """
        return self._eval_1d_current(flux, n_electrons=n_electrons)

area property

area: float

Active electrode surface area \(A\).

laplacian

laplacian(c: ndarray, grid: ndarray | Mesh1D, bc_left: Any = None, bc_right: Any = None, t: float = 0.0) -> np.ndarray

Evaluate the 1D Cartesian Laplacian \(\nabla^2 c = \partial^2 c / \partial x^2\).

Parameters:

Name	Type	Description	Default
c	ndarray	1D concentration array of shape (N,).	required
grid	ndarray or Mesh1D	1D spatial grid \(x\) of shape (N,).	required
bc_left	BoundaryCondition	Left boundary condition at \(x_0\).	None
bc_right	BoundaryCondition	Right boundary condition at \(x_{N-1}\).	None
t	float	Time in seconds for evaluating time-dependent boundary conditions.	0.0

Returns:

Type	Description
ndarray	1D array of \(\partial^2 c / \partial x^2\) values of shape (N,).

Source code in src/softpotato/geometry/planar.py

def laplacian(
    self,
    c: np.ndarray,
    grid: np.ndarray | Mesh1D,
    bc_left: Any = None,
    bc_right: Any = None,
    t: float = 0.0,
) -> np.ndarray:
    r"""Evaluate the 1D Cartesian Laplacian :math:`\nabla^2 c = \partial^2 c / \partial x^2`.

    Parameters
    ----------
    c : np.ndarray
        1D concentration array of shape ``(N,)``.
    grid : np.ndarray or Mesh1D
        1D spatial grid :math:`x` of shape ``(N,)``.
    bc_left : BoundaryCondition, optional
        Left boundary condition at :math:`x_0`.
    bc_right : BoundaryCondition, optional
        Right boundary condition at :math:`x_{N-1}`.
    t : float, default 0.0
        Time in seconds for evaluating time-dependent boundary conditions.

    Returns
    -------
    np.ndarray
        1D array of :math:`\partial^2 c / \partial x^2` values of shape ``(N,)``.
    """
    return self._eval_1d_laplacian(c, grid, bc_left=bc_left, bc_right=bc_right, t=t)

laplacian_matrix

laplacian_matrix(grid: ndarray | Mesh1D, bc_left: Any = None, bc_right: Any = None) -> tuple[np.ndarray, np.ndarray]

Assemble banded matrix representation of \(\partial^2 / \partial x^2\).

Parameters:

Name	Type	Description	Default
grid	ndarray or Mesh1D	1D spatial grid \(x\) of shape (N,).	required
bc_left	BoundaryCondition	Left boundary condition at \(x_0\).	None
bc_right	BoundaryCondition	Right boundary condition at \(x_{N-1}\).	None

Returns:

Name	Type	Description
L_banded	ndarray	Banded tridiagonal operator matrix of shape (3, N).
bc_scale	ndarray	Boundary gradient scaling vector of shape (N,).

Source code in src/softpotato/geometry/planar.py

def laplacian_matrix(
    self,
    grid: np.ndarray | Mesh1D,
    bc_left: Any = None,
    bc_right: Any = None,
) -> tuple[np.ndarray, np.ndarray]:
    r"""Assemble banded matrix representation of :math:`\partial^2 / \partial x^2`.

    Parameters
    ----------
    grid : np.ndarray or Mesh1D
        1D spatial grid :math:`x` of shape ``(N,)``.
    bc_left : BoundaryCondition, optional
        Left boundary condition at :math:`x_0`.
    bc_right : BoundaryCondition, optional
        Right boundary condition at :math:`x_{N-1}`.

    Returns
    -------
    L_banded : np.ndarray
        Banded tridiagonal operator matrix of shape ``(3, N)``.
    bc_scale : np.ndarray
        Boundary gradient scaling vector of shape ``(N,)``.
    """
    return self._build_1d_laplacian_matrix(grid, bc_left=bc_left, bc_right=bc_right)

surface_flux

surface_flux(c: ndarray, grid: ndarray | Mesh1D, diffusivity: float, side: str = 'left') -> float | np.ndarray

Calculate normal interfacial diffusion flux \(J = -D (\partial c / \partial n)\).

Parameters:

Name	Type	Description	Default
c	ndarray	Concentration array of shape (N_x,) or (N_t, N_x).	required
grid	ndarray or Mesh1D	1D spatial grid \(x\) of shape (N_x,).	required
diffusivity	float	Species diffusion coefficient \(D > 0\).	required
side	('left', 'right')	Boundary side at which to evaluate the normal flux.	'left'

Returns:

Type	Description
float or ndarray	Normal interfacial diffusion flux \(J\).

Source code in src/softpotato/geometry/planar.py

def surface_flux(
    self,
    c: np.ndarray,
    grid: np.ndarray | Mesh1D,
    diffusivity: float,
    side: str = "left",
) -> float | np.ndarray:
    r"""Calculate normal interfacial diffusion flux :math:`J = -D (\partial c / \partial n)`.

    Parameters
    ----------
    c : np.ndarray
        Concentration array of shape ``(N_x,)`` or ``(N_t, N_x)``.
    grid : np.ndarray or Mesh1D
        1D spatial grid :math:`x` of shape ``(N_x,)``.
    diffusivity : float
        Species diffusion coefficient :math:`D > 0`.
    side : {'left', 'right'}, default 'left'
        Boundary side at which to evaluate the normal flux.

    Returns
    -------
    float or np.ndarray
        Normal interfacial diffusion flux :math:`J`.
    """
    return self._eval_1d_surface_flux(c, grid, diffusivity=diffusivity, side=side)

current

current(flux: float | ndarray, n_electrons: int = 1, side: Literal['left', 'right'] = 'left', grid: ndarray | Mesh1D | None = None) -> float | np.ndarray

Compute Faradaic current \(I = n F A J\).

Parameters:

Name	Type	Description	Default
flux	float or ndarray	Interfacial diffusion flux \(J\).	required
n_electrons	int	Number of electrons transferred per reaction (\(n > 0\)).	1
side	('left', 'right')	Domain boundary at which to compute the Faradaic current.	'left'
grid	ndarray	Spatial grid coordinates (unused for planar geometry).	None

Returns:

Type	Description
float or ndarray	Faradaic current \(I\) in Amperes (\(\text{A}\)).

Source code in src/softpotato/geometry/planar.py

def current(
    self,
    flux: float | np.ndarray,
    n_electrons: int = 1,
    side: Literal["left", "right"] = "left",
    grid: np.ndarray | Mesh1D | None = None,
) -> float | np.ndarray:
    r"""Compute Faradaic current :math:`I = n F A J`.

    Parameters
    ----------
    flux : float or np.ndarray
        Interfacial diffusion flux :math:`J`.
    n_electrons : int, default 1
        Number of electrons transferred per reaction (:math:`n > 0`).
    side : {'left', 'right'}, default 'left'
        Domain boundary at which to compute the Faradaic current.
    grid : np.ndarray, optional
        Spatial grid coordinates (unused for planar geometry).

    Returns
    -------
    float or np.ndarray
        Faradaic current :math:`I` in Amperes (:math:`\text{A}`).
    """
    return self._eval_1d_current(flux, n_electrons=n_electrons)

SphericalGeometry

Bases: Geometry

1D spherical radial coordinate geometry (\(p = 2\)).

Models radial diffusion to or from a spherical microelectrode, hanging mercury drop electrode (HMDE), or hemispherical microelectrode on an insulating plane across radial domain \(r \in [r_0, r_{\max}]\) with \(r_0 > 0\):

\[ \nabla^2 c = \frac{\partial^2 c}{\partial r^2} + \frac{2}{r} \frac{\partial c}{\partial r} = \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 \frac{\partial c}{\partial r} \right) \]

The active electrode surface area at \(r = r_0\) is:

• Full sphere (hemisphere=False): \(A = 4 \pi r_0^2\)
• Hemisphere (hemisphere=True): \(A = 2 \pi r_0^2\)

Parameters:

Name	Type	Description	Default
r0	float	Electrode radius \(r_0 > 0\) in meters or cm. Defaults to 1.0 until bound to a radial coordinate grid.	None
hemisphere	bool	If True, computes surface area for a hemisphere (\(A = 2 \pi r_0^2\)) instead of a full sphere (\(A = 4 \pi r_0^2\)).	False
radius	float	Keyword alias for r0.	None
mode	('sphere', 'hemisphere')	Alternative specification of full sphere vs. hemisphere geometry.	'sphere'
area	float	Explicit surface area override \(A > 0\). If None, computed analytically from \(r_0\) and hemisphere.	None

Attributes:

Name	Type	Description
name	str	Geometry identifier ('spherical').
dimension	int	Spatial dimension (1).
coordinate_name	str	Primary coordinate symbol ('r').
p	int	Symmetry factor (2).
r0	float	Electrode radius \(r_0 > 0\).
hemisphere	bool	Whether the electrode is a hemisphere on an insulating substrate.

Raises:

Type	Description
ValueError	If r0 <= 0, area <= 0, or mode is invalid.

Examples:

>>> import numpy as np
>>> from softpotato.geometry import SphericalGeometry
>>> geom = SphericalGeometry(r0=1e-4)
>>> r = np.linspace(1.0, 5.0, 201)
>>> c = 1.0 / r  # Harmonic profile in spherical coordinates: ∇²(1/r) = 0
>>> lap = geom.laplacian(c, r)
>>> bool(np.max(np.abs(lap[1:-1])) < 1e-4)
True

Source code in src/softpotato/geometry/spherical.py

class SphericalGeometry(Geometry, name="spherical"):
    r"""1D spherical radial coordinate geometry (:math:`p = 2`).

    Models radial diffusion to or from a spherical microelectrode, hanging mercury
    drop electrode (HMDE), or hemispherical microelectrode on an insulating plane
    across radial domain :math:`r \in [r_0, r_{\max}]` with :math:`r_0 > 0`:

    .. math::

        \nabla^2 c = \frac{\partial^2 c}{\partial r^2} + \frac{2}{r} \frac{\partial c}{\partial r}
                   = \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 \frac{\partial c}{\partial r} \right)

    The active electrode surface area at :math:`r = r_0` is:

    - **Full sphere** (``hemisphere=False``): :math:`A = 4 \pi r_0^2`
    - **Hemisphere** (``hemisphere=True``): :math:`A = 2 \pi r_0^2`

    Parameters
    ----------
    r0 : float, optional
        Electrode radius :math:`r_0 > 0` in meters or cm. Defaults to ``1.0`` until
        bound to a radial coordinate grid.
    hemisphere : bool, default False
        If True, computes surface area for a hemisphere (:math:`A = 2 \pi r_0^2`)
        instead of a full sphere (:math:`A = 4 \pi r_0^2`).
    radius : float, optional
        Keyword alias for ``r0``.
    mode : {'sphere', 'hemisphere'}, optional
        Alternative specification of full sphere vs. hemisphere geometry.
    area : float, optional
        Explicit surface area override :math:`A > 0`. If None, computed analytically
        from :math:`r_0` and ``hemisphere``.

    Attributes
    ----------
    name : str
        Geometry identifier (``'spherical'``).
    dimension : int
        Spatial dimension (``1``).
    coordinate_name : str
        Primary coordinate symbol (``'r'``).
    p : int
        Symmetry factor (``2``).
    r0 : float
        Electrode radius :math:`r_0 > 0`.
    hemisphere : bool
        Whether the electrode is a hemisphere on an insulating substrate.

    Raises
    ------
    ValueError
        If ``r0 <= 0``, ``area <= 0``, or ``mode`` is invalid.

    Examples
    --------
    >>> import numpy as np
    >>> from softpotato.geometry import SphericalGeometry
    >>> geom = SphericalGeometry(r0=1e-4)
    >>> r = np.linspace(1.0, 5.0, 201)
    >>> c = 1.0 / r  # Harmonic profile in spherical coordinates: ∇²(1/r) = 0
    >>> lap = geom.laplacian(c, r)
    >>> bool(np.max(np.abs(lap[1:-1])) < 1e-4)
    True
    """

    name: str = "spherical"
    dimension: int = 1
    coordinate_name: str = "r"
    p: int = 2

    def __init__(
        self,
        r0: float | None = None,
        hemisphere: bool = False,
        *,
        radius: float | None = None,
        mode: Literal["sphere", "hemisphere"] | None = None,
        area: float | None = None,
    ) -> None:
        if r0 is not None and radius is not None and r0 != radius:
            raise ValueError(
                "Cannot specify both 'r0' and 'radius' with different values."
            )

        resolved_r0 = radius if radius is not None else r0
        if resolved_r0 is not None:
            if resolved_r0 <= 0.0:
                raise ValueError(
                    f"Electrode radius r0 must be positive (r0 > 0), got {resolved_r0}."
                )
            self.r0 = float(resolved_r0)
            self._r0_explicit = True
        else:
            self.r0 = 1.0
            self._r0_explicit = False

        if mode is not None:
            if mode not in ("sphere", "hemisphere"):
                raise ValueError(
                    f"mode must be 'sphere' or 'hemisphere', got {mode!r}."
                )
            self.hemisphere = (mode == "hemisphere") or bool(hemisphere)
        else:
            self.hemisphere = bool(hemisphere)

        if area is not None:
            if area <= 0.0:
                raise ValueError(
                    f"Electrode area must be positive (area > 0), got {area}."
                )
            self._explicit_area: float | None = float(area)
        else:
            self._explicit_area = None

    @property
    def radius(self) -> float:
        r"""Electrode radius :math:`r_0` (alias for :attr:`r0`)."""
        return self.r0

    @property
    def area(self) -> float:
        r"""Active electrode surface area :math:`A = 4\pi r_0^2` (or :math:`2\pi r_0^2` for hemisphere)."""
        if self._explicit_area is not None:
            return self._explicit_area
        factor = 2.0 * math.pi if self.hemisphere else 4.0 * math.pi
        return factor * (self.r0**2)

    def validate_grid(self, grid: np.ndarray | Mesh1D) -> np.ndarray:
        """Validate radial grid bounds (:math:`r_0 > 0`) and infer ``r0`` if not explicitly set."""
        arr = super().validate_grid(grid)
        if not self._r0_explicit:
            self.r0 = float(arr[0])
        return arr

    def laplacian(
        self,
        c: np.ndarray,
        grid: np.ndarray | Mesh1D,
        bc_left: Any = None,
        bc_right: Any = None,
        t: float = 0.0,
    ) -> np.ndarray:
        r"""Evaluate the spherical Laplacian :math:`\nabla^2 c = \frac{\partial^2 c}{\partial r^2} + \frac{2}{r}\frac{\partial c}{\partial r}`.

        Parameters
        ----------
        c : np.ndarray
            1D concentration array of shape ``(N,)``.
        grid : np.ndarray or Mesh1D
            1D radial coordinate grid :math:`r \in [r_0, r_{\max}]` of shape ``(N,)`` with :math:`r_0 > 0`.
        bc_left : BoundaryCondition, optional
            Inner boundary condition at :math:`r = r_0`.
        bc_right : BoundaryCondition, optional
            Outer boundary condition at :math:`r = r_{\max}`.
        t : float, default 0.0
            Time in seconds for evaluating time-dependent boundary conditions.

        Returns
        -------
        np.ndarray
            1D array of :math:`\nabla^2 c` values of shape ``(N,)``.
        """
        return self._eval_1d_laplacian(c, grid, bc_left=bc_left, bc_right=bc_right, t=t)

    def laplacian_matrix(
        self,
        grid: np.ndarray | Mesh1D,
        bc_left: Any = None,
        bc_right: Any = None,
    ) -> tuple[np.ndarray, np.ndarray]:
        r"""Assemble banded matrix representation of the spherical Laplacian :math:`\nabla^2`.

        Parameters
        ----------
        grid : np.ndarray or Mesh1D
            1D radial grid :math:`r` of shape ``(N,)`` with :math:`r_0 > 0`.
        bc_left : BoundaryCondition, optional
            Inner boundary condition at :math:`r = r_0`.
        bc_right : BoundaryCondition, optional
            Outer boundary condition at :math:`r = r_{\max}`.

        Returns
        -------
        L_banded : np.ndarray
            Banded tridiagonal operator matrix of shape ``(3, N)``.
        bc_scale : np.ndarray
            Boundary gradient scaling vector of shape ``(N,)``.
        """
        return self._build_1d_laplacian_matrix(grid, bc_left=bc_left, bc_right=bc_right)

    def surface_flux(
        self,
        c: np.ndarray,
        grid: np.ndarray | Mesh1D,
        diffusivity: float,
        side: str = "left",
    ) -> float | np.ndarray:
        r"""Calculate normal interfacial diffusion flux :math:`J = -D (\partial c / \partial n)`.

        Parameters
        ----------
        c : np.ndarray
            Concentration array of shape ``(N_r,)`` or ``(N_t, N_r)``.
        grid : np.ndarray or Mesh1D
            1D radial grid :math:`r` of shape ``(N_r,)`` with :math:`r_0 > 0`.
        diffusivity : float
            Species diffusion coefficient :math:`D > 0`.
        side : {'left', 'right'}, default 'left'
            Boundary side at which to evaluate the normal flux.

        Returns
        -------
        float or np.ndarray
            Normal interfacial diffusion flux :math:`J`.
        """
        return self._eval_1d_surface_flux(c, grid, diffusivity=diffusivity, side=side)

    def current(
        self,
        flux: float | np.ndarray,
        n_electrons: int = 1,
        side: Literal["left", "right"] = "left",
        grid: np.ndarray | Mesh1D | None = None,
    ) -> float | np.ndarray:
        r"""Integrate normal flux over the spherical or hemispherical surface :math:`I = n F A J`.

        Parameters
        ----------
        flux : float or np.ndarray
            Interfacial diffusion flux :math:`J`.
        n_electrons : int, default 1
            Number of electrons transferred per reaction (:math:`n > 0`).
        side : {'left', 'right'}, default 'left'
            Domain boundary at which to compute the Faradaic current.
        grid : np.ndarray or Mesh1D, optional
            Spatial grid coordinates used to determine outer boundary area when ``side='right'``.

        Returns
        -------
        float or np.ndarray
            Faradaic current :math:`I` in Amperes (:math:`\text{A}`).
        """
        if side == "right" and grid is not None and self._explicit_area is None:
            r_right = float(np.asarray(grid)[-1])
            factor = 2.0 * math.pi if self.hemisphere else 4.0 * math.pi
            area = factor * (r_right**2)
            return self._eval_1d_current(flux, n_electrons=n_electrons, area=area)
        return self._eval_1d_current(flux, n_electrons=n_electrons)

radius property

radius: float

Electrode radius \(r_0\) (alias for r0).

area property

area: float

Active electrode surface area \(A = 4\pi r_0^2\) (or \(2\pi r_0^2\) for hemisphere).

validate_grid

validate_grid(grid: ndarray | Mesh1D) -> np.ndarray

Validate radial grid bounds (\(r_0 > 0\)) and infer r0 if not explicitly set.

Source code in src/softpotato/geometry/spherical.py

def validate_grid(self, grid: np.ndarray | Mesh1D) -> np.ndarray:
    """Validate radial grid bounds (:math:`r_0 > 0`) and infer ``r0`` if not explicitly set."""
    arr = super().validate_grid(grid)
    if not self._r0_explicit:
        self.r0 = float(arr[0])
    return arr

laplacian

laplacian(c: ndarray, grid: ndarray | Mesh1D, bc_left: Any = None, bc_right: Any = None, t: float = 0.0) -> np.ndarray

Evaluate the spherical Laplacian \(\nabla^2 c = \frac{\partial^2 c}{\partial r^2} + \frac{2}{r}\frac{\partial c}{\partial r}\).

Parameters:

Name	Type	Description	Default
c	ndarray	1D concentration array of shape (N,).	required
grid	ndarray or Mesh1D	1D radial coordinate grid \(r \in [r_0, r_{\max}]\) of shape (N,) with \(r_0 > 0\).	required
bc_left	BoundaryCondition	Inner boundary condition at \(r = r_0\).	None
bc_right	BoundaryCondition	Outer boundary condition at \(r = r_{\max}\).	None
t	float	Time in seconds for evaluating time-dependent boundary conditions.	0.0

Returns:

Type	Description
ndarray	1D array of \(\nabla^2 c\) values of shape (N,).

Source code in src/softpotato/geometry/spherical.py

def laplacian(
    self,
    c: np.ndarray,
    grid: np.ndarray | Mesh1D,
    bc_left: Any = None,
    bc_right: Any = None,
    t: float = 0.0,
) -> np.ndarray:
    r"""Evaluate the spherical Laplacian :math:`\nabla^2 c = \frac{\partial^2 c}{\partial r^2} + \frac{2}{r}\frac{\partial c}{\partial r}`.

    Parameters
    ----------
    c : np.ndarray
        1D concentration array of shape ``(N,)``.
    grid : np.ndarray or Mesh1D
        1D radial coordinate grid :math:`r \in [r_0, r_{\max}]` of shape ``(N,)`` with :math:`r_0 > 0`.
    bc_left : BoundaryCondition, optional
        Inner boundary condition at :math:`r = r_0`.
    bc_right : BoundaryCondition, optional
        Outer boundary condition at :math:`r = r_{\max}`.
    t : float, default 0.0
        Time in seconds for evaluating time-dependent boundary conditions.

    Returns
    -------
    np.ndarray
        1D array of :math:`\nabla^2 c` values of shape ``(N,)``.
    """
    return self._eval_1d_laplacian(c, grid, bc_left=bc_left, bc_right=bc_right, t=t)

laplacian_matrix

laplacian_matrix(grid: ndarray | Mesh1D, bc_left: Any = None, bc_right: Any = None) -> tuple[np.ndarray, np.ndarray]

Assemble banded matrix representation of the spherical Laplacian \(\nabla^2\).

Parameters:

Name	Type	Description	Default
grid	ndarray or Mesh1D	1D radial grid \(r\) of shape (N,) with \(r_0 > 0\).	required
bc_left	BoundaryCondition	Inner boundary condition at \(r = r_0\).	None
bc_right	BoundaryCondition	Outer boundary condition at \(r = r_{\max}\).	None

Returns:

Name	Type	Description
L_banded	ndarray	Banded tridiagonal operator matrix of shape (3, N).
bc_scale	ndarray	Boundary gradient scaling vector of shape (N,).

Source code in src/softpotato/geometry/spherical.py

def laplacian_matrix(
    self,
    grid: np.ndarray | Mesh1D,
    bc_left: Any = None,
    bc_right: Any = None,
) -> tuple[np.ndarray, np.ndarray]:
    r"""Assemble banded matrix representation of the spherical Laplacian :math:`\nabla^2`.

    Parameters
    ----------
    grid : np.ndarray or Mesh1D
        1D radial grid :math:`r` of shape ``(N,)`` with :math:`r_0 > 0`.
    bc_left : BoundaryCondition, optional
        Inner boundary condition at :math:`r = r_0`.
    bc_right : BoundaryCondition, optional
        Outer boundary condition at :math:`r = r_{\max}`.

    Returns
    -------
    L_banded : np.ndarray
        Banded tridiagonal operator matrix of shape ``(3, N)``.
    bc_scale : np.ndarray
        Boundary gradient scaling vector of shape ``(N,)``.
    """
    return self._build_1d_laplacian_matrix(grid, bc_left=bc_left, bc_right=bc_right)

surface_flux

surface_flux(c: ndarray, grid: ndarray | Mesh1D, diffusivity: float, side: str = 'left') -> float | np.ndarray

Calculate normal interfacial diffusion flux \(J = -D (\partial c / \partial n)\).

Parameters:

Name	Type	Description	Default
c	ndarray	Concentration array of shape (N_r,) or (N_t, N_r).	required
grid	ndarray or Mesh1D	1D radial grid \(r\) of shape (N_r,) with \(r_0 > 0\).	required
diffusivity	float	Species diffusion coefficient \(D > 0\).	required
side	('left', 'right')	Boundary side at which to evaluate the normal flux.	'left'

Returns:

Type	Description
float or ndarray	Normal interfacial diffusion flux \(J\).

Source code in src/softpotato/geometry/spherical.py

def surface_flux(
    self,
    c: np.ndarray,
    grid: np.ndarray | Mesh1D,
    diffusivity: float,
    side: str = "left",
) -> float | np.ndarray:
    r"""Calculate normal interfacial diffusion flux :math:`J = -D (\partial c / \partial n)`.

    Parameters
    ----------
    c : np.ndarray
        Concentration array of shape ``(N_r,)`` or ``(N_t, N_r)``.
    grid : np.ndarray or Mesh1D
        1D radial grid :math:`r` of shape ``(N_r,)`` with :math:`r_0 > 0`.
    diffusivity : float
        Species diffusion coefficient :math:`D > 0`.
    side : {'left', 'right'}, default 'left'
        Boundary side at which to evaluate the normal flux.

    Returns
    -------
    float or np.ndarray
        Normal interfacial diffusion flux :math:`J`.
    """
    return self._eval_1d_surface_flux(c, grid, diffusivity=diffusivity, side=side)

current

current(flux: float | ndarray, n_electrons: int = 1, side: Literal['left', 'right'] = 'left', grid: ndarray | Mesh1D | None = None) -> float | np.ndarray

Integrate normal flux over the spherical or hemispherical surface \(I = n F A J\).

Parameters:

Name	Type	Description	Default
flux	float or ndarray	Interfacial diffusion flux \(J\).	required
n_electrons	int	Number of electrons transferred per reaction (\(n > 0\)).	1
side	('left', 'right')	Domain boundary at which to compute the Faradaic current.	'left'
grid	ndarray or Mesh1D	Spatial grid coordinates used to determine outer boundary area when side='right'.	None

Returns:

Type	Description
float or ndarray	Faradaic current \(I\) in Amperes (\(\text{A}\)).

Source code in src/softpotato/geometry/spherical.py

def current(
    self,
    flux: float | np.ndarray,
    n_electrons: int = 1,
    side: Literal["left", "right"] = "left",
    grid: np.ndarray | Mesh1D | None = None,
) -> float | np.ndarray:
    r"""Integrate normal flux over the spherical or hemispherical surface :math:`I = n F A J`.

    Parameters
    ----------
    flux : float or np.ndarray
        Interfacial diffusion flux :math:`J`.
    n_electrons : int, default 1
        Number of electrons transferred per reaction (:math:`n > 0`).
    side : {'left', 'right'}, default 'left'
        Domain boundary at which to compute the Faradaic current.
    grid : np.ndarray or Mesh1D, optional
        Spatial grid coordinates used to determine outer boundary area when ``side='right'``.

    Returns
    -------
    float or np.ndarray
        Faradaic current :math:`I` in Amperes (:math:`\text{A}`).
    """
    if side == "right" and grid is not None and self._explicit_area is None:
        r_right = float(np.asarray(grid)[-1])
        factor = 2.0 * math.pi if self.hemisphere else 4.0 * math.pi
        area = factor * (r_right**2)
        return self._eval_1d_current(flux, n_electrons=n_electrons, area=area)
    return self._eval_1d_current(flux, n_electrons=n_electrons)

get_geometry

get_geometry(name_or_instance: str | Geometry = 'planar', **kwargs: Any) -> Geometry

Instantiate or resolve a spatial coordinate geometry by name or instance.

Parameters:

Name	Type	Description	Default
name_or_instance	str or Geometry	Registered geometry identifier (case-insensitive) such as 'planar', 'spherical', 'hemispherical', or 'cylindrical', or an already-instantiated Geometry object.	'planar'
**kwargs	Any	Configuration parameters forwarded to the geometry constructor: For 'planar' / 'cartesian': area (float, default 1.0): Electrode surface area \(A\). For 'spherical' / 'sphere' / 'hemispherical' / 'hemisphere': r0 or radius (float, default 1.0): Electrode radius \(r_0 > 0\). hemisphere (bool, default False): If True, uses hemispherical surface area \(A = 2\pi r_0^2\) instead of full sphere \(A = 4\pi r_0^2\). area (float, optional): Explicit electrode surface area override. For 'cylindrical' / 'cylinder': r0 or radius (float, default 1.0): Cylinder radius \(r_0 > 0\). length or L (float, default 1.0): Cylinder length \(L > 0\). area (float, optional): Explicit electrode surface area override.	{}

Returns:

Type	Description
Geometry	Resolved Geometry instance.

Raises:

Type	Description
UnknownGeometryError	If name_or_instance is an unrecognized string identifier (inherits from both KeyError and ValueError).
TypeError	If name_or_instance is neither a string nor a Geometry.
ValueError	If keyword arguments are supplied when name_or_instance is already a Geometry instance.

Examples:

>>> from softpotato.geometry import get_geometry
>>> planar = get_geometry("planar", area=0.05)
>>> planar.area
0.05
>>> sphere = get_geometry("spherical", r0=1e-4)
>>> sphere.coordinate_name
'r'
>>> get_geometry(planar) is planar
True

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

def get_geometry(
    name_or_instance: str | Geometry = "planar",
    **kwargs: Any,
) -> Geometry:
    """Instantiate or resolve a spatial coordinate geometry by name or instance.

    Parameters
    ----------
    name_or_instance : str or Geometry, default 'planar'
        Registered geometry identifier (case-insensitive) such as ``'planar'``,
        ``'spherical'``, ``'hemispherical'``, or ``'cylindrical'``, or an
        already-instantiated :class:`~softpotato.geometry.Geometry` object.
    **kwargs : Any
        Configuration parameters forwarded to the geometry constructor:

        - For ``'planar'`` / ``'cartesian'``:
          - ``area`` (float, default 1.0): Electrode surface area :math:`A`.
        - For ``'spherical'`` / ``'sphere'`` / ``'hemispherical'`` / ``'hemisphere'``:
          - ``r0`` or ``radius`` (float, default 1.0): Electrode radius :math:`r_0 > 0`.
          - ``hemisphere`` (bool, default False): If True, uses hemispherical surface
            area :math:`A = 2\\pi r_0^2` instead of full sphere :math:`A = 4\\pi r_0^2`.
          - ``area`` (float, optional): Explicit electrode surface area override.
        - For ``'cylindrical'`` / ``'cylinder'``:
          - ``r0`` or ``radius`` (float, default 1.0): Cylinder radius :math:`r_0 > 0`.
          - ``length`` or ``L`` (float, default 1.0): Cylinder length :math:`L > 0`.
          - ``area`` (float, optional): Explicit electrode surface area override.

    Returns
    -------
    Geometry
        Resolved :class:`~softpotato.geometry.Geometry` instance.

    Raises
    ------
    UnknownGeometryError
        If ``name_or_instance`` is an unrecognized string identifier (inherits from
        both :exc:`KeyError` and :exc:`ValueError`).
    TypeError
        If ``name_or_instance`` is neither a string nor a :class:`~softpotato.geometry.Geometry`.
    ValueError
        If keyword arguments are supplied when ``name_or_instance`` is already a
        :class:`~softpotato.geometry.Geometry` instance.

    Examples
    --------
    >>> from softpotato.geometry import get_geometry
    >>> planar = get_geometry("planar", area=0.05)
    >>> planar.area
    0.05
    >>> sphere = get_geometry("spherical", r0=1e-4)
    >>> sphere.coordinate_name
    'r'
    >>> get_geometry(planar) is planar
    True
    """
    if isinstance(name_or_instance, Geometry):
        if kwargs:
            raise ValueError(
                "Cannot pass constructor keyword arguments when name_or_instance is already a Geometry instance."
            )
        return name_or_instance

    if not isinstance(name_or_instance, str):
        raise TypeError(
            f"name_or_instance must be a string or Geometry instance, got {type(name_or_instance).__name__}."
        )

    normalized = name_or_instance.strip().lower()
    if normalized in ("hemisphere", "hemispherical"):
        kwargs.setdefault("hemisphere", True)
        normalized = "spherical"

    if normalized not in Geometry._registry:
        available = ", ".join(sorted(Geometry._registry.keys()))
        raise UnknownGeometryError(
            f"Unknown geometry '{name_or_instance}'. Available geometries: {available}"
        )
    return Geometry._registry[normalized](**kwargs)

list_geometries

list_geometries() -> list[str]

Return a sorted list of all registered geometry names and aliases.

Returns:

Type	Description
list[str]	Alphabetically sorted list of registered geometry identifiers.

Examples:

>>> from softpotato.geometry import list_geometries
>>> geoms = list_geometries()
>>> "planar" in geoms and "spherical" in geoms and "cylindrical" in geoms
True

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

def list_geometries() -> list[str]:
    """Return a sorted list of all registered geometry names and aliases.

    Returns
    -------
    list[str]
        Alphabetically sorted list of registered geometry identifiers.

    Examples
    --------
    >>> from softpotato.geometry import list_geometries
    >>> geoms = list_geometries()
    >>> "planar" in geoms and "spherical" in geoms and "cylindrical" in geoms
    True
    """
    return sorted(Geometry._registry.keys())

register_geometry

register_geometry(name: str) -> Callable[[type[Geometry]], type[Geometry]]

Decorator to register a custom geometry class under a given alias.

Parameters:

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

Returns:

Type	Description
Callable[[type[Geometry]], type[Geometry]]	Class decorator registering the subclass in Geometry._registry.

Examples:

>>> from softpotato.geometry import Geometry, register_geometry, get_geometry
>>> @register_geometry("custom_planar_demo")
... class CustomPlanar(Geometry):
...     @property
...     def area(self) -> float:
...         return 2.0
...     def laplacian(self, c, grid, bc_left=None, bc_right=None, t=0.0):
...         return self._eval_1d_laplacian(c, grid, bc_left, bc_right, t)
...     def laplacian_matrix(self, grid, bc_left=None, bc_right=None):
...         return self._build_1d_laplacian_matrix(grid, bc_left, bc_right)
...     def surface_flux(self, c, grid, diffusivity, side="left"):
...         return self._eval_1d_surface_flux(c, grid, diffusivity, side)
...     def current(self, flux, n_electrons=1):
...         return self._eval_1d_current(flux, n_electrons)
>>> geom = get_geometry("custom_planar_demo")
>>> geom.area
2.0

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

def register_geometry(name: str) -> Callable[[type[Geometry]], type[Geometry]]:
    """Decorator to register a custom geometry class under a given alias.

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

    Returns
    -------
    Callable[[type[Geometry]], type[Geometry]]
        Class decorator registering the subclass in ``Geometry._registry``.

    Examples
    --------
    >>> from softpotato.geometry import Geometry, register_geometry, get_geometry
    >>> @register_geometry("custom_planar_demo")
    ... class CustomPlanar(Geometry):
    ...     @property
    ...     def area(self) -> float:
    ...         return 2.0
    ...     def laplacian(self, c, grid, bc_left=None, bc_right=None, t=0.0):
    ...         return self._eval_1d_laplacian(c, grid, bc_left, bc_right, t)
    ...     def laplacian_matrix(self, grid, bc_left=None, bc_right=None):
    ...         return self._build_1d_laplacian_matrix(grid, bc_left, bc_right)
    ...     def surface_flux(self, c, grid, diffusivity, side="left"):
    ...         return self._eval_1d_surface_flux(c, grid, diffusivity, side)
    ...     def current(self, flux, n_electrons=1):
    ...         return self._eval_1d_current(flux, n_electrons)
    >>> geom = get_geometry("custom_planar_demo")
    >>> geom.area
    2.0
    """

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

    return decorator