Simulation Benchmarking Suite (softpotato.analytical.benchmark)

The softpotato.analytical.benchmark module provides automated validation harnesses, exact analytical benchmarks, and convergence testing utilities.

benchmark

Automated numerical solver validation suite and analytical benchmarking harness.

This submodule provides a comprehensive benchmarking suite that validates numerical PDE solvers (softpotato.solver) against closed-form analytical ground truths (softpotato.analytical).

Features

1. Standardized Error Metrics:

   • Pointwise absolute error (\(L_\infty\) norm, compute_error_metrics).
   • Root-Mean-Square Error (RMSE).
   • Discrete relative \(L_2\) error norm.
   • Mean absolute error (MAE).
   • Mean and maximum relative errors with zero-division protection.

2. Order-of-Convergence Verification:

   • Spatial discretization convergence verification (\(\mathcal{O}(\Delta x^2)\)).
   • Temporal discretization convergence verification (\(\mathcal{O}(\Delta t^2)\) for Crank–Nicolson, \(\mathcal{O}(\Delta t)\) for backward-Euler BTCS).
   • Empirical order estimation via log-log regression (estimate_convergence_order).

3. Canonical Electrochemical Benchmark Problems:

   • Fourier sine diffusion decay (run_fourier_decay_benchmark).
   • Planar Cottrell potential step chronoamperometry (run_cottrell_benchmark).
   • Total mass conservation under zero-flux Neumann BCs (run_mass_conservation_benchmark).

4. Automated Suite Execution:

   • Unified single-case runner (run_benchmark).
   • Multi-solver comparative benchmark suite (run_benchmark_suite) with formatted tabular ASCII reporting.

Functions & Classes

Name	Type	Description
compute_error_metrics	Function	Compute standardized error norms (\(L_\infty\), RMSE, relative \(L_2\), MAE) between arrays.
estimate_convergence_order	Function	Estimate empirical numerical order of accuracy \(p\) via log-log linear regression.
verify_spatial_convergence	Function	Verify spatial discretization convergence order across a sequence of grid refinements.
verify_temporal_convergence	Function	Verify temporal discretization convergence order across a sequence of time step refinements.
run_fourier_decay_benchmark	Function	Benchmark a solver against the exact 1D Fourier sine-mode diffusion decay solution.
run_cottrell_benchmark	Function	Benchmark a solver against the exact planar Cottrell concentration profile and surface flux.
run_spherical_cottrell_benchmark	Function	Benchmark a solver against the exact spherical Cottrell transient current and profile.
run_mass_conservation_benchmark	Function	Benchmark total mass conservation in an insulated domain with zero-flux Neumann boundaries.
run_grid_efficiency_benchmark	Function	Evaluate runtime, memory scaling, and \(L_2\) error norms across uniform vs. expanding meshes.
list_benchmark_cases	Function	Return the list of registered canonical benchmark case identifiers.
run_benchmark	Function	Run a single named benchmark problem for a specified solver.
run_benchmark_suite	Function	Execute the full validation suite across one or more solvers and generate a summary report.
BenchmarkMetrics	Dataclass	Container for pointwise and integral numerical error metrics.
ConvergenceResult	Dataclass	Container for grid/time-step refinement errors and observed convergence order.
BenchmarkCaseResult	Dataclass	Container for the outcome, timing, and error metrics of a single benchmark case.
BenchmarkSuiteResult	Dataclass	Container aggregating results across multiple solvers and benchmark cases.
GridComparisonEntry	Dataclass	Container for a single mesh configuration in a uniform vs. expanding grid benchmark.
GridComparisonResult	Dataclass	Container comparing accuracy, runtime speedup, and memory scaling across meshes.

BenchmarkMetrics dataclass

Standardized error metrics between numerical and analytical solutions.

Parameters:

Name	Type	Description	Default
rmse	float	Root-mean-square error: \[ \text{RMSE} = \sqrt{\frac{1}{N} \sum_{i=1}^N (u_i^{\text{num}} - u_i^{\text{exact}})^2} \]	required
max_abs_error	float	Maximum pointwise absolute deviation (\(L_\infty\) norm): \[ \|u^{\text{num}} - u^{\text{exact}}\|_\infty = \max_i |u_i^{\text{num}} - u_i^{\text{exact}}| \]	required
l2_error	float	Discrete relative \(L_2\) error norm: \[ L_2 = \frac{\|u^{\text{num}} - u^{\text{exact}}\|_2}{\|u^{\text{exact}}\|_2} \] (evaluates to absolute \(L_2\) norm if \(\|u^{\text{exact}}\|_2 = 0\)).	required
mean_abs_error	float	Mean absolute error (MAE): \[ \text{MAE} = \frac{1}{N} \sum_{i=1}^N |u_i^{\text{num}} - u_i^{\text{exact}}| \]	required
mean_rel_error	float	Mean relative error across all evaluated points: \[ \text{MRE} = \frac{1}{N} \sum_{i=1}^N \frac{|u_i^{\text{num}} - u_i^{\text{exact}}|}{|u_i^{\text{exact}}| + \epsilon} \]	required
max_rel_error	float	Maximum relative error: \[ \text{MaxRE} = \max_i \frac{|u_i^{\text{num}} - u_i^{\text{exact}}|}{|u_i^{\text{exact}}| + \epsilon} \]	required

Examples:

>>> import numpy as np
>>> from softpotato.analytical.benchmark import compute_error_metrics
>>> num = np.array([1.0, 2.01, 3.0])
>>> exact = np.array([1.0, 2.0, 3.0])
>>> metrics = compute_error_metrics(num, exact)
>>> round(metrics.max_abs_error, 4)
0.01

Source code in src/softpotato/analytical/benchmark.py

@dataclass(frozen=True, kw_only=True)
class BenchmarkMetrics:
    r"""Standardized error metrics between numerical and analytical solutions.

    Parameters
    ----------
    rmse : float
        Root-mean-square error:

        .. math::

            \text{RMSE} = \sqrt{\frac{1}{N} \sum_{i=1}^N (u_i^{\text{num}} - u_i^{\text{exact}})^2}

    max_abs_error : float
        Maximum pointwise absolute deviation (:math:`L_\infty` norm):

        .. math::

            \|u^{\text{num}} - u^{\text{exact}}\|_\infty = \max_i |u_i^{\text{num}} - u_i^{\text{exact}}|

    l2_error : float
        Discrete relative :math:`L_2` error norm:

        .. math::

            L_2 = \frac{\|u^{\text{num}} - u^{\text{exact}}\|_2}{\|u^{\text{exact}}\|_2}

        (evaluates to absolute :math:`L_2` norm if :math:`\|u^{\text{exact}}\|_2 = 0`).
    mean_abs_error : float
        Mean absolute error (MAE):

        .. math::

            \text{MAE} = \frac{1}{N} \sum_{i=1}^N |u_i^{\text{num}} - u_i^{\text{exact}}|

    mean_rel_error : float
        Mean relative error across all evaluated points:

        .. math::

            \text{MRE} = \frac{1}{N} \sum_{i=1}^N \frac{|u_i^{\text{num}} - u_i^{\text{exact}}|}{|u_i^{\text{exact}}| + \epsilon}

    max_rel_error : float
        Maximum relative error:

        .. math::

            \text{MaxRE} = \max_i \frac{|u_i^{\text{num}} - u_i^{\text{exact}}|}{|u_i^{\text{exact}}| + \epsilon}

    Examples
    --------
    >>> import numpy as np
    >>> from softpotato.analytical.benchmark import compute_error_metrics
    >>> num = np.array([1.0, 2.01, 3.0])
    >>> exact = np.array([1.0, 2.0, 3.0])
    >>> metrics = compute_error_metrics(num, exact)
    >>> round(metrics.max_abs_error, 4)
    0.01
    """

    rmse: float
    max_abs_error: float
    l2_error: float
    mean_abs_error: float
    mean_rel_error: float
    max_rel_error: float

    def to_dict(self) -> dict[str, float]:
        """Convert metrics to a standard Python dictionary."""
        return {
            "rmse": self.rmse,
            "max_abs_error": self.max_abs_error,
            "l2_error": self.l2_error,
            "mean_abs_error": self.mean_abs_error,
            "mean_rel_error": self.mean_rel_error,
            "max_rel_error": self.max_rel_error,
        }

to_dict

to_dict() -> dict[str, float]

Convert metrics to a standard Python dictionary.

Source code in src/softpotato/analytical/benchmark.py

def to_dict(self) -> dict[str, float]:
    """Convert metrics to a standard Python dictionary."""
    return {
        "rmse": self.rmse,
        "max_abs_error": self.max_abs_error,
        "l2_error": self.l2_error,
        "mean_abs_error": self.mean_abs_error,
        "mean_rel_error": self.mean_rel_error,
        "max_rel_error": self.max_rel_error,
    }

ConvergenceResult dataclass

Results from an order-of-convergence refinement study.

Parameters:

Name	Type	Description	Default
spacings	ndarray	Discretization spacings tested (\(\Delta x\) for spatial, \(\Delta t\) for temporal).	required
errors	ndarray	Corresponding error norm evaluated for each spacing.	required
observed_order	float	Empirical convergence order \(p\) extracted from linear regression: \[ \ln(\text{error}) = p \cdot \ln(h) + C \]	required
expected_order	float	Theoretical convergence order (e.g., \(2.0\) for second-order schemes).	required
r_squared	float	Goodness-of-fit \(R^2\) of the log-log linear regression.	required
passed	bool	Whether observed_order meets expected_order within tolerance and r_squared is high (\(R^2 \ge 0.95\)).	required
spacing_type	('spatial', 'temporal')	Discretization domain refined during the study.	'spatial'
details	dict[str, Any]	Additional context such as solver name, test case, and tolerance thresholds.	dict()

Examples:

>>> import numpy as np
>>> from softpotato.analytical.benchmark import estimate_convergence_order
>>> h = np.array([0.1, 0.05, 0.025])
>>> err = 0.5 * h**2  # pure 2nd-order error
>>> p, r2 = estimate_convergence_order(h, err)
>>> round(p, 2)
2.0
>>> round(r2, 4)
1.0

Source code in src/softpotato/analytical/benchmark.py

@dataclass(frozen=True, kw_only=True)
class ConvergenceResult:
    r"""Results from an order-of-convergence refinement study.

    Parameters
    ----------
    spacings : np.ndarray
        Discretization spacings tested (:math:`\Delta x` for spatial,
        :math:`\Delta t` for temporal).
    errors : np.ndarray
        Corresponding error norm evaluated for each spacing.
    observed_order : float
        Empirical convergence order :math:`p` extracted from linear regression:

        .. math::

            \ln(\text{error}) = p \cdot \ln(h) + C

    expected_order : float
        Theoretical convergence order (e.g., :math:`2.0` for second-order schemes).
    r_squared : float
        Goodness-of-fit :math:`R^2` of the log-log linear regression.
    passed : bool
        Whether :attr:`observed_order` meets :attr:`expected_order` within
        tolerance and :attr:`r_squared` is high (:math:`R^2 \ge 0.95`).
    spacing_type : {'spatial', 'temporal'}
        Discretization domain refined during the study.
    details : dict[str, Any]
        Additional context such as solver name, test case, and tolerance thresholds.

    Examples
    --------
    >>> import numpy as np
    >>> from softpotato.analytical.benchmark import estimate_convergence_order
    >>> h = np.array([0.1, 0.05, 0.025])
    >>> err = 0.5 * h**2  # pure 2nd-order error
    >>> p, r2 = estimate_convergence_order(h, err)
    >>> round(p, 2)
    2.0
    >>> round(r2, 4)
    1.0
    """

    spacings: np.ndarray
    errors: np.ndarray
    observed_order: float
    expected_order: float
    r_squared: float
    passed: bool
    spacing_type: Literal["spatial", "temporal"]
    details: dict[str, Any] = field(default_factory=dict)

BenchmarkCaseResult dataclass

Execution and validation outcome for a single benchmark problem.

Parameters:

Name	Type	Description	Default
case_name	str	Name of the benchmark problem (e.g. 'fourier_decay', 'cottrell').	required
solver_name	str	Name of the numerical solver evaluated.	required
concentration_metrics	BenchmarkMetrics	Standardized error metrics for the spatial/temporal concentration field.	required
flux_metrics	BenchmarkMetrics or None	Standardized error metrics for interfacial boundary flux or current, or None if not applicable for the case.	None
execution_time	float	Wall-clock simulation execution time in seconds.	0.0
passed	bool	Whether all error metrics satisfied the predefined tolerance criteria.	True
details	dict[str, Any]	Diagnostic metadata and parameters used for the benchmark run.	dict()

Source code in src/softpotato/analytical/benchmark.py

@dataclass(frozen=True, kw_only=True)
class BenchmarkCaseResult:
    """Execution and validation outcome for a single benchmark problem.

    Parameters
    ----------
    case_name : str
        Name of the benchmark problem (e.g. ``'fourier_decay'``, ``'cottrell'``).
    solver_name : str
        Name of the numerical solver evaluated.
    concentration_metrics : BenchmarkMetrics
        Standardized error metrics for the spatial/temporal concentration field.
    flux_metrics : BenchmarkMetrics or None
        Standardized error metrics for interfacial boundary flux or current,
        or None if not applicable for the case.
    execution_time : float
        Wall-clock simulation execution time in seconds.
    passed : bool
        Whether all error metrics satisfied the predefined tolerance criteria.
    details : dict[str, Any]
        Diagnostic metadata and parameters used for the benchmark run.
    """

    case_name: str
    solver_name: str
    concentration_metrics: BenchmarkMetrics
    flux_metrics: BenchmarkMetrics | None = None
    execution_time: float = 0.0
    passed: bool = True
    details: dict[str, Any] = field(default_factory=dict)

BenchmarkSuiteResult dataclass

Aggregated outcome of running a benchmark suite across solvers and test cases.

Parameters:

Name	Type	Description	Default
cases	list[BenchmarkCaseResult]	Collection of all individual benchmark run results.	required
total_time	float	Total wall-clock execution time for the entire suite in seconds.	required
all_passed	bool	Flag indicating if all individual benchmark cases passed their tolerances.	required

Source code in src/softpotato/analytical/benchmark.py

@dataclass(frozen=True, kw_only=True)
class BenchmarkSuiteResult:
    """Aggregated outcome of running a benchmark suite across solvers and test cases.

    Parameters
    ----------
    cases : list[BenchmarkCaseResult]
        Collection of all individual benchmark run results.
    total_time : float
        Total wall-clock execution time for the entire suite in seconds.
    all_passed : bool
        Flag indicating if all individual benchmark cases passed their tolerances.
    """

    cases: list[BenchmarkCaseResult]
    total_time: float
    all_passed: bool

    def summary(self) -> str:
        """Generate a formatted ASCII comparison table of all benchmark runs."""
        headers = [
            "Benchmark Case",
            "Solver",
            "Status",
            "Conc. RMSE",
            "Conc. L_inf",
            "Flux RMSE",
            "Time (s)",
        ]
        rows: list[list[str]] = []
        for res in self.cases:
            status = "PASSED" if res.passed else "FAILED"
            c_rmse = f"{res.concentration_metrics.rmse:.3e}"
            c_linf = f"{res.concentration_metrics.max_abs_error:.3e}"
            f_rmse = (
                f"{res.flux_metrics.rmse:.3e}"
                if res.flux_metrics is not None
                else "N/A"
            )
            time_str = f"{res.execution_time:.4f}"
            rows.append(
                [
                    res.case_name,
                    res.solver_name,
                    status,
                    c_rmse,
                    c_linf,
                    f_rmse,
                    time_str,
                ]
            )

        col_widths = [len(h) for h in headers]
        for row in rows:
            for i, val in enumerate(row):
                col_widths[i] = max(col_widths[i], len(val))

        sep = "+-" + "-+-".join("-" * w for w in col_widths) + "-+"
        header_row = (
            "| "
            + " | ".join(h.ljust(w) for h, w in zip(headers, col_widths, strict=True))
            + " |"
        )
        data_rows = [
            "| "
            + " | ".join(val.ljust(w) for val, w in zip(row, col_widths, strict=True))
            + " |"
            for row in rows
        ]

        title = f"Soft Potato Solver Validation Suite (Total Time: {self.total_time:.3f}s | All Passed: {self.all_passed})"
        return "\n".join([title, sep, header_row, sep, *data_rows, sep])

summary

summary() -> str

Generate a formatted ASCII comparison table of all benchmark runs.

Source code in src/softpotato/analytical/benchmark.py

def summary(self) -> str:
    """Generate a formatted ASCII comparison table of all benchmark runs."""
    headers = [
        "Benchmark Case",
        "Solver",
        "Status",
        "Conc. RMSE",
        "Conc. L_inf",
        "Flux RMSE",
        "Time (s)",
    ]
    rows: list[list[str]] = []
    for res in self.cases:
        status = "PASSED" if res.passed else "FAILED"
        c_rmse = f"{res.concentration_metrics.rmse:.3e}"
        c_linf = f"{res.concentration_metrics.max_abs_error:.3e}"
        f_rmse = (
            f"{res.flux_metrics.rmse:.3e}"
            if res.flux_metrics is not None
            else "N/A"
        )
        time_str = f"{res.execution_time:.4f}"
        rows.append(
            [
                res.case_name,
                res.solver_name,
                status,
                c_rmse,
                c_linf,
                f_rmse,
                time_str,
            ]
        )

    col_widths = [len(h) for h in headers]
    for row in rows:
        for i, val in enumerate(row):
            col_widths[i] = max(col_widths[i], len(val))

    sep = "+-" + "-+-".join("-" * w for w in col_widths) + "-+"
    header_row = (
        "| "
        + " | ".join(h.ljust(w) for h, w in zip(headers, col_widths, strict=True))
        + " |"
    )
    data_rows = [
        "| "
        + " | ".join(val.ljust(w) for val, w in zip(row, col_widths, strict=True))
        + " |"
        for row in rows
    ]

    title = f"Soft Potato Solver Validation Suite (Total Time: {self.total_time:.3f}s | All Passed: {self.all_passed})"
    return "\n".join([title, sep, header_row, sep, *data_rows, sep])

GridComparisonEntry dataclass

Performance and accuracy metrics for a single spatial mesh configuration.

Parameters:

Name	Type	Description	Default
label	str	Human-readable label describing the mesh configuration (e.g. 'geometric (N=60)').	required
mesh_kind	str	Mesh family identifier ('uniform', 'geometric', or 'exponential').	required
n_points	int	Total number of spatial grid nodes \(N_x\).	required
dx_min	float	Minimum interval spacing \(\Delta x_{\min}\) across the mesh.	required
dx_max	float	Maximum interval spacing \(\Delta x_{\max}\) across the mesh.	required
execution_time	float	Wall-clock solver runtime in seconds.	required
memory_bytes	int	Memory footprint in bytes of the stored concentration state array(s).	required
current_metrics	BenchmarkMetrics	Standardized error metrics (\(L_2\), \(L_\infty\), RMSE, relative error) for the Faradaic current transient or voltammogram.	required
concentration_metrics	BenchmarkMetrics or None	Standardized error metrics for the spatial concentration field, if available.	None
details	dict[str, Any]	Additional diagnostic metadata (e.g., peak current error in cyclic voltammetry).	dict()

Source code in src/softpotato/analytical/benchmark.py

@dataclass(frozen=True, kw_only=True)
class GridComparisonEntry:
    r"""Performance and accuracy metrics for a single spatial mesh configuration.

    Parameters
    ----------
    label : str
        Human-readable label describing the mesh configuration (e.g. ``'geometric (N=60)'``).
    mesh_kind : str
        Mesh family identifier (``'uniform'``, ``'geometric'``, or ``'exponential'``).
    n_points : int
        Total number of spatial grid nodes :math:`N_x`.
    dx_min : float
        Minimum interval spacing :math:`\Delta x_{\min}` across the mesh.
    dx_max : float
        Maximum interval spacing :math:`\Delta x_{\max}` across the mesh.
    execution_time : float
        Wall-clock solver runtime in seconds.
    memory_bytes : int
        Memory footprint in bytes of the stored concentration state array(s).
    current_metrics : BenchmarkMetrics
        Standardized error metrics (:math:`L_2`, :math:`L_\infty`, RMSE, relative error)
        for the Faradaic current transient or voltammogram.
    concentration_metrics : BenchmarkMetrics or None
        Standardized error metrics for the spatial concentration field, if available.
    details : dict[str, Any]
        Additional diagnostic metadata (e.g., peak current error in cyclic voltammetry).
    """

    label: str
    mesh_kind: str
    n_points: int
    dx_min: float
    dx_max: float
    execution_time: float
    memory_bytes: int
    current_metrics: BenchmarkMetrics
    concentration_metrics: BenchmarkMetrics | None = None
    details: dict[str, Any] = field(default_factory=dict)

    @property
    def l2_error(self) -> float:
        r"""Discrete relative :math:`L_2` error norm of the Faradaic current."""
        return self.current_metrics.l2_error

    @property
    def max_rel_error(self) -> float:
        r"""Maximum relative error of the Faradaic current (or peak current error for CV)."""
        if "peak_current_rel_error" in self.details:
            return float(self.details["peak_current_rel_error"])
        return self.current_metrics.max_rel_error

l2_error property

l2_error: float

Discrete relative \(L_2\) error norm of the Faradaic current.

max_rel_error property

max_rel_error: float

Maximum relative error of the Faradaic current (or peak current error for CV).

GridComparisonResult dataclass

Comparative benchmark outcome evaluating uniform vs. expanding meshes.

Parameters:

Name	Type	Description	Default
case_name	str	Electrochemical benchmark problem evaluated (e.g. 'cottrell', 'spherical_cottrell', 'cyclic_voltammetry').	required
solver_name	str	Name of the numerical solver evaluated.	required
entries	list[GridComparisonEntry]	Ordered list of evaluated mesh configurations and their error/timing metrics.	required
point_reduction_factor	float	Ratio of uniform grid points to expanding grid points required for target precision.	required
speedup_factor	float	Wall-clock runtime speedup of the expanding mesh relative to the finest uniform mesh.	required
memory_reduction_factor	float	Memory footprint reduction ratio of the expanding mesh relative to the finest uniform mesh.	required
passed	bool	Whether the expanding mesh met the target precision and efficiency criteria.	required

Source code in src/softpotato/analytical/benchmark.py

@dataclass(frozen=True, kw_only=True)
class GridComparisonResult:
    r"""Comparative benchmark outcome evaluating uniform vs. expanding meshes.

    Parameters
    ----------
    case_name : str
        Electrochemical benchmark problem evaluated (e.g. ``'cottrell'``,
        ``'spherical_cottrell'``, ``'cyclic_voltammetry'``).
    solver_name : str
        Name of the numerical solver evaluated.
    entries : list[GridComparisonEntry]
        Ordered list of evaluated mesh configurations and their error/timing metrics.
    point_reduction_factor : float
        Ratio of uniform grid points to expanding grid points required for target precision.
    speedup_factor : float
        Wall-clock runtime speedup of the expanding mesh relative to the finest uniform mesh.
    memory_reduction_factor : float
        Memory footprint reduction ratio of the expanding mesh relative to the finest uniform mesh.
    passed : bool
        Whether the expanding mesh met the target precision and efficiency criteria.
    """

    case_name: str
    solver_name: str
    entries: list[GridComparisonEntry]
    point_reduction_factor: float
    speedup_factor: float
    memory_reduction_factor: float
    passed: bool

    def summary(self) -> str:
        """Generate a formatted ASCII table comparing runtime, memory, and accuracy across meshes."""
        headers = [
            "Mesh Configuration",
            "Kind",
            "N_points",
            "dx_min (m)",
            "Current L_2",
            "Max Rel Err (%)",
            "Memory (KB)",
            "Time (ms)",
        ]
        rows: list[list[str]] = []
        for entry in self.entries:
            rows.append(
                [
                    entry.label,
                    entry.mesh_kind,
                    str(entry.n_points),
                    f"{entry.dx_min:.3e}",
                    f"{entry.current_metrics.l2_error:.3e}",
                    f"{entry.current_metrics.max_rel_error * 100.0:.4f}",
                    f"{entry.memory_bytes / 1024.0:.2f}",
                    f"{entry.execution_time * 1000.0:.2f}",
                ]
            )

        col_widths = [len(h) for h in headers]
        for row in rows:
            for i, val in enumerate(row):
                col_widths[i] = max(col_widths[i], len(val))

        sep = "+-" + "-+-".join("-" * w for w in col_widths) + "-+"
        header_row = (
            "| "
            + " | ".join(h.ljust(w) for h, w in zip(headers, col_widths, strict=True))
            + " |"
        )
        data_rows = [
            "| "
            + " | ".join(val.ljust(w) for val, w in zip(row, col_widths, strict=True))
            + " |"
            for row in rows
        ]

        title = (
            f"Grid Efficiency Benchmark: {self.case_name} (Solver: {self.solver_name} | "
            f"Point Reduction: {self.point_reduction_factor:.1f}x | "
            f"Speedup: {self.speedup_factor:.2f}x | Passed: {self.passed})"
        )
        return "\n".join([title, sep, header_row, sep, *data_rows, sep])

summary

summary() -> str

Generate a formatted ASCII table comparing runtime, memory, and accuracy across meshes.

Source code in src/softpotato/analytical/benchmark.py

def summary(self) -> str:
    """Generate a formatted ASCII table comparing runtime, memory, and accuracy across meshes."""
    headers = [
        "Mesh Configuration",
        "Kind",
        "N_points",
        "dx_min (m)",
        "Current L_2",
        "Max Rel Err (%)",
        "Memory (KB)",
        "Time (ms)",
    ]
    rows: list[list[str]] = []
    for entry in self.entries:
        rows.append(
            [
                entry.label,
                entry.mesh_kind,
                str(entry.n_points),
                f"{entry.dx_min:.3e}",
                f"{entry.current_metrics.l2_error:.3e}",
                f"{entry.current_metrics.max_rel_error * 100.0:.4f}",
                f"{entry.memory_bytes / 1024.0:.2f}",
                f"{entry.execution_time * 1000.0:.2f}",
            ]
        )

    col_widths = [len(h) for h in headers]
    for row in rows:
        for i, val in enumerate(row):
            col_widths[i] = max(col_widths[i], len(val))

    sep = "+-" + "-+-".join("-" * w for w in col_widths) + "-+"
    header_row = (
        "| "
        + " | ".join(h.ljust(w) for h, w in zip(headers, col_widths, strict=True))
        + " |"
    )
    data_rows = [
        "| "
        + " | ".join(val.ljust(w) for val, w in zip(row, col_widths, strict=True))
        + " |"
        for row in rows
    ]

    title = (
        f"Grid Efficiency Benchmark: {self.case_name} (Solver: {self.solver_name} | "
        f"Point Reduction: {self.point_reduction_factor:.1f}x | "
        f"Speedup: {self.speedup_factor:.2f}x | Passed: {self.passed})"
    )
    return "\n".join([title, sep, header_row, sep, *data_rows, sep])

compute_error_metrics

compute_error_metrics(numerical: ArrayLike, exact: ArrayLike, *, eps: float = 1e-12) -> BenchmarkMetrics

Compute standardized error norms between numerical results and exact ground truths.

Parameters:

Name	Type	Description	Default
numerical	array - like	Simulated numerical values (scalar, 1D, or multidimensional array).	required
exact	array - like	Exact analytical reference values corresponding to numerical.	required
eps	float	Small positive constant added to denominators to prevent division by zero when evaluating relative errors against zero-valued ground truths.	1e-12

Returns:

Type	Description
BenchmarkMetrics	Dataclass containing RMSE, \(L_\infty\) (max_abs_error), discrete \(L_2\) relative error, MAE, mean relative error, and maximum relative error.

Raises:

Type	Description
ValueError	If input arrays cannot be broadcast together or have different shapes.

Examples:

>>> import numpy as np
>>> from softpotato.analytical.benchmark import compute_error_metrics
>>> y_num = np.array([1.0, 2.05, 3.0])
>>> y_exact = np.array([1.0, 2.0, 3.0])
>>> m = compute_error_metrics(y_num, y_exact)
>>> round(m.max_abs_error, 4)
0.05
>>> round(m.rmse, 4)
0.0289

Source code in src/softpotato/analytical/benchmark.py

def compute_error_metrics(
    numerical: ArrayLike,
    exact: ArrayLike,
    *,
    eps: float = 1e-12,
) -> BenchmarkMetrics:
    r"""Compute standardized error norms between numerical results and exact ground truths.

    Parameters
    ----------
    numerical : array-like
        Simulated numerical values (scalar, 1D, or multidimensional array).
    exact : array-like
        Exact analytical reference values corresponding to ``numerical``.
    eps : float, default 1e-12
        Small positive constant added to denominators to prevent division by zero
        when evaluating relative errors against zero-valued ground truths.

    Returns
    -------
    BenchmarkMetrics
        Dataclass containing RMSE, :math:`L_\infty` (:attr:`max_abs_error`),
        discrete :math:`L_2` relative error, MAE, mean relative error, and
        maximum relative error.

    Raises
    ------
    ValueError
        If input arrays cannot be broadcast together or have different shapes.

    Examples
    --------
    >>> import numpy as np
    >>> from softpotato.analytical.benchmark import compute_error_metrics
    >>> y_num = np.array([1.0, 2.05, 3.0])
    >>> y_exact = np.array([1.0, 2.0, 3.0])
    >>> m = compute_error_metrics(y_num, y_exact)
    >>> round(m.max_abs_error, 4)
    0.05
    >>> round(m.rmse, 4)
    0.0289
    """
    num_arr = np.asarray(numerical, dtype=float)
    exact_arr = np.asarray(exact, dtype=float)

    if num_arr.shape != exact_arr.shape:
        try:
            num_arr, exact_arr = np.broadcast_arrays(num_arr, exact_arr)
        except ValueError as err:
            raise ValueError(
                f"Shape mismatch: numerical shape {num_arr.shape} cannot be broadcast "
                f"with exact reference shape {exact_arr.shape}."
            ) from err

    diff = num_arr - exact_arr
    abs_diff = np.abs(diff)

    max_abs = float(np.max(abs_diff))
    mean_abs = float(np.mean(abs_diff))
    rmse = float(np.sqrt(np.mean(diff**2)))

    norm_exact = float(np.linalg.norm(exact_arr))
    if norm_exact > eps:
        l2_err = float(np.linalg.norm(diff) / norm_exact)
    else:
        l2_err = float(np.linalg.norm(diff) / math.sqrt(exact_arr.size))

    rel_diff = abs_diff / (np.abs(exact_arr) + eps)
    mean_rel = float(np.mean(rel_diff))
    max_rel = float(np.max(rel_diff))

    return BenchmarkMetrics(
        rmse=rmse,
        max_abs_error=max_abs,
        l2_error=l2_err,
        mean_abs_error=mean_abs,
        mean_rel_error=mean_rel,
        max_rel_error=max_rel,
    )

estimate_convergence_order

estimate_convergence_order(spacings: ArrayLike, errors: ArrayLike) -> tuple[float, float]

Estimate empirical order of convergence \(p\) and \(R^2\) via log-log regression.

Fits the power law \(\text{error} = C \cdot h^p\) in logarithmic space:

\[ \ln(\text{error}) = p \cdot \ln(h) + \ln(C) \]

Parameters:

Name	Type	Description	Default
spacings	array - like	Array of strictly positive discretization increments \(h\) (\(\Delta x\) or \(\Delta t\)).	required
errors	array - like	Array of strictly positive corresponding error norms.	required

Returns:

Name	Type	Description
order	float	Empirical convergence order \(p\) (slope of log-log regression).
r_squared	float	Coefficient of determination \(R^2 \in [0, 1]\) indicating the linearity and reliability of the power-law convergence regime.

Raises:

Type	Description
ValueError	If inputs have fewer than 2 data points or contain non-positive values.

Examples:

>>> import numpy as np
>>> from softpotato.analytical.benchmark import estimate_convergence_order
>>> h = np.array([0.1, 0.05, 0.025, 0.0125])
>>> err = 2.5 * h**2
>>> p, r2 = estimate_convergence_order(h, err)
>>> round(p, 4)
2.0
>>> round(r2, 4)
1.0

Source code in src/softpotato/analytical/benchmark.py

def estimate_convergence_order(
    spacings: ArrayLike,
    errors: ArrayLike,
) -> tuple[float, float]:
    r"""Estimate empirical order of convergence :math:`p` and :math:`R^2` via log-log regression.

    Fits the power law :math:`\text{error} = C \cdot h^p` in logarithmic space:

    .. math::

        \ln(\text{error}) = p \cdot \ln(h) + \ln(C)

    Parameters
    ----------
    spacings : array-like
        Array of strictly positive discretization increments :math:`h`
        (:math:`\Delta x` or :math:`\Delta t`).
    errors : array-like
        Array of strictly positive corresponding error norms.

    Returns
    -------
    order : float
        Empirical convergence order :math:`p` (slope of log-log regression).
    r_squared : float
        Coefficient of determination :math:`R^2 \in [0, 1]` indicating the linearity
        and reliability of the power-law convergence regime.

    Raises
    ------
    ValueError
        If inputs have fewer than 2 data points or contain non-positive values.

    Examples
    --------
    >>> import numpy as np
    >>> from softpotato.analytical.benchmark import estimate_convergence_order
    >>> h = np.array([0.1, 0.05, 0.025, 0.0125])
    >>> err = 2.5 * h**2
    >>> p, r2 = estimate_convergence_order(h, err)
    >>> round(p, 4)
    2.0
    >>> round(r2, 4)
    1.0
    """
    h_arr = np.asarray(spacings, dtype=float)
    e_arr = np.asarray(errors, dtype=float)

    if h_arr.ndim != 1 or e_arr.ndim != 1:
        raise ValueError("spacings and errors must be 1D arrays.")
    if len(h_arr) != len(e_arr):
        raise ValueError(
            f"Length mismatch: len(spacings)={len(h_arr)} != len(errors)={len(e_arr)}."
        )
    if len(h_arr) < 2:
        raise ValueError(
            f"At least 2 points are required to compute convergence order, got {len(h_arr)}."
        )
    if np.any(h_arr <= 0):
        raise ValueError("All spacings must be strictly positive (h > 0).")
    if np.any(e_arr <= 0):
        raise ValueError("All errors must be strictly positive (error > 0).")

    log_h = np.log(h_arr)
    log_e = np.log(e_arr)

    # Ordinary least-squares regression: log_e = p * log_h + intercept
    cov = np.cov(log_h, log_e)
    var_h = cov[0, 0]
    cov_he = cov[0, 1]

    if var_h <= 1e-15:
        raise ValueError(
            "Spacings do not have sufficient variation to determine slope."
        )

    order = float(cov_he / var_h)
    intercept = float(np.mean(log_e) - order * np.mean(log_h))

    # Calculate R^2
    fitted = order * log_h + intercept
    ss_tot = float(np.sum((log_e - np.mean(log_e)) ** 2))
    ss_res = float(np.sum((log_e - fitted) ** 2))

    r_squared = 1.0 - (ss_res / ss_tot) if ss_tot > 1e-15 else 1.0
    r_squared = float(np.clip(r_squared, 0.0, 1.0))

    return order, r_squared

run_fourier_decay_benchmark

run_fourier_decay_benchmark(solver: BaseSolver | str = 'crank_nicolson', *, n_points: int = 51, grid: Mesh1D | ndarray | None = None, D: float = 0.05, L: float = 1.0, t_end: float = 0.2, n_times: int = 11, max_error_threshold: float = 0.005, **solver_kwargs: Any) -> BenchmarkCaseResult

Run 1D Fourier sine diffusion decay benchmark comparing against exact solution.

Governing PDE:

\[ \frac{\partial c}{\partial t} = D \frac{\partial^2 c}{\partial x^2}, \quad x \in [0, L] \]

Initial and boundary conditions:

\[ c(x, 0) = \sin\left(\frac{\pi x}{L}\right), \quad c(0, t) = c(L, t) = 0 \]

Exact analytical solution:

\[ c(x, t) = \sin\left(\frac{\pi x}{L}\right) \exp\left(-D \left(\frac{\pi}{L}\right)^2 t\right) \]

Exact surface flux at \(x = 0\):

\[ J(0, t) = -D \left.\frac{\partial c}{\partial x}\right|_{x=0} = -D \frac{\pi}{L} \exp\left(-D \left(\frac{\pi}{L}\right)^2 t\right) \]

Parameters:

Name	Type	Description	Default
solver	BaseSolver or str	Numerical solver instance or registered string name.	"crank_nicolson"
n_points	int	Number of spatial grid points \(N_x\) (used when grid is None).	51
grid	Mesh1D or ndarray	Custom 1D spatial mesh or coordinate array. If None, uses a uniform grid spanning \([0, L]\) with n_points nodes.	None
D	float	Diffusion coefficient \(D > 0\).	0.05
L	float	Domain length \(L > 0\).	1.0
t_end	float	Simulation end time \(t_{\text{end}} > 0\).	0.2
n_times	int	Number of time evaluation points across \([0, t_{\text{end}}]\).	11
max_error_threshold	float	Maximum allowable \(L_\infty\) concentration error to consider the run passed.	5e-3
**solver_kwargs	Any	Additional keyword arguments forwarded to the solver (e.g. dt).	{}

Returns:

Type	Description
BenchmarkCaseResult	Outcome containing concentration metrics, interfacial flux metrics, wall-clock runtime, and pass/fail indicator.

Examples:

>>> from softpotato.analytical.benchmark import run_fourier_decay_benchmark
>>> res = run_fourier_decay_benchmark("crank_nicolson", n_points=41, t_end=0.1)
>>> res.passed
True
>>> res.concentration_metrics.max_abs_error < 5e-3
True

Source code in src/softpotato/analytical/benchmark.py

def run_fourier_decay_benchmark(
    solver: BaseSolver | str = "crank_nicolson",
    *,
    n_points: int = 51,
    grid: Mesh1D | np.ndarray | None = None,
    D: float = 0.05,
    L: float = 1.0,
    t_end: float = 0.2,
    n_times: int = 11,
    max_error_threshold: float = 5e-3,
    **solver_kwargs: Any,
) -> BenchmarkCaseResult:
    r"""Run 1D Fourier sine diffusion decay benchmark comparing against exact solution.

    Governing PDE:

    .. math::

        \frac{\partial c}{\partial t} = D \frac{\partial^2 c}{\partial x^2}, \quad x \in [0, L]

    Initial and boundary conditions:

    .. math::

        c(x, 0) = \sin\left(\frac{\pi x}{L}\right), \quad c(0, t) = c(L, t) = 0

    Exact analytical solution:

    .. math::

        c(x, t) = \sin\left(\frac{\pi x}{L}\right) \exp\left(-D \left(\frac{\pi}{L}\right)^2 t\right)

    Exact surface flux at :math:`x = 0`:

    .. math::

        J(0, t) = -D \left.\frac{\partial c}{\partial x}\right|_{x=0} = -D \frac{\pi}{L} \exp\left(-D \left(\frac{\pi}{L}\right)^2 t\right)

    Parameters
    ----------
    solver : BaseSolver or str, default "crank_nicolson"
        Numerical solver instance or registered string name.
    n_points : int, default 51
        Number of spatial grid points :math:`N_x` (used when ``grid`` is None).
    grid : Mesh1D or np.ndarray, optional
        Custom 1D spatial mesh or coordinate array. If None, uses a uniform grid
        spanning :math:`[0, L]` with ``n_points`` nodes.
    D : float, default 0.05
        Diffusion coefficient :math:`D > 0`.
    L : float, default 1.0
        Domain length :math:`L > 0`.
    t_end : float, default 0.2
        Simulation end time :math:`t_{\text{end}} > 0`.
    n_times : int, default 11
        Number of time evaluation points across :math:`[0, t_{\text{end}}]`.
    max_error_threshold : float, default 5e-3
        Maximum allowable :math:`L_\infty` concentration error to consider the run passed.
    **solver_kwargs : Any
        Additional keyword arguments forwarded to the solver (e.g. ``dt``).

    Returns
    -------
    BenchmarkCaseResult
        Outcome containing concentration metrics, interfacial flux metrics,
        wall-clock runtime, and pass/fail indicator.

    Examples
    --------
    >>> from softpotato.analytical.benchmark import run_fourier_decay_benchmark
    >>> res = run_fourier_decay_benchmark("crank_nicolson", n_points=41, t_end=0.1)
    >>> res.passed
    True
    >>> res.concentration_metrics.max_abs_error < 5e-3
    True
    """
    solver_inst, solver_name = _resolve_solver_instance(solver)

    mesh_or_grid: Mesh1D | np.ndarray = (
        np.linspace(0.0, L, n_points) if grid is None else grid
    )
    grid_nodes = (
        mesh_or_grid.nodes
        if isinstance(mesh_or_grid, Mesh1D)
        else np.asarray(mesh_or_grid, dtype=float)
    )
    domain_len = float(grid_nodes[-1] - grid_nodes[0])
    x0 = float(grid_nodes[0])
    t_eval = np.linspace(0.0, t_end, n_times)

    problem = DiffusionProblem(
        grid=mesh_or_grid,
        diffusivity={"c": D},
        boundary_conditions={"c": (DirichletBC(0.0), DirichletBC(0.0))},
        initial_conditions={
            "c": lambda x_pts: np.sin(np.pi * (x_pts - x0) / domain_len)
        },
    )

    t0 = time.perf_counter()
    result = solver_inst.solve(
        problem, t_span=(0.0, t_end), t_eval=t_eval, **solver_kwargs
    )
    wall_time = time.perf_counter() - t0

    if not result.success:
        failed_metrics = BenchmarkMetrics(
            rmse=float("inf"),
            max_abs_error=float("inf"),
            l2_error=float("inf"),
            mean_abs_error=float("inf"),
            mean_rel_error=float("inf"),
            max_rel_error=float("inf"),
        )
        return BenchmarkCaseResult(
            case_name="fourier_decay",
            solver_name=solver_name,
            concentration_metrics=failed_metrics,
            flux_metrics=None,
            execution_time=wall_time,
            passed=False,
            details={"error_message": result.message},
        )

    # Exact 2D grid evaluation
    t_grid, x_grid = np.meshgrid(t_eval, grid_nodes, indexing="ij")
    decay_rate = D * (np.pi / domain_len) ** 2
    c_exact = np.sin(np.pi * (x_grid - x0) / domain_len) * np.exp(-decay_rate * t_grid)
    c_num = result["c"]

    conc_metrics = compute_error_metrics(c_num, c_exact)

    # Exact interfacial flux J(0, t)
    j_exact = -D * (np.pi / domain_len) * np.exp(-decay_rate * t_eval)
    j_num = result.fluxes.get("c", None)
    flux_metrics = compute_error_metrics(j_num, j_exact) if j_num is not None else None

    passed = bool(result.success and conc_metrics.max_abs_error <= max_error_threshold)

    return BenchmarkCaseResult(
        case_name="fourier_decay",
        solver_name=solver_name,
        concentration_metrics=conc_metrics,
        flux_metrics=flux_metrics,
        execution_time=wall_time,
        passed=passed,
        details={
            "n_points": len(grid_nodes),
            "D": D,
            "L": domain_len,
            "t_end": t_end,
            "dx": float(grid_nodes[1] - grid_nodes[0]),
            "dx_min": problem.dx_min,
            "dx_max": problem.dx_max,
            "is_uniform": problem.is_uniform,
            "max_error_threshold": max_error_threshold,
        },
    )

run_cottrell_benchmark

run_cottrell_benchmark(solver: BaseSolver | str = 'crank_nicolson', *, n_points: int = 101, grid: Mesh1D | ndarray | None = None, mesh_kind: Literal['uniform', 'geometric', 'exponential'] = 'uniform', gamma: float = 1.07, beta: float = 3.0, D: float = 1e-05, c_bulk: float = 0.001, n_electrons: int = 1, area: float = 0.0001, t_span: tuple[float, float] = (0.01, 0.2), n_times: int = 20, domain_factor: float = 6.0, max_error_threshold: float = 0.05, **solver_kwargs: Any) -> BenchmarkCaseResult

Run Cottrell potential step benchmark comparing concentration profiles and flux.

Simulates chronoamperometric semi-infinite linear diffusion following a step to zero surface concentration (\(c(0, t) = 0\)). The domain length is automatically dimensioned as \(L = \text{domain\_factor} \times \sqrt{D \cdot t_{\text{end}}}\) to ensure the right boundary satisfies semi-infinite boundary conditions (\(c(L, t) \approx c^*\)).

Exact analytical spatial concentration profile:

\[ c(x, t) = c^* \text{erf}\left(\frac{x}{2 \sqrt{D t}}\right) \]

Exact Cottrell current transient:

\[ I(t) = \frac{n F A \sqrt{D} c^*}{\sqrt{\pi t}} \]

Parameters:

Name	Type	Description	Default
solver	BaseSolver or str	Numerical solver instance or registered string name.	"crank_nicolson"
n_points	int	Number of spatial grid points \(N_x\).	101
grid	Mesh1D or ndarray	Custom 1D spatial mesh or coordinate array. Overrides mesh_kind and n_points.	None
mesh_kind	('uniform', 'geometric', 'exponential')	Mesh generator used when grid is None.	'uniform'
gamma	float	Expansion ratio \(\gamma > 1\) when mesh_kind='geometric'.	1.07
beta	float	Stretching parameter \(\beta > 0\) when mesh_kind='exponential'.	3.0
D	float	Diffusion coefficient \(D > 0\).	1e-5
c_bulk	float	Bulk concentration \(c^* > 0\).	1e-3
n_electrons	int	Number of electrons transferred per redox event \(n \ge 1\).	1
area	float	Electrode area in \(\text{m}^2\) (or \(\text{cm}^2\)).	1e-4
t_span	tuple[float, float]	Evaluation time interval \((t_{\min}, t_{\max})\) with \(t_{\min} > 0\) to evaluate beyond the initial Cottrell step discontinuity.	(0.01, 0.2)
n_times	int	Number of evaluation time points across t_span.	20
domain_factor	float	Diffusion layer multiplier setting domain length \(L = \text{domain\_factor} \sqrt{D t_{\max}}\).	6.0
max_error_threshold	float	Maximum allowable \(L_\infty\) concentration error to consider the run passed.	5e-2
**solver_kwargs	Any	Additional keyword arguments forwarded to the solver.	{}

Returns:

Type	Description
BenchmarkCaseResult	Outcome containing concentration metrics, Cottrell current metrics, runtime, and pass/fail indicator.

Examples:

>>> from softpotato.analytical.benchmark import run_cottrell_benchmark
>>> res = run_cottrell_benchmark("scipy_ivp", n_points=81, t_span=(0.02, 0.1))
>>> res.passed
True
>>> res.concentration_metrics.max_abs_error < 0.05
True

Source code in src/softpotato/analytical/benchmark.py

def run_cottrell_benchmark(
    solver: BaseSolver | str = "crank_nicolson",
    *,
    n_points: int = 101,
    grid: Mesh1D | np.ndarray | None = None,
    mesh_kind: Literal["uniform", "geometric", "exponential"] = "uniform",
    gamma: float = 1.07,
    beta: float = 3.0,
    D: float = 1e-5,
    c_bulk: float = 1e-3,
    n_electrons: int = 1,
    area: float = 1e-4,
    t_span: tuple[float, float] = (0.01, 0.2),
    n_times: int = 20,
    domain_factor: float = 6.0,
    max_error_threshold: float = 5e-2,
    **solver_kwargs: Any,
) -> BenchmarkCaseResult:
    r"""Run Cottrell potential step benchmark comparing concentration profiles and flux.

    Simulates chronoamperometric semi-infinite linear diffusion following a step to
    zero surface concentration (:math:`c(0, t) = 0`). The domain length is automatically
    dimensioned as :math:`L = \text{domain\_factor} \times \sqrt{D \cdot t_{\text{end}}}`
    to ensure the right boundary satisfies semi-infinite boundary conditions
    (:math:`c(L, t) \approx c^*`).

    Exact analytical spatial concentration profile:

    .. math::

        c(x, t) = c^* \text{erf}\left(\frac{x}{2 \sqrt{D t}}\right)

    Exact Cottrell current transient:

    .. math::

        I(t) = \frac{n F A \sqrt{D} c^*}{\sqrt{\pi t}}

    Parameters
    ----------
    solver : BaseSolver or str, default "crank_nicolson"
        Numerical solver instance or registered string name.
    n_points : int, default 101
        Number of spatial grid points :math:`N_x`.
    grid : Mesh1D or np.ndarray, optional
        Custom 1D spatial mesh or coordinate array. Overrides ``mesh_kind`` and ``n_points``.
    mesh_kind : {'uniform', 'geometric', 'exponential'}, default 'uniform'
        Mesh generator used when ``grid`` is None.
    gamma : float, default 1.07
        Expansion ratio :math:`\gamma > 1` when ``mesh_kind='geometric'``.
    beta : float, default 3.0
        Stretching parameter :math:`\beta > 0` when ``mesh_kind='exponential'``.
    D : float, default 1e-5
        Diffusion coefficient :math:`D > 0`.
    c_bulk : float, default 1e-3
        Bulk concentration :math:`c^* > 0`.
    n_electrons : int, default 1
        Number of electrons transferred per redox event :math:`n \ge 1`.
    area : float, default 1e-4
        Electrode area in :math:`\text{m}^2` (or :math:`\text{cm}^2`).
    t_span : tuple[float, float], default (0.01, 0.2)
        Evaluation time interval :math:`(t_{\min}, t_{\max})` with :math:`t_{\min} > 0`
        to evaluate beyond the initial Cottrell step discontinuity.
    n_times : int, default 20
        Number of evaluation time points across :attr:`t_span`.
    domain_factor : float, default 6.0
        Diffusion layer multiplier setting domain length :math:`L = \text{domain\_factor} \sqrt{D t_{\max}}`.
    max_error_threshold : float, default 5e-2
        Maximum allowable :math:`L_\infty` concentration error to consider the run passed.
    **solver_kwargs : Any
        Additional keyword arguments forwarded to the solver.

    Returns
    -------
    BenchmarkCaseResult
        Outcome containing concentration metrics, Cottrell current metrics,
        runtime, and pass/fail indicator.

    Examples
    --------
    >>> from softpotato.analytical.benchmark import run_cottrell_benchmark
    >>> res = run_cottrell_benchmark("scipy_ivp", n_points=81, t_span=(0.02, 0.1))
    >>> res.passed
    True
    >>> res.concentration_metrics.max_abs_error < 0.05
    True
    """
    solver_inst, solver_name = _resolve_solver_instance(solver)

    t_start, t_end = t_span
    if t_start <= 0.0:
        raise ValueError(
            f"t_span start must be strictly positive (t_start > 0), got {t_start}."
        )
    if t_start >= t_end:
        raise ValueError(
            f"t_span end must be strictly greater than start, got {t_span}."
        )

    if grid is not None:
        mesh_or_grid: Mesh1D | np.ndarray = grid
    elif mesh_kind == "uniform":
        l_domain = domain_factor * math.sqrt(D * t_end)
        mesh_or_grid = np.linspace(0.0, l_domain, n_points)
    else:
        mesh_or_grid = create_diffusion_mesh(
            D=D,
            t_max=t_end,
            kind=mesh_kind,
            factor=domain_factor,
            n_points=n_points,
            gamma=gamma,
            beta=beta,
        )

    grid_nodes = (
        mesh_or_grid.nodes
        if isinstance(mesh_or_grid, Mesh1D)
        else np.asarray(mesh_or_grid, dtype=float)
    )
    l_domain = float(grid_nodes[-1] - grid_nodes[0])

    # Solve from 0.0 to t_end so transient starts at t=0
    # Include 0.0 in integration grid, then evaluate metrics for t >= t_start
    dt_eval = (t_end - t_start) / max(1, n_times - 1)
    n_pre_steps = max(1, math.ceil(t_start / dt_eval))
    pre_t = np.linspace(0.0, t_start, n_pre_steps + 1)[:-1]
    t_eval_window = np.linspace(t_start, t_end, n_times)
    full_t_eval = np.concatenate([pre_t, t_eval_window])

    problem = DiffusionProblem(
        grid=mesh_or_grid,
        diffusivity={"c": D},
        boundary_conditions={"c": (DirichletBC(0.0), DirichletBC(c_bulk))},
        initial_conditions={"c": c_bulk},
    )

    t0 = time.perf_counter()
    result = solver_inst.solve(
        problem, t_span=(0.0, t_end), t_eval=full_t_eval, **solver_kwargs
    )
    wall_time = time.perf_counter() - t0

    if not result.success:
        failed_metrics = BenchmarkMetrics(
            rmse=float("inf"),
            max_abs_error=float("inf"),
            l2_error=float("inf"),
            mean_abs_error=float("inf"),
            mean_rel_error=float("inf"),
            max_rel_error=float("inf"),
        )
        return BenchmarkCaseResult(
            case_name="cottrell_step",
            solver_name=solver_name,
            concentration_metrics=failed_metrics,
            flux_metrics=None,
            execution_time=wall_time,
            passed=False,
            details={"error_message": result.message},
        )

    # Filter to evaluation time window [t_start, t_end]
    mask = result.t >= (t_start - 1e-12)
    t_eval = result.t[mask]
    c_num = result["c"][mask, :]

    # Exact concentration field across all evaluated times and grid points
    c_exact = np.zeros_like(c_num)
    for k, t_k in enumerate(t_eval):
        c_exact[k, :] = step_concentration_profile(
            grid_nodes, t=float(t_k), D=D, c_bulk=c_bulk
        )

    conc_metrics = compute_error_metrics(c_num, c_exact)

    # Compare Cottrell currents: for an oxidation, I_num = -n * F * A * J(0) > 0
    # J(0) = -D * dc/dx < 0 when reactant is consumed at x=0, so I_num > 0 matches IUPAC anodic sign.
    flux_metrics = None
    if "c" in result.fluxes:
        j_num = result.fluxes["c"][mask]
        i_num = -n_electrons * FARADAY * area * j_num
        i_exact = cottrell(t_eval, n=n_electrons, D=D, c_bulk=c_bulk, area=area)
        flux_metrics = compute_error_metrics(i_num, i_exact)

    passed = bool(result.success and conc_metrics.max_abs_error <= max_error_threshold)

    return BenchmarkCaseResult(
        case_name="cottrell_step",
        solver_name=solver_name,
        concentration_metrics=conc_metrics,
        flux_metrics=flux_metrics,
        execution_time=wall_time,
        passed=passed,
        details={
            "n_points": len(grid_nodes),
            "D": D,
            "c_bulk": c_bulk,
            "t_span": t_span,
            "domain_length": l_domain,
            "dx_min": problem.dx_min,
            "dx_max": problem.dx_max,
            "is_uniform": problem.is_uniform,
            "max_error_threshold": max_error_threshold,
        },
    )

run_spherical_cottrell_benchmark

run_spherical_cottrell_benchmark(solver: BaseSolver | str = 'scipy_ivp', *, n_points: int = 60, grid: Mesh1D | ndarray | None = None, mesh_kind: Literal['geometric', 'exponential', 'uniform'] = 'geometric', gamma: float = 1.07, beta: float = 3.0, r0: float = 2.5e-05, D: float = 1e-09, c_bulk: float = 1.0, n_electrons: int = 1, hemisphere: bool = False, t_span: tuple[float, float] = (0.02, 1.0), n_times: int = 20, domain_factor: float = 6.0, max_rel_error_threshold: float = 0.01, **solver_kwargs: Any) -> BenchmarkCaseResult

Run spherical Cottrell potential step benchmark against exact analytical theory.

Simulates radial diffusion to a spherical or hemispherical electrode of radius \(r_0\) following a potential step to zero surface concentration (\(c(r_0, t) = 0\)):

\[ \frac{\partial c}{\partial t} = D \left( \frac{\partial^2 c}{\partial r^2} + \frac{2}{r} \frac{\partial c}{\partial r} \right), \quad r \in [r_0, r_{\max}] \]

Exact analytical concentration profile:

\[ c(r, t) = c^* \left[ 1 - \frac{r_0}{r} \operatorname{erfc}\left( \frac{r - r_0}{2 \sqrt{D t}} \right) \right] \]

Exact analytical spherical Cottrell current:

\[ I(t) = n F A D c^* \left( \frac{1}{\sqrt{\pi D t}} + \frac{1}{r_0} \right) \]

Parameters:

Name	Type	Description	Default
solver	BaseSolver or str	Numerical solver instance or registered string name.	"scipy_ivp"
n_points	int	Number of spatial grid points \(N_r\).	60
grid	Mesh1D or ndarray	Custom 1D radial mesh or coordinate array starting at \(r_0 > 0\).	None
mesh_kind	('geometric', 'exponential', 'uniform')	Mesh generator used when grid is None.	'geometric'
gamma	float	Expansion ratio \(\gamma > 1\) when mesh_kind='geometric'.	1.07
beta	float	Stretching parameter \(\beta > 0\) when mesh_kind='exponential'.	3.0
r0	float	Electrode radius \(r_0 > 0\) in meters.	25e-6
D	float	Diffusion coefficient \(D > 0\) in \(\text{m}^2/\text{s}\).	1e-9
c_bulk	float	Bulk concentration \(c^* > 0\) in \(\text{mol}/\text{m}^3\).	1.0
n_electrons	int	Number of electrons transferred per redox event \(n \ge 1\).	1
hemisphere	bool	If True, uses hemispherical surface area \(A = 2\pi r_0^2\) instead of full sphere \(A = 4\pi r_0^2\).	False
t_span	tuple[float, float]	Evaluation time window \((t_{\min}, t_{\max})\) with \(t_{\min} > 0\).	(0.02, 1.0)
n_times	int	Number of evaluation time points across t_span.	20
domain_factor	float	Diffusion layer multiplier setting outer boundary \(r_{\max} = r_0 + \text{domain\_factor} \sqrt{D t_{\max}}\).	6.0
max_rel_error_threshold	float	Maximum allowable relative current error to consider the run passed.	1e-2
**solver_kwargs	Any	Additional keyword arguments forwarded to the solver.	{}

Returns:

Type	Description
BenchmarkCaseResult	Outcome containing concentration metrics, spherical Cottrell current metrics, wall-clock runtime, and pass/fail indicator.

Source code in src/softpotato/analytical/benchmark.py

def run_spherical_cottrell_benchmark(
    solver: BaseSolver | str = "scipy_ivp",
    *,
    n_points: int = 60,
    grid: Mesh1D | np.ndarray | None = None,
    mesh_kind: Literal["geometric", "exponential", "uniform"] = "geometric",
    gamma: float = 1.07,
    beta: float = 3.0,
    r0: float = 25e-6,
    D: float = 1e-9,
    c_bulk: float = 1.0,
    n_electrons: int = 1,
    hemisphere: bool = False,
    t_span: tuple[float, float] = (0.02, 1.0),
    n_times: int = 20,
    domain_factor: float = 6.0,
    max_rel_error_threshold: float = 1e-2,
    **solver_kwargs: Any,
) -> BenchmarkCaseResult:
    r"""Run spherical Cottrell potential step benchmark against exact analytical theory.

    Simulates radial diffusion to a spherical or hemispherical electrode of radius
    :math:`r_0` following a potential step to zero surface concentration
    (:math:`c(r_0, t) = 0`):

    .. math::

        \frac{\partial c}{\partial t} = D \left( \frac{\partial^2 c}{\partial r^2} + \frac{2}{r} \frac{\partial c}{\partial r} \right), \quad r \in [r_0, r_{\max}]

    Exact analytical concentration profile:

    .. math::

        c(r, t) = c^* \left[ 1 - \frac{r_0}{r} \operatorname{erfc}\left( \frac{r - r_0}{2 \sqrt{D t}} \right) \right]

    Exact analytical spherical Cottrell current:

    .. math::

        I(t) = n F A D c^* \left( \frac{1}{\sqrt{\pi D t}} + \frac{1}{r_0} \right)

    Parameters
    ----------
    solver : BaseSolver or str, default "scipy_ivp"
        Numerical solver instance or registered string name.
    n_points : int, default 60
        Number of spatial grid points :math:`N_r`.
    grid : Mesh1D or np.ndarray, optional
        Custom 1D radial mesh or coordinate array starting at :math:`r_0 > 0`.
    mesh_kind : {'geometric', 'exponential', 'uniform'}, default 'geometric'
        Mesh generator used when ``grid`` is None.
    gamma : float, default 1.07
        Expansion ratio :math:`\gamma > 1` when ``mesh_kind='geometric'``.
    beta : float, default 3.0
        Stretching parameter :math:`\beta > 0` when ``mesh_kind='exponential'``.
    r0 : float, default 25e-6
        Electrode radius :math:`r_0 > 0` in meters.
    D : float, default 1e-9
        Diffusion coefficient :math:`D > 0` in :math:`\text{m}^2/\text{s}`.
    c_bulk : float, default 1.0
        Bulk concentration :math:`c^* > 0` in :math:`\text{mol}/\text{m}^3`.
    n_electrons : int, default 1
        Number of electrons transferred per redox event :math:`n \ge 1`.
    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`.
    t_span : tuple[float, float], default (0.02, 1.0)
        Evaluation time window :math:`(t_{\min}, t_{\max})` with :math:`t_{\min} > 0`.
    n_times : int, default 20
        Number of evaluation time points across :attr:`t_span`.
    domain_factor : float, default 6.0
        Diffusion layer multiplier setting outer boundary :math:`r_{\max} = r_0 + \text{domain\_factor} \sqrt{D t_{\max}}`.
    max_rel_error_threshold : float, default 1e-2
        Maximum allowable relative current error to consider the run passed.
    **solver_kwargs : Any
        Additional keyword arguments forwarded to the solver.

    Returns
    -------
    BenchmarkCaseResult
        Outcome containing concentration metrics, spherical Cottrell current metrics,
        wall-clock runtime, and pass/fail indicator.
    """
    solver_inst, solver_name = _resolve_solver_instance(solver)

    t_start, t_end = t_span
    if t_start <= 0.0:
        raise ValueError(
            f"t_span start must be strictly positive (t_start > 0), got {t_start}."
        )
    if t_start >= t_end:
        raise ValueError(
            f"t_span end must be strictly greater than start, got {t_span}."
        )

    if grid is not None:
        mesh_or_grid: Mesh1D | np.ndarray = grid
        grid_nodes = (
            mesh_or_grid.nodes
            if isinstance(mesh_or_grid, Mesh1D)
            else np.asarray(mesh_or_grid, dtype=float)
        )
        r0_eff = float(grid_nodes[0])
    else:
        r0_eff = float(r0)
        mesh_or_grid = create_diffusion_mesh(
            D=D,
            t_max=t_end,
            kind=mesh_kind,
            x_min=r0_eff,
            factor=domain_factor,
            n_points=n_points,
            gamma=gamma,
            beta=beta,
        )
        grid_nodes = mesh_or_grid.nodes

    geom = SphericalGeometry(r0=r0_eff, hemisphere=hemisphere)
    t_eval = np.linspace(t_start, t_end, n_times)

    problem = DiffusionProblem(
        grid=mesh_or_grid,
        diffusivity={"c": D},
        boundary_conditions={"c": (DirichletBC(0.0), DirichletBC(c_bulk))},
        initial_conditions={"c": c_bulk},
        geometry=geom,
        n_electrons=n_electrons,
    )

    t0 = time.perf_counter()
    result = solver_inst.solve(
        problem, t_span=(0.0, t_end), t_eval=t_eval, **solver_kwargs
    )
    wall_time = time.perf_counter() - t0

    if not result.success:
        failed_metrics = BenchmarkMetrics(
            rmse=float("inf"),
            max_abs_error=float("inf"),
            l2_error=float("inf"),
            mean_abs_error=float("inf"),
            mean_rel_error=float("inf"),
            max_rel_error=float("inf"),
        )
        return BenchmarkCaseResult(
            case_name="spherical_cottrell",
            solver_name=solver_name,
            concentration_metrics=failed_metrics,
            flux_metrics=None,
            execution_time=wall_time,
            passed=False,
            details={"error_message": result.message},
        )

    c_num = result["c"]
    c_exact = np.empty_like(c_num)
    for k, t_k in enumerate(t_eval):
        eta = (grid_nodes - r0_eff) / (2.0 * math.sqrt(D * float(t_k)))
        c_exact[k, :] = c_bulk * (1.0 - (r0_eff / grid_nodes) * erfc(eta))

    conc_metrics = compute_error_metrics(c_num, c_exact)

    i_num = np.abs(result.currents["c"])
    i_exact = cottrell_spherical(
        t_eval, r=r0_eff, n=n_electrons, D=D, c_bulk=c_bulk, area=geom.area
    )
    flux_metrics = compute_error_metrics(i_num, i_exact)

    passed = bool(
        result.success and flux_metrics.max_rel_error <= max_rel_error_threshold
    )

    return BenchmarkCaseResult(
        case_name="spherical_cottrell",
        solver_name=solver_name,
        concentration_metrics=conc_metrics,
        flux_metrics=flux_metrics,
        execution_time=wall_time,
        passed=passed,
        details={
            "n_points": len(grid_nodes),
            "r0": r0_eff,
            "D": D,
            "c_bulk": c_bulk,
            "t_span": t_span,
            "dx_min": problem.dx_min,
            "dx_max": problem.dx_max,
            "is_uniform": problem.is_uniform,
            "max_rel_error_threshold": max_rel_error_threshold,
        },
    )

run_mass_conservation_benchmark

run_mass_conservation_benchmark(solver: BaseSolver | str = 'crank_nicolson', *, n_points: int = 61, grid: Mesh1D | ndarray | None = None, mesh_kind: Literal['uniform', 'geometric', 'exponential'] = 'uniform', gamma: float = 1.05, beta: float = 2.0, two_sided: bool = True, D: float = 0.02, L: float = 1.0, t_span: tuple[float, float] = (0.0, 0.2), n_times: int = 11, max_drift_threshold: float = 0.001, **solver_kwargs: Any) -> BenchmarkCaseResult

Run mass conservation benchmark under zero-flux Neumann boundary conditions.

Solves 1D diffusion on a closed interval with zero-flux boundaries:

\[ \left.\frac{\partial c}{\partial x}\right|_{x=0} = 0, \quad \left.\frac{\partial c}{\partial x}\right|_{x=L} = 0 \]

initialized with an internal Gaussian concentration pulse:

\[ c(x, 0) = \exp\left(-50 \left(x - \frac{L}{2}\right)^2\right) \]

Analytical conservation invariant:

\[ M(t) = \int_0^L c(x, t) \, dx = M(0) = \text{constant}, \quad \forall t \ge 0 \]

Parameters:

Name	Type	Description	Default
solver	BaseSolver or str	Numerical solver instance or registered string name.	"crank_nicolson"
n_points	int	Number of spatial grid points \(N_x\).	61
grid	Mesh1D or ndarray	Custom 1D spatial mesh or coordinate array. Overrides mesh_kind and n_points.	None
mesh_kind	('uniform', 'geometric', 'exponential')	Mesh generator used when grid is None.	'uniform'
gamma	float	Expansion ratio \(\gamma > 1\) when mesh_kind='geometric'.	1.05
beta	float	Stretching parameter \(\beta > 0\) when mesh_kind='exponential'.	2.0
two_sided	bool	Whether to use symmetric two-sided node clustering for non-uniform meshes.	True
D	float	Diffusion coefficient \(D > 0\).	0.02
L	float	Domain length \(L > 0\).	1.0
t_span	tuple[float, float]	Simulation time span \((t_0, t_{\text{end}})\).	(0.0, 0.2)
n_times	int	Number of evaluation time points.	11
max_drift_threshold	float	Maximum relative mass drift \(|M(t) - M(0)| / M(0)\) allowed.	1e-3
**solver_kwargs	Any	Additional keyword arguments forwarded to the solver.	{}

Returns:

Type	Description
BenchmarkCaseResult	Outcome containing mass error metrics, execution time, and pass/fail indicator.

Examples:

>>> from softpotato.analytical.benchmark import run_mass_conservation_benchmark
>>> res = run_mass_conservation_benchmark("crank_nicolson", n_points=51)
>>> res.passed
True
>>> res.details["max_relative_drift"] < 1e-3
True

Source code in src/softpotato/analytical/benchmark.py

def run_mass_conservation_benchmark(
    solver: BaseSolver | str = "crank_nicolson",
    *,
    n_points: int = 61,
    grid: Mesh1D | np.ndarray | None = None,
    mesh_kind: Literal["uniform", "geometric", "exponential"] = "uniform",
    gamma: float = 1.05,
    beta: float = 2.0,
    two_sided: bool = True,
    D: float = 0.02,
    L: float = 1.0,
    t_span: tuple[float, float] = (0.0, 0.2),
    n_times: int = 11,
    max_drift_threshold: float = 1e-3,
    **solver_kwargs: Any,
) -> BenchmarkCaseResult:
    r"""Run mass conservation benchmark under zero-flux Neumann boundary conditions.

    Solves 1D diffusion on a closed interval with zero-flux boundaries:

    .. math::

        \left.\frac{\partial c}{\partial x}\right|_{x=0} = 0, \quad \left.\frac{\partial c}{\partial x}\right|_{x=L} = 0

    initialized with an internal Gaussian concentration pulse:

    .. math::

        c(x, 0) = \exp\left(-50 \left(x - \frac{L}{2}\right)^2\right)

    Analytical conservation invariant:

    .. math::

        M(t) = \int_0^L c(x, t) \, dx = M(0) = \text{constant}, \quad \forall t \ge 0

    Parameters
    ----------
    solver : BaseSolver or str, default "crank_nicolson"
        Numerical solver instance or registered string name.
    n_points : int, default 61
        Number of spatial grid points :math:`N_x`.
    grid : Mesh1D or np.ndarray, optional
        Custom 1D spatial mesh or coordinate array. Overrides ``mesh_kind`` and ``n_points``.
    mesh_kind : {'uniform', 'geometric', 'exponential'}, default 'uniform'
        Mesh generator used when ``grid`` is None.
    gamma : float, default 1.05
        Expansion ratio :math:`\gamma > 1` when ``mesh_kind='geometric'``.
    beta : float, default 2.0
        Stretching parameter :math:`\beta > 0` when ``mesh_kind='exponential'``.
    two_sided : bool, default True
        Whether to use symmetric two-sided node clustering for non-uniform meshes.
    D : float, default 0.02
        Diffusion coefficient :math:`D > 0`.
    L : float, default 1.0
        Domain length :math:`L > 0`.
    t_span : tuple[float, float], default (0.0, 0.2)
        Simulation time span :math:`(t_0, t_{\text{end}})`.
    n_times : int, default 11
        Number of evaluation time points.
    max_drift_threshold : float, default 1e-3
        Maximum relative mass drift :math:`|M(t) - M(0)| / M(0)` allowed.
    **solver_kwargs : Any
        Additional keyword arguments forwarded to the solver.

    Returns
    -------
    BenchmarkCaseResult
        Outcome containing mass error metrics, execution time, and pass/fail indicator.

    Examples
    --------
    >>> from softpotato.analytical.benchmark import run_mass_conservation_benchmark
    >>> res = run_mass_conservation_benchmark("crank_nicolson", n_points=51)
    >>> res.passed
    True
    >>> res.details["max_relative_drift"] < 1e-3
    True
    """
    solver_inst, solver_name = _resolve_solver_instance(solver)

    if grid is not None:
        mesh_or_grid: Mesh1D | np.ndarray = grid
    elif mesh_kind == "uniform":
        mesh_or_grid = np.linspace(0.0, L, n_points)
    elif mesh_kind == "geometric":
        mesh_or_grid = GeometricExpandingMesh(
            0.0, L, n_points=n_points, gamma=gamma, two_sided=two_sided
        )
    elif mesh_kind == "exponential":
        mesh_or_grid = ExponentialMesh(
            0.0, L, n_points=n_points, beta=beta, two_sided=two_sided
        )
    else:
        raise ValueError(
            f"Unknown mesh_kind '{mesh_kind}'. Expected 'uniform', 'geometric', or 'exponential'."
        )

    grid_nodes = (
        mesh_or_grid.nodes
        if isinstance(mesh_or_grid, Mesh1D)
        else np.asarray(mesh_or_grid, dtype=float)
    )
    t_eval = np.linspace(t_span[0], t_span[1], n_times)

    pulse_center = 0.5 * float(grid_nodes[0] + grid_nodes[-1])
    init_func = lambda x_pts: np.exp(-50.0 * (x_pts - pulse_center) ** 2)

    problem = DiffusionProblem(
        grid=mesh_or_grid,
        diffusivity={"c": D},
        boundary_conditions={"c": (NeumannBC(0.0), NeumannBC(0.0))},
        initial_conditions={"c": init_func},
    )

    t0 = time.perf_counter()
    result = solver_inst.solve(problem, t_span=t_span, t_eval=t_eval, **solver_kwargs)
    wall_time = time.perf_counter() - t0

    if not result.success:
        failed_metrics = BenchmarkMetrics(
            rmse=float("inf"),
            max_abs_error=float("inf"),
            l2_error=float("inf"),
            mean_abs_error=float("inf"),
            mean_rel_error=float("inf"),
            max_rel_error=float("inf"),
        )
        return BenchmarkCaseResult(
            case_name="mass_conservation",
            solver_name=solver_name,
            concentration_metrics=failed_metrics,
            flux_metrics=None,
            execution_time=wall_time,
            passed=False,
            details={"error_message": result.message},
        )

    # Compute numerical total mass M(t) = integral c(x, t) dx
    c = result["c"]
    masses = np.array([trapezoid(c[k, :], grid_nodes) for k in range(len(t_eval))])
    m0 = masses[0]

    exact_masses = np.full_like(masses, m0)
    mass_metrics = compute_error_metrics(masses, exact_masses)

    rel_drifts = np.abs(masses - m0) / (abs(m0) + 1e-15)
    max_drift = float(np.max(rel_drifts))

    passed = bool(result.success and max_drift <= max_drift_threshold)

    return BenchmarkCaseResult(
        case_name="mass_conservation",
        solver_name=solver_name,
        concentration_metrics=mass_metrics,
        flux_metrics=None,
        execution_time=wall_time,
        passed=passed,
        details={
            "n_points": len(grid_nodes),
            "initial_mass": float(m0),
            "final_mass": float(masses[-1]),
            "max_relative_drift": max_drift,
            "max_drift_threshold": max_drift_threshold,
            "is_uniform": problem.is_uniform,
        },
    )

verify_spatial_convergence

verify_spatial_convergence(solver: BaseSolver | str = 'crank_nicolson', *, n_points_list: Sequence[int] = (21, 41, 81, 161), mesh_kind: Literal['uniform', 'geometric', 'exponential'] = 'uniform', expansion_span: float = 4.0, beta: float = 1.5, two_sided: bool = False, mapping: str = 'exponential', grid_factory: Callable[[int], Mesh1D | ndarray] | None = None, target: Literal['concentration', 'flux'] = 'concentration', D: float = 0.05, L: float = 1.0, t_fixed: float = 0.1, dt_fine: float = 0.0001, expected_order: float = 2.0, order_tolerance: float = 0.25, metric: Literal['max_abs_error', 'rmse', 'l2_error'] = 'max_abs_error', **solver_kwargs: Any) -> ConvergenceResult

Verify spatial discretization order across grid refinement \(\Delta x\).

Executes a sequence of spatial grid refinements on the Fourier sine decay problem, using a fine time step dt_fine so spatial error dominates. Computes the empirical convergence order \(p\) via log-log regression:

\[ \ln(\text{error}) = p \cdot \ln(\Delta x) + C \]

Supports uniform meshes (Uniform1DMesh), geometrically expanding meshes (GeometricExpandingMesh with fixed total stretching span \(\Delta x_{\max} / \Delta x_{\min} = \text{expansion\_span}\) so that \(\gamma_N = \text{expansion\_span}^{1/(N-1)}\) and adjacent step differences scale as \(\Delta x_{j+1} - \Delta x_j = \mathcal{O}(\Delta x^2)\)), and smooth coordinate-mapped exponential meshes (ExponentialMesh).

Parameters:

Name	Type	Description	Default
solver	BaseSolver or str	Numerical solver instance or registered string name.	"crank_nicolson"
n_points_list	sequence of int	List of spatial grid point counts to evaluate.	(21, 41, 81, 161)
mesh_kind	('uniform', 'geometric', 'exponential')	Spatial mesh family to refine when grid_factory is None.	'uniform'
expansion_span	float	Global ratio \(\Delta x_{\max} / \Delta x_{\min} > 1\) maintained across refinements when mesh_kind='geometric'.	4.0
beta	float	Stretching parameter \(\beta > 0\) when mesh_kind='exponential'.	1.5
two_sided	bool	Whether to use symmetric two-sided clustering for expanding meshes.	False
mapping	str	Coordinate transformation ('exponential', 'tanh', or 'sinh') when mesh_kind='exponential'.	'exponential'
grid_factory	Callable[[int], Mesh1D | ndarray]	Custom callable f(n_points) returning a 1D mesh or array for each refinement level.	None
target	('concentration', 'flux')	Quantity whose discretization error is regressed against characteristic spacing: spatial concentration profile or interfacial surface flux \(J(0, t)\).	'concentration'
D	float	Diffusion coefficient \(D > 0\).	0.05
L	float	Domain length \(L > 0\).	1.0
t_fixed	float	Fixed time at which error is evaluated.	0.1
dt_fine	float	Small time step size ensuring temporal error is negligible.	1e-4
expected_order	float	Theoretical spatial order of convergence (default 2.0 for central differences).	2.0
order_tolerance	float	Permissible deficit below expected_order to consider the test passed.	0.25
metric	('max_abs_error', 'rmse', 'l2_error')	Error norm to regress against \(\Delta x\).	'max_abs_error'
**solver_kwargs	Any	Additional keyword arguments forwarded to the solver.	{}

Returns:

Type	Description
ConvergenceResult	Dataclass containing grid spacings \(\Delta x\), error array, empirical order \(p\), \(R^2\), and pass/fail indicator.

Examples:

>>> from softpotato.analytical.benchmark import verify_spatial_convergence
>>> res = verify_spatial_convergence("crank_nicolson", n_points_list=(21, 41, 81))
>>> res.passed
True
>>> round(res.observed_order, 1)
2.0

Source code in src/softpotato/analytical/benchmark.py

def verify_spatial_convergence(
    solver: BaseSolver | str = "crank_nicolson",
    *,
    n_points_list: Sequence[int] = (21, 41, 81, 161),
    mesh_kind: Literal["uniform", "geometric", "exponential"] = "uniform",
    expansion_span: float = 4.0,
    beta: float = 1.5,
    two_sided: bool = False,
    mapping: str = "exponential",
    grid_factory: Callable[[int], Mesh1D | np.ndarray] | None = None,
    target: Literal["concentration", "flux"] = "concentration",
    D: float = 0.05,
    L: float = 1.0,
    t_fixed: float = 0.1,
    dt_fine: float = 1e-4,
    expected_order: float = 2.0,
    order_tolerance: float = 0.25,
    metric: Literal["max_abs_error", "rmse", "l2_error"] = "max_abs_error",
    **solver_kwargs: Any,
) -> ConvergenceResult:
    r"""Verify spatial discretization order across grid refinement :math:`\Delta x`.

    Executes a sequence of spatial grid refinements on the Fourier sine decay problem,
    using a fine time step :attr:`dt_fine` so spatial error dominates. Computes the
    empirical convergence order :math:`p` via log-log regression:

    .. math::

        \ln(\text{error}) = p \cdot \ln(\Delta x) + C

    Supports uniform meshes (:class:`~softpotato.grid.Uniform1DMesh`), geometrically
    expanding meshes (:class:`~softpotato.grid.GeometricExpandingMesh` with fixed
    total stretching span :math:`\Delta x_{\max} / \Delta x_{\min} = \text{expansion\_span}`
    so that :math:`\gamma_N = \text{expansion\_span}^{1/(N-1)}` and adjacent step differences
    scale as :math:`\Delta x_{j+1} - \Delta x_j = \mathcal{O}(\Delta x^2)`), and smooth
    coordinate-mapped exponential meshes (:class:`~softpotato.grid.ExponentialMesh`).

    Parameters
    ----------
    solver : BaseSolver or str, default "crank_nicolson"
        Numerical solver instance or registered string name.
    n_points_list : sequence of int, default (21, 41, 81, 161)
        List of spatial grid point counts to evaluate.
    mesh_kind : {'uniform', 'geometric', 'exponential'}, default 'uniform'
        Spatial mesh family to refine when ``grid_factory`` is None.
    expansion_span : float, default 4.0
        Global ratio :math:`\Delta x_{\max} / \Delta x_{\min} > 1` maintained across
        refinements when ``mesh_kind='geometric'``.
    beta : float, default 1.5
        Stretching parameter :math:`\beta > 0` when ``mesh_kind='exponential'``.
    two_sided : bool, default False
        Whether to use symmetric two-sided clustering for expanding meshes.
    mapping : str, default 'exponential'
        Coordinate transformation (``'exponential'``, ``'tanh'``, or ``'sinh'``)
        when ``mesh_kind='exponential'``.
    grid_factory : Callable[[int], Mesh1D | np.ndarray], optional
        Custom callable ``f(n_points)`` returning a 1D mesh or array for each refinement level.
    target : {'concentration', 'flux'}, default 'concentration'
        Quantity whose discretization error is regressed against characteristic spacing:
        spatial concentration profile or interfacial surface flux :math:`J(0, t)`.
    D : float, default 0.05
        Diffusion coefficient :math:`D > 0`.
    L : float, default 1.0
        Domain length :math:`L > 0`.
    t_fixed : float, default 0.1
        Fixed time at which error is evaluated.
    dt_fine : float, default 1e-4
        Small time step size ensuring temporal error is negligible.
    expected_order : float, default 2.0
        Theoretical spatial order of convergence (default 2.0 for central differences).
    order_tolerance : float, default 0.25
        Permissible deficit below :attr:`expected_order` to consider the test passed.
    metric : {'max_abs_error', 'rmse', 'l2_error'}, default 'max_abs_error'
        Error norm to regress against :math:`\Delta x`.
    **solver_kwargs : Any
        Additional keyword arguments forwarded to the solver.

    Returns
    -------
    ConvergenceResult
        Dataclass containing grid spacings :math:`\Delta x`, error array,
        empirical order :math:`p`, :math:`R^2`, and pass/fail indicator.

    Examples
    --------
    >>> from softpotato.analytical.benchmark import verify_spatial_convergence
    >>> res = verify_spatial_convergence("crank_nicolson", n_points_list=(21, 41, 81))
    >>> res.passed
    True
    >>> round(res.observed_order, 1)
    2.0
    """
    solver_inst, solver_name = _resolve_solver_instance(solver)
    if target not in ("concentration", "flux"):
        raise ValueError(
            f"Unknown target '{target}'. Expected 'concentration' or 'flux'."
        )
    if mesh_kind == "geometric" and expansion_span <= 1.0:
        raise ValueError(
            f"expansion_span must be > 1.0 for geometric meshes, got {expansion_span}."
        )

    dx_list: list[float] = []
    error_list: list[float] = []

    decay_factor = math.exp(-D * (math.pi / L) ** 2 * t_fixed)
    j_exact = -D * (math.pi / L) * decay_factor

    for n_pts in n_points_list:
        if grid_factory is not None:
            mesh_or_grid = grid_factory(n_pts)
        elif mesh_kind == "uniform":
            mesh_or_grid = np.linspace(0.0, L, n_pts)
        elif mesh_kind == "geometric":
            m_intervals = max(1, (n_pts - 1) // 2 if two_sided else (n_pts - 1))
            gamma_n = float(expansion_span) ** (1.0 / m_intervals)
            mesh_or_grid = GeometricExpandingMesh(
                0.0, L, n_points=n_pts, gamma=gamma_n, two_sided=two_sided
            )
        elif mesh_kind == "exponential":
            mesh_or_grid = ExponentialMesh(
                0.0,
                L,
                n_points=n_pts,
                beta=beta,
                two_sided=two_sided,
                mapping=mapping,
            )
        else:
            raise ValueError(
                f"Unknown mesh_kind '{mesh_kind}'. Expected 'uniform', 'geometric', or 'exponential'."
            )

        grid_nodes = (
            mesh_or_grid.nodes
            if isinstance(mesh_or_grid, Mesh1D)
            else np.asarray(mesh_or_grid, dtype=float)
        )
        dx = float((grid_nodes[-1] - grid_nodes[0]) / (len(grid_nodes) - 1))

        problem = DiffusionProblem(
            grid=mesh_or_grid,
            diffusivity={"c": D},
            boundary_conditions={"c": (DirichletBC(0.0), DirichletBC(0.0))},
            initial_conditions={"c": lambda x_pts: np.sin(np.pi * x_pts / L)},
        )

        kwargs = dict(solver_kwargs)
        if "dt" not in kwargs:
            kwargs["dt"] = dt_fine

        res = solver_inst.solve(
            problem,
            t_span=(0.0, t_fixed),
            t_eval=np.array([0.0, t_fixed]),
            **kwargs,
        )
        if not res.success:
            raise RuntimeError(
                f"Solver {solver_name} failed during spatial convergence run: {res.message}"
            )

        if target == "flux":
            j_num = float(res.fluxes["c"][-1])
            err_val = abs(j_num - j_exact)
        else:
            c_num = res["c"][-1, :]
            c_exact = np.sin(np.pi * grid_nodes / L) * decay_factor
            m = compute_error_metrics(c_num, c_exact)
            err_val = getattr(m, metric)

        dx_list.append(dx)
        error_list.append(err_val)

    dx_arr = np.array(dx_list, dtype=float)
    err_arr = np.array(error_list, dtype=float)

    order, r2 = estimate_convergence_order(dx_arr, err_arr)
    passed = bool(order >= (expected_order - order_tolerance) and r2 >= 0.95)

    return ConvergenceResult(
        spacings=dx_arr,
        errors=err_arr,
        observed_order=order,
        expected_order=expected_order,
        r_squared=r2,
        passed=passed,
        spacing_type="spatial",
        details={
            "solver_name": solver_name,
            "mesh_kind": mesh_kind,
            "target": target,
            "metric": metric,
            "n_points_list": list(n_points_list),
            "order_tolerance": order_tolerance,
        },
    )

verify_temporal_convergence

verify_temporal_convergence(solver: BaseSolver | str = 'crank_nicolson', *, dt_list: Sequence[float] = (0.04, 0.02, 0.01, 0.005), n_points: int = 51, D: float = 0.05, L: float = 1.0, t_end: float = 0.2, reference: Literal['fine_step', 'analytical'] = 'fine_step', dt_ref: float = 1e-05, expected_order: float = 2.0, order_tolerance: float = 0.3, metric: Literal['max_abs_error', 'rmse', 'l2_error'] = 'max_abs_error', **solver_kwargs: Any) -> ConvergenceResult

Verify temporal discretization order across time-step refinement \(\Delta t\).

Executes a sequence of temporal step refinements on the Fourier sine decay problem. To prevent spatial discretization errors from masking temporal convergence, the reference solution can be computed either analytically or via a high-resolution self-consistent numerical baseline (reference='fine_step').

Linear regression in logarithmic coordinates determines the temporal order \(p\):

\[ \ln(\text{error}) = p \cdot \ln(\Delta t) + C \]

Parameters:

Name	Type	Description	Default
solver	BaseSolver or str	Numerical solver instance or registered string name.	"crank_nicolson"
dt_list	sequence of float	Sequence of time-step sizes \(\Delta t\) to evaluate.	(0.04, 0.02, 0.01, 0.005)
n_points	int	Number of spatial grid points \(N_x\).	51
D	float	Diffusion coefficient \(D > 0\).	0.05
L	float	Domain length \(L > 0\).	1.0
t_end	float	Simulation end time \(t_{\text{end}}\).	0.2
reference	('fine_step', 'analytical')	Reference standard for error computation. 'fine_step' (recommended) runs the solver on the exact same spatial grid with dt_ref, cleanly canceling spatial discretization error and isolating the temporal order.	'fine_step'
dt_ref	float	Reference time-step used when reference='fine_step'.	1e-5
expected_order	float	Expected theoretical temporal order (e.g. 2.0 for Crank–Nicolson, 1.0 for BTCS).	2.0
order_tolerance	float	Permissible deficit below expected_order to consider the test passed.	0.3
metric	('max_abs_error', 'rmse', 'l2_error')	Error norm to regress against \(\Delta t\).	'max_abs_error'
**solver_kwargs	Any	Additional keyword arguments forwarded to the solver.	{}

Returns:

Type	Description
ConvergenceResult	Dataclass containing time steps \(\Delta t\), error array, empirical temporal order \(p\), \(R^2\), and pass/fail indicator.

Examples:

>>> from softpotato.analytical.benchmark import verify_temporal_convergence
>>> # Crank-Nicolson should exhibit 2nd-order temporal convergence (order ~ 2)
>>> res = verify_temporal_convergence("crank_nicolson", dt_list=(0.04, 0.02, 0.01), expected_order=2.0)
>>> res.passed
True
>>> round(res.observed_order, 1)
2.0

Source code in src/softpotato/analytical/benchmark.py

def verify_temporal_convergence(
    solver: BaseSolver | str = "crank_nicolson",
    *,
    dt_list: Sequence[float] = (0.04, 0.02, 0.01, 0.005),
    n_points: int = 51,
    D: float = 0.05,
    L: float = 1.0,
    t_end: float = 0.2,
    reference: Literal["fine_step", "analytical"] = "fine_step",
    dt_ref: float = 1e-5,
    expected_order: float = 2.0,
    order_tolerance: float = 0.3,
    metric: Literal["max_abs_error", "rmse", "l2_error"] = "max_abs_error",
    **solver_kwargs: Any,
) -> ConvergenceResult:
    r"""Verify temporal discretization order across time-step refinement :math:`\Delta t`.

    Executes a sequence of temporal step refinements on the Fourier sine decay problem.
    To prevent spatial discretization errors from masking temporal convergence, the
    reference solution can be computed either analytically or via a high-resolution
    self-consistent numerical baseline (:attr:`reference='fine_step'`).

    Linear regression in logarithmic coordinates determines the temporal order :math:`p`:

    .. math::

        \ln(\text{error}) = p \cdot \ln(\Delta t) + C

    Parameters
    ----------
    solver : BaseSolver or str, default "crank_nicolson"
        Numerical solver instance or registered string name.
    dt_list : sequence of float, default (0.04, 0.02, 0.01, 0.005)
        Sequence of time-step sizes :math:`\Delta t` to evaluate.
    n_points : int, default 51
        Number of spatial grid points :math:`N_x`.
    D : float, default 0.05
        Diffusion coefficient :math:`D > 0`.
    L : float, default 1.0
        Domain length :math:`L > 0`.
    t_end : float, default 0.2
        Simulation end time :math:`t_{\text{end}}`.
    reference : {'fine_step', 'analytical'}, default 'fine_step'
        Reference standard for error computation. ``'fine_step'`` (recommended)
        runs the solver on the exact same spatial grid with :attr:`dt_ref`,
        cleanly canceling spatial discretization error and isolating the temporal order.
    dt_ref : float, default 1e-5
        Reference time-step used when ``reference='fine_step'``.
    expected_order : float, default 2.0
        Expected theoretical temporal order (e.g. 2.0 for Crank–Nicolson, 1.0 for BTCS).
    order_tolerance : float, default 0.3
        Permissible deficit below :attr:`expected_order` to consider the test passed.
    metric : {'max_abs_error', 'rmse', 'l2_error'}, default 'max_abs_error'
        Error norm to regress against :math:`\Delta t`.
    **solver_kwargs : Any
        Additional keyword arguments forwarded to the solver.

    Returns
    -------
    ConvergenceResult
        Dataclass containing time steps :math:`\Delta t`, error array,
        empirical temporal order :math:`p`, :math:`R^2`, and pass/fail indicator.

    Examples
    --------
    >>> from softpotato.analytical.benchmark import verify_temporal_convergence
    >>> # Crank-Nicolson should exhibit 2nd-order temporal convergence (order ~ 2)
    >>> res = verify_temporal_convergence("crank_nicolson", dt_list=(0.04, 0.02, 0.01), expected_order=2.0)
    >>> res.passed
    True
    >>> round(res.observed_order, 1)
    2.0
    """
    solver_inst, solver_name = _resolve_solver_instance(solver)

    grid = np.linspace(0.0, L, n_points)
    problem = DiffusionProblem(
        grid=grid,
        diffusivity={"c": D},
        boundary_conditions={"c": (DirichletBC(0.0), DirichletBC(0.0))},
        initial_conditions={"c": lambda x_pts: np.sin(np.pi * x_pts / L)},
    )

    if reference == "fine_step":
        ref_kwargs = dict(solver_kwargs)
        ref_kwargs["dt"] = dt_ref
        ref_res = solver_inst.solve(
            problem,
            t_span=(0.0, t_end),
            t_eval=np.array([0.0, t_end]),
            **ref_kwargs,
        )
        if not ref_res.success:
            raise RuntimeError(
                f"Reference solver failed during temporal convergence run: {ref_res.message}"
            )
        c_ref = ref_res["c"][-1, :]
    else:
        decay_factor = math.exp(-D * (math.pi / L) ** 2 * t_end)
        c_ref = np.sin(np.pi * grid / L) * decay_factor

    dt_arr = np.array(dt_list, dtype=float)
    error_list: list[float] = []

    for dt_val in dt_arr:
        kwargs = dict(solver_kwargs)
        kwargs["dt"] = float(dt_val)

        res = solver_inst.solve(
            problem,
            t_span=(0.0, t_end),
            t_eval=np.array([0.0, t_end]),
            **kwargs,
        )
        if not res.success:
            raise RuntimeError(
                f"Solver {solver_name} failed at dt={dt_val}: {res.message}"
            )

        c_num = res["c"][-1, :]
        m = compute_error_metrics(c_num, c_ref)
        error_list.append(getattr(m, metric))

    err_arr = np.array(error_list, dtype=float)
    order, r2 = estimate_convergence_order(dt_arr, err_arr)
    passed = bool(order >= (expected_order - order_tolerance) and r2 >= 0.95)

    return ConvergenceResult(
        spacings=dt_arr,
        errors=err_arr,
        observed_order=order,
        expected_order=expected_order,
        r_squared=r2,
        passed=passed,
        spacing_type="temporal",
        details={
            "solver_name": solver_name,
            "metric": metric,
            "reference": reference,
            "dt_list": list(dt_list),
            "order_tolerance": order_tolerance,
        },
    )

list_benchmark_cases

list_benchmark_cases() -> list[str]

Return list of standard benchmark case names.

Source code in src/softpotato/analytical/benchmark.py

def list_benchmark_cases() -> list[str]:
    """Return list of standard benchmark case names."""
    return [
        "fourier_decay",
        "cottrell_step",
        "spherical_cottrell",
        "mass_conservation",
    ]

run_benchmark

run_benchmark(case: Literal['fourier_decay', 'cottrell_step', 'cottrell', 'spherical_cottrell', 'mass_conservation'], solver: BaseSolver | str = 'crank_nicolson', **kwargs: Any) -> BenchmarkCaseResult

Unified entry point to run a specific benchmark case with a solver.

Parameters:

Name	Type	Description	Default
case	('fourier_decay', 'cottrell_step', 'cottrell', 'spherical_cottrell', 'mass_conservation')	Benchmark problem to execute.	'fourier_decay'
solver	BaseSolver or str	Numerical solver instance or registered string name.	"crank_nicolson"
**kwargs	Any	Parameters passed to the specific benchmark runner function.	{}

Returns:

Type	Description
BenchmarkCaseResult	Outcome of the benchmark run.

Raises:

Type	Description
ValueError	If an unrecognized benchmark case name is provided.

Examples:

>>> from softpotato.analytical.benchmark import run_benchmark
>>> res = run_benchmark("fourier_decay", "crank_nicolson")
>>> res.passed
True

Source code in src/softpotato/analytical/benchmark.py

def run_benchmark(
    case: Literal[
        "fourier_decay",
        "cottrell_step",
        "cottrell",
        "spherical_cottrell",
        "mass_conservation",
    ],
    solver: BaseSolver | str = "crank_nicolson",
    **kwargs: Any,
) -> BenchmarkCaseResult:
    """Unified entry point to run a specific benchmark case with a solver.

    Parameters
    ----------
    case : {'fourier_decay', 'cottrell_step', 'cottrell', 'spherical_cottrell', 'mass_conservation'}
        Benchmark problem to execute.
    solver : BaseSolver or str, default "crank_nicolson"
        Numerical solver instance or registered string name.
    **kwargs : Any
        Parameters passed to the specific benchmark runner function.

    Returns
    -------
    BenchmarkCaseResult
        Outcome of the benchmark run.

    Raises
    ------
    ValueError
        If an unrecognized benchmark case name is provided.

    Examples
    --------
    >>> from softpotato.analytical.benchmark import run_benchmark
    >>> res = run_benchmark("fourier_decay", "crank_nicolson")
    >>> res.passed
    True
    """
    case_norm = case.lower().strip()
    if case_norm == "fourier_decay":
        return run_fourier_decay_benchmark(solver, **kwargs)
    if case_norm in ("cottrell", "cottrell_step"):
        return run_cottrell_benchmark(solver, **kwargs)
    if case_norm == "spherical_cottrell":
        return run_spherical_cottrell_benchmark(solver, **kwargs)
    if case_norm == "mass_conservation":
        return run_mass_conservation_benchmark(solver, **kwargs)

    raise ValueError(
        f"Unknown benchmark case '{case}'. Available cases: {list_benchmark_cases()}."
    )

run_benchmark_suite

run_benchmark_suite(solvers: Sequence[BaseSolver | str] | None = None, cases: Sequence[Literal['fourier_decay', 'cottrell_step', 'cottrell', 'spherical_cottrell', 'mass_conservation']] | None = None, **kwargs: Any) -> BenchmarkSuiteResult

Execute a comprehensive benchmark suite across multiple solvers and cases.

Parameters:

Name	Type	Description	Default
solvers	sequence of (BaseSolver or str)	List of solvers to evaluate. Defaults to unconditionally stable solvers: ('crank_nicolson', 'scipy_ivp', 'implicit').	None
cases	sequence of str	List of benchmark case names to run. Defaults to all standard cases: ('fourier_decay', 'cottrell_step', 'mass_conservation').	None
**kwargs	Any	Additional keyword arguments forwarded to benchmark runs.	{}

Returns:

Type	Description
BenchmarkSuiteResult	Aggregated outcomes of all benchmark runs with summary reporting.

Examples:

>>> from softpotato.analytical.benchmark import run_benchmark_suite
>>> suite = run_benchmark_suite(solvers=["crank_nicolson", "scipy_ivp"], cases=["fourier_decay"])
>>> suite.all_passed
True
>>> "crank_nicolson" in suite.summary()
True

Source code in src/softpotato/analytical/benchmark.py

def run_benchmark_suite(
    solvers: Sequence[BaseSolver | str] | None = None,
    cases: (
        Sequence[
            Literal[
                "fourier_decay",
                "cottrell_step",
                "cottrell",
                "spherical_cottrell",
                "mass_conservation",
            ]
        ]
        | None
    ) = None,
    **kwargs: Any,
) -> BenchmarkSuiteResult:
    """Execute a comprehensive benchmark suite across multiple solvers and cases.

    Parameters
    ----------
    solvers : sequence of (BaseSolver or str), optional
        List of solvers to evaluate. Defaults to unconditionally stable solvers:
        ``('crank_nicolson', 'scipy_ivp', 'implicit')``.
    cases : sequence of str, optional
        List of benchmark case names to run. Defaults to all standard cases:
        ``('fourier_decay', 'cottrell_step', 'mass_conservation')``.
    **kwargs : Any
        Additional keyword arguments forwarded to benchmark runs.

    Returns
    -------
    BenchmarkSuiteResult
        Aggregated outcomes of all benchmark runs with summary reporting.

    Examples
    --------
    >>> from softpotato.analytical.benchmark import run_benchmark_suite
    >>> suite = run_benchmark_suite(solvers=["crank_nicolson", "scipy_ivp"], cases=["fourier_decay"])
    >>> suite.all_passed
    True
    >>> "crank_nicolson" in suite.summary()
    True
    """
    eval_solvers: Sequence[BaseSolver | str] = (
        ("crank_nicolson", "scipy_ivp", "implicit") if solvers is None else solvers
    )

    eval_cases: Sequence[
        Literal[
            "fourier_decay",
            "cottrell_step",
            "cottrell",
            "spherical_cottrell",
            "mass_conservation",
        ]
    ] = (
        ("fourier_decay", "cottrell_step", "mass_conservation")
        if cases is None
        else cases
    )

    all_results: list[BenchmarkCaseResult] = []
    t_start = time.perf_counter()

    for s in eval_solvers:
        for c in eval_cases:
            res = run_benchmark(c, s, **kwargs)
            all_results.append(res)

    total_time = time.perf_counter() - t_start
    all_passed = bool(all(res.passed for res in all_results))

    return BenchmarkSuiteResult(
        cases=all_results,
        total_time=total_time,
        all_passed=all_passed,
    )

run_grid_efficiency_benchmark

run_grid_efficiency_benchmark(case: Literal['cottrell', 'cottrell_step', 'spherical_cottrell', 'cyclic_voltammetry'] = 'cottrell', solver: BaseSolver | str = 'scipy_ivp', *, expanding_points: int = 60, uniform_points_list: Sequence[int] = (60, 300, 600), gamma: float = 1.08, beta: float = 4.0, include_exponential: bool = True, D: float = 1e-09, c_bulk: float = 1.0, n_electrons: int = 1, area: float = 1.0, r0: float = 5e-05, t_span: tuple[float, float] = (0.02, 1.0), n_times: int = 25, scan_rate: float = 0.1, domain_factor: float = 6.0, target_rel_error: float = 0.001, **solver_kwargs: Any) -> GridComparisonResult

Evaluate runtime, memory scaling, and \(L_2\) error norms across uniform vs. expanding meshes.

Compares GeometricExpandingMesh (and optionally ExponentialMesh) against Uniform1DMesh at matched and refined spatial point counts (\(5\times\) to \(10\times\) more nodes) for planar Cottrell chronoamperometry, spherical Cottrell chronoamperometry, or reversible cyclic voltammetry.

Parameters:

Name	Type	Description	Default
case	('cottrell', 'cottrell_step', 'spherical_cottrell', 'cyclic_voltammetry')	Electrochemical benchmark problem to run.	'cottrell'
solver	BaseSolver or str	Numerical solver instance or registered string name.	"scipy_ivp"
expanding_points	int	Number of spatial nodes \(N_x\) for the expanding mesh(es).	60
uniform_points_list	sequence of int	Sequence of node counts \(N_x\) for the uniform meshes.	(60, 300, 600)
gamma	float	Expansion ratio \(\gamma > 1\) for GeometricExpandingMesh.	1.08
beta	float	Stretching parameter \(\beta > 0\) for ExponentialMesh.	4.0
include_exponential	bool	Whether to also benchmark an ExponentialMesh of size expanding_points.	True
D	float	Diffusion coefficient \(D > 0\) (\(\text{m}^2/\text{s}\)).	1e-9
c_bulk	float	Bulk concentration \(c^* > 0\) (\(\text{mol}/\text{m}^3\)).	1.0
n_electrons	int	Number of electrons transferred \(n \ge 1\).	1
area	float	Planar electrode area \(A\) (\(\text{m}^2\)).	1.0
r0	float	Spherical electrode radius \(r_0 > 0\) (\(\text{m}\)) when case='spherical_cottrell'.	50e-6
t_span	tuple[float, float]	Evaluation time interval \((t_{\min}, t_{\max})\) for Cottrell benchmarks.	(0.02, 1.0)
n_times	int	Number of time evaluation points.	25
scan_rate	float	Potential sweep rate \(v > 0\) (\(\text{V}/\text{s}\)) when case='cyclic_voltammetry'.	0.1
domain_factor	float	Diffusion layer thickness multiplier \(\delta = \text{domain\_factor} \sqrt{D t_{\max}}\).	6.0
target_rel_error	float	Target maximum relative error threshold (default \(10^{-3} = 0.1\%\)) used to determine passed.	1e-3
**solver_kwargs	Any	Additional keyword arguments forwarded to the solver.	{}

Returns:

Type	Description
GridComparisonResult	Structured comparison of runtime, memory footprint, \(L_2\) error norm, and relative error across uniform and expanding meshes.

Source code in src/softpotato/analytical/benchmark.py

def run_grid_efficiency_benchmark(
    case: Literal[
        "cottrell", "cottrell_step", "spherical_cottrell", "cyclic_voltammetry"
    ] = "cottrell",
    solver: BaseSolver | str = "scipy_ivp",
    *,
    expanding_points: int = 60,
    uniform_points_list: Sequence[int] = (60, 300, 600),
    gamma: float = 1.08,
    beta: float = 4.0,
    include_exponential: bool = True,
    D: float = 1e-9,
    c_bulk: float = 1.0,
    n_electrons: int = 1,
    area: float = 1.0,
    r0: float = 50e-6,
    t_span: tuple[float, float] = (0.02, 1.0),
    n_times: int = 25,
    scan_rate: float = 0.1,
    domain_factor: float = 6.0,
    target_rel_error: float = 1e-3,
    **solver_kwargs: Any,
) -> GridComparisonResult:
    r"""Evaluate runtime, memory scaling, and :math:`L_2` error norms across uniform vs. expanding meshes.

    Compares :class:`~softpotato.grid.GeometricExpandingMesh` (and optionally
    :class:`~softpotato.grid.ExponentialMesh`) against :class:`~softpotato.grid.Uniform1DMesh`
    at matched and refined spatial point counts (:math:`5\times` to :math:`10\times` more nodes)
    for planar Cottrell chronoamperometry, spherical Cottrell chronoamperometry, or
    reversible cyclic voltammetry.

    Parameters
    ----------
    case : {'cottrell', 'cottrell_step', 'spherical_cottrell', 'cyclic_voltammetry'}, default 'cottrell'
        Electrochemical benchmark problem to run.
    solver : BaseSolver or str, default "scipy_ivp"
        Numerical solver instance or registered string name.
    expanding_points : int, default 60
        Number of spatial nodes :math:`N_x` for the expanding mesh(es).
    uniform_points_list : sequence of int, default (60, 300, 600)
        Sequence of node counts :math:`N_x` for the uniform meshes.
    gamma : float, default 1.08
        Expansion ratio :math:`\gamma > 1` for :class:`~softpotato.grid.GeometricExpandingMesh`.
    beta : float, default 4.0
        Stretching parameter :math:`\beta > 0` for :class:`~softpotato.grid.ExponentialMesh`.
    include_exponential : bool, default True
        Whether to also benchmark an :class:`~softpotato.grid.ExponentialMesh` of size
        ``expanding_points``.
    D : float, default 1e-9
        Diffusion coefficient :math:`D > 0` (:math:`\text{m}^2/\text{s}`).
    c_bulk : float, default 1.0
        Bulk concentration :math:`c^* > 0` (:math:`\text{mol}/\text{m}^3`).
    n_electrons : int, default 1
        Number of electrons transferred :math:`n \ge 1`.
    area : float, default 1.0
        Planar electrode area :math:`A` (:math:`\text{m}^2`).
    r0 : float, default 50e-6
        Spherical electrode radius :math:`r_0 > 0` (:math:`\text{m}`) when ``case='spherical_cottrell'``.
    t_span : tuple[float, float], default (0.02, 1.0)
        Evaluation time interval :math:`(t_{\min}, t_{\max})` for Cottrell benchmarks.
    n_times : int, default 25
        Number of time evaluation points.
    scan_rate : float, default 0.1
        Potential sweep rate :math:`v > 0` (:math:`\text{V}/\text{s}`) when ``case='cyclic_voltammetry'``.
    domain_factor : float, default 6.0
        Diffusion layer thickness multiplier :math:`\delta = \text{domain\_factor} \sqrt{D t_{\max}}`.
    target_rel_error : float, default 1e-3
        Target maximum relative error threshold (default :math:`10^{-3} = 0.1\%`) used
        to determine ``passed``.
    **solver_kwargs : Any
        Additional keyword arguments forwarded to the solver.

    Returns
    -------
    GridComparisonResult
        Structured comparison of runtime, memory footprint, :math:`L_2` error norm,
        and relative error across uniform and expanding meshes.
    """
    solver_inst, solver_name = _resolve_solver_instance(solver)
    case_norm = case.lower().strip()
    if case_norm == "cottrell_step":
        case_norm = "cottrell"
    if case_norm not in ("cottrell", "spherical_cottrell", "cyclic_voltammetry"):
        raise ValueError(
            f"Unknown grid efficiency benchmark case '{case}'. "
            f"Expected 'cottrell', 'spherical_cottrell', or 'cyclic_voltammetry'."
        )

    kwargs = dict(solver_kwargs)
    if solver_name in ("scipy_ivp", "solve_ivp"):
        kwargs.setdefault("rtol", 1e-8)
        kwargs.setdefault("atol", 1e-10)

    x_min = float(r0) if case_norm == "spherical_cottrell" else 0.0
    if case_norm == "cyclic_voltammetry":
        e_start = 0.3
        e_vertex = -0.3
        e0 = 0.0
        t_half = abs(e_start - e_vertex) / scan_rate
        t_max = 2.0 * t_half
        t_eval = np.linspace(0.05, t_max, max(n_times, 120))
    else:
        t_start, t_max = t_span
        if t_start <= 0.0 or t_start >= t_max:
            raise ValueError(f"t_span must satisfy 0 < t_start < t_max, got {t_span}.")
        t_eval = np.linspace(t_start, t_max, n_times)

    # Build mesh list: GeometricExpandingMesh, optional ExponentialMesh, and Uniform1DMesh sequence
    meshes_to_run: list[tuple[str, str, Mesh1D]] = [
        (
            f"GeometricExpandingMesh (N={expanding_points}, gamma={gamma})",
            "geometric",
            GeometricExpandingMesh.from_diffusion_layer(
                D=D,
                t_max=t_max,
                x_min=x_min,
                factor=domain_factor,
                n_points=expanding_points,
                gamma=gamma,
            ),
        )
    ]
    if include_exponential:
        meshes_to_run.append(
            (
                f"ExponentialMesh (N={expanding_points}, beta={beta})",
                "exponential",
                ExponentialMesh.from_diffusion_layer(
                    D=D,
                    t_max=t_max,
                    x_min=x_min,
                    factor=domain_factor,
                    n_points=expanding_points,
                    beta=beta,
                ),
            )
        )
    for n_uni in uniform_points_list:
        meshes_to_run.append(
            (
                f"Uniform1DMesh (N={n_uni})",
                "uniform",
                Uniform1DMesh.from_diffusion_layer(
                    D=D,
                    t_max=t_max,
                    x_min=x_min,
                    factor=domain_factor,
                    n_points=int(n_uni),
                ),
            )
        )

    entries: list[GridComparisonEntry] = []

    if case_norm in ("cottrell", "spherical_cottrell"):
        geom: Any = (
            SphericalGeometry(r0=x_min)
            if case_norm == "spherical_cottrell"
            else "planar"
        )
        eff_area = geom.area if case_norm == "spherical_cottrell" else area
        i_exact = (
            cottrell_spherical(
                t_eval, r=x_min, n=n_electrons, D=D, c_bulk=c_bulk, area=eff_area
            )
            if case_norm == "spherical_cottrell"
            else cottrell(t_eval, n=n_electrons, D=D, c_bulk=c_bulk, area=eff_area)
        )

        for label, m_kind, mesh in meshes_to_run:
            problem = DiffusionProblem(
                grid=mesh,
                diffusivity={"A": D},
                boundary_conditions={"A": (DirichletBC(0.0), DirichletBC(c_bulk))},
                initial_conditions={"A": c_bulk},
                geometry=geom,
                n_electrons=n_electrons,
            )
            t0 = time.perf_counter()
            res = solver_inst.solve(
                problem, t_span=(0.0, t_max), t_eval=t_eval, **kwargs
            )
            elapsed = time.perf_counter() - t0

            c_num = res["A"]
            c_exact = np.empty_like(c_num)
            for k, t_k in enumerate(t_eval):
                if case_norm == "spherical_cottrell":
                    eta = (mesh.nodes - x_min) / (2.0 * math.sqrt(D * float(t_k)))
                    c_exact[k, :] = c_bulk * (1.0 - (x_min / mesh.nodes) * erfc(eta))
                else:
                    c_exact[k, :] = step_concentration_profile(
                        mesh.nodes, t=float(t_k), D=D, c_bulk=c_bulk
                    )

            conc_metrics = compute_error_metrics(c_num, c_exact)
            i_num = np.abs(res.currents["A"]) * (
                eff_area if case_norm == "cottrell" else 1.0
            )
            curr_metrics = compute_error_metrics(i_num, i_exact)

            entries.append(
                GridComparisonEntry(
                    label=label,
                    mesh_kind=m_kind,
                    n_points=mesh.n_points,
                    dx_min=mesh.dx_min,
                    dx_max=mesh.dx_max,
                    execution_time=elapsed,
                    memory_bytes=int(c_num.nbytes),
                    current_metrics=curr_metrics,
                    concentration_metrics=conc_metrics,
                )
            )
    else:
        # Reversible Cyclic Voltammetry (O + n e- <-> R)
        e_start = 0.3
        e_vertex = -0.3
        e0 = 0.0
        t_half = abs(e_start - e_vertex) / scan_rate

        def _e_waveform(t: float) -> float:
            if t <= t_half:
                return e_start - scan_rate * t
            return e_vertex + scan_rate * (t - t_half)

        def _c_surf_ox(t: float) -> float:
            c_ox, _ = nernst_equilibrium_concentrations(
                _e_waveform(t), E0=e0, c_total=c_bulk, n=n_electrons
            )
            return float(c_ox)

        ip_theory = abs(
            float(
                randles_sevcik(scan_rate, n=n_electrons, area=area, D=D, c_bulk=c_bulk)
            )
        )

        # High-resolution expanding reference solution for full waveform L_2 comparison
        ref_mesh = GeometricExpandingMesh.from_diffusion_layer(
            D=D, t_max=t_max, factor=domain_factor, n_points=250, gamma=1.025
        )
        ref_prob = DiffusionProblem(
            grid=ref_mesh,
            diffusivity={"Ox": D},
            boundary_conditions={"Ox": (DirichletBC(_c_surf_ox), DirichletBC(c_bulk))},
            initial_conditions={"Ox": c_bulk},
            n_electrons=n_electrons,
        )
        ref_res = solver_inst.solve(
            ref_prob, t_span=(0.0, t_max), t_eval=t_eval, **kwargs
        )
        i_ref = ref_res.currents["Ox"] * area

        for label, m_kind, mesh in meshes_to_run:
            problem = DiffusionProblem(
                grid=mesh,
                diffusivity={"Ox": D},
                boundary_conditions={
                    "Ox": (DirichletBC(_c_surf_ox), DirichletBC(c_bulk))
                },
                initial_conditions={"Ox": c_bulk},
                n_electrons=n_electrons,
            )
            t0 = time.perf_counter()
            res = solver_inst.solve(
                problem, t_span=(0.0, t_max), t_eval=t_eval, **kwargs
            )
            elapsed = time.perf_counter() - t0

            c_num = res["Ox"]
            i_num = res.currents["Ox"] * area
            ip_num = float(np.max(np.abs(i_num)))
            ip_rel_err = abs(ip_num - ip_theory) / ip_theory
            curr_metrics = compute_error_metrics(i_num, i_ref)

            entries.append(
                GridComparisonEntry(
                    label=label,
                    mesh_kind=m_kind,
                    n_points=mesh.n_points,
                    dx_min=mesh.dx_min,
                    dx_max=mesh.dx_max,
                    execution_time=elapsed,
                    memory_bytes=int(c_num.nbytes),
                    current_metrics=curr_metrics,
                    concentration_metrics=None,
                    details={
                        "peak_current_sim": ip_num,
                        "peak_current_theory": ip_theory,
                        "peak_current_rel_error": ip_rel_err,
                    },
                )
            )

    geom_entry = entries[0]
    finest_uni_entry = max(
        (e for e in entries if e.mesh_kind == "uniform"),
        key=lambda e: e.n_points,
    )

    point_reduction = float(finest_uni_entry.n_points) / float(geom_entry.n_points)
    speedup = float(finest_uni_entry.execution_time) / max(
        1e-12, float(geom_entry.execution_time)
    )
    mem_reduction = float(finest_uni_entry.memory_bytes) / float(
        geom_entry.memory_bytes
    )

    if case_norm == "cyclic_voltammetry":
        geom_err = float(geom_entry.details["peak_current_rel_error"])
    else:
        geom_err = geom_entry.max_rel_error

    passed = bool(geom_err <= target_rel_error and point_reduction >= 5.0)

    return GridComparisonResult(
        case_name=case_norm,
        solver_name=solver_name,
        entries=entries,
        point_reduction_factor=point_reduction,
        speedup_factor=speedup,
        memory_reduction_factor=mem_reduction,
        passed=passed,
    )