Skip to content

Grid Module (softpotato.grid)

The softpotato.grid subpackage provides 1D spatial discretization meshes and grid generators for numerical PDE simulations, including uniform, geometrically expanding, and exponential/conformal-mapping meshes, as well as automatic domain sizing from the semi-infinite diffusion layer thickness \(\delta = \text{factor} \cdot \sqrt{D \, t_{\max}}\).

Why Non-Uniform Spatial Discretization?

In electrochemical simulations, concentration perturbations and diffusive fluxes are concentrated within a narrow diffusion boundary layer adjacent to the working electrode surface (\(x = x_{\min}\)), while the bulk solution (\(\nabla^2 c \approx 0\)) experiences negligible gradients.

  • Uniform grids (Uniform1DMesh) distribute spatial nodes equally across the entire domain, requiring hundreds or thousands of nodes to resolve steep interfacial gradients \(\left.\partial c / \partial x\right|_{x_{\min}}\) while oversampling the quiescent bulk.
  • Expanding and stretched grids (GeometricExpandingMesh and ExponentialMesh) cluster fine intervals (\(\Delta x_0 \ll \delta\)) near active boundaries and smoothly coarsen into the bulk solution, achieving 5–20× faster simulation runtimes and 80–90% fewer spatial nodes at equal or superior flux accuracy while preserving second-order spatial convergence (\(\mathcal{O}(\Delta x^2)\)).

Mesh Classes and Generators

Class / Function Spacing Profile Key Parameters Description
Mesh1D Abstract nodes, dx Abstract base class implementing the 1D mesh interface and NumPy array protocol (__array__, __array_ufunc__, slicing, arithmetic).
Uniform1DMesh \(\Delta x_j = \Delta x\) x_min, x_max, n_points or dx Equidistant 1D spatial mesh with constant node spacing \(\Delta x = (x_{\max} - x_{\min}) / (N - 1)\).
GeometricExpandingMesh \(\Delta x_j = \Delta x_0 \gamma^j\) x_min, x_max, any two of (n_points, gamma, dx_min), two_sided Non-uniform mesh with constant interval expansion ratio \(\gamma > 1\).
ExponentialMesh \(x(\xi)\) via exponential, tanh, or sinh x_min, x_max, n_points, beta, two_sided, mapping Smooth non-uniform mesh using exponential or hyperbolic coordinate stretching of uniform \(\xi \in [0, 1]\).
create_diffusion_mesh 'geometric', 'exponential', or 'uniform' D, t_max, kind, factor High-level factory constructing a 1D mesh automatically sized to \(\delta = \text{factor} \cdot \sqrt{D \, t_{\max}}\).
diffusion_layer_thickness Scalar \(\delta\) D, t_max, factor Estimates the characteristic semi-infinite diffusion layer thickness \(\delta = \text{factor} \cdot \sqrt{D \, t_{\max}}\).

Mathematical Formulations

1. Physics-Informed Domain Sizing

For semi-infinite diffusion over duration \(t_{\max}\) with diffusivity \(D\), the characteristic penetration thickness is:

\[ \delta = \text{factor} \cdot \sqrt{D \, t_{\max}}, \quad x_{\max} = x_{\min} + \delta \]

At the default factor=6.0, the complementary error function decays to \(\operatorname{erfc}(6 / 2) \approx 2.2 \times 10^{-5}\), ensuring outer boundary perturbations remain negligible.

2. Geometrically Expanding Meshes (GeometricExpandingMesh)

Consecutive cell intervals expand by a constant geometric factor \(\gamma > 1\):

\[ \Delta x_j = \Delta x_0 \cdot \gamma^j \quad (j = 0, 1, \dots, N - 2), \quad x_{j+1} = x_j + \Delta x_j \]

where the total domain length satisfies the geometric series:

\[ L = x_{\max} - x_{\min} = \Delta x_0 \sum_{j=0}^{N-2} \gamma^j = \Delta x_0 \frac{\gamma^{N-1} - 1}{\gamma - 1} \]

GeometricExpandingMesh supports three flexible parameterizations by specifying any two of n_points (\(N\)), gamma (\(\gamma\)), and dx_min (\(\Delta x_0\)): 1. (n_points, gamma): Solves analytically for the initial boundary step \(\Delta x_0\). 2. (dx_min, gamma): Computes the required node count \(N\) so the boundary step is at most dx_min. 3. (n_points, dx_min): Solves numerically via Brent's root-finding method (scipy.optimize.brentq) for the exact expansion ratio \(\gamma > 1\).

3. Exponential & Conformal-Mapping Meshes (ExponentialMesh)

Maps uniform computational coordinates \(\xi_j = j / (N - 1) \in [0, 1]\) to physical coordinates \(x(\xi) \in [x_{\min}, x_{\max}]\) parameterized by a dimensionless stretching factor \(\beta > 0\):

  • Exponential mapping (mapping="exponential", default):

$$ x(\xi) = x_{\min} + (x_{\max} - x_{\min}) \frac{\exp(\beta \xi) - 1}{\exp(\beta) - 1} $$

  • Hyperbolic tangent mapping (mapping="tanh" or ExponentialMesh.from_tanh):

$$ x(\xi) = x_{\min} + (x_{\max} - x_{\min}) \left(1 - \frac{\tanh(\beta (1 - \xi))}{\tanh(\beta)}\right) $$

  • Hyperbolic sine mapping (mapping="sinh" or ExponentialMesh.from_sinh):

$$ x(\xi) = x_{\min} + (x_{\max} - x_{\min}) \frac{\sinh(\beta \xi)}{\sinh(\beta)} $$

4. Two-Sided Symmetric Clustering (two_sided=True)

Both GeometricExpandingMesh and ExponentialMesh support two_sided=True, clustering fine spatial intervals symmetrically at both \(x_{\min}\) and \(x_{\max}\) while coarsening toward the domain midpoint \(x_{\text{mid}} = \frac{1}{2}(x_{\min} + x_{\max})\). This is ideal for thin-layer electrochemical cells, SECM tip–substrate gaps (generation–collection / feedback), and dual-working-electrode (electrodes="both") simulations.

User Guide & Examples

1. Automatic Diffusion-Layer Mesh Generation

Use create_diffusion_mesh or the .from_diffusion_layer() classmethods to automatically size the domain to \(\delta = 6\sqrt{D t_{\max}}\):

from softpotato.grid import (
    ExponentialMesh,
    GeometricExpandingMesh,
    Uniform1DMesh,
    create_diffusion_mesh,
)

D = 1e-9  # m^2/s
t_max = 1.0  # s

# High-level factory (defaults to kind="geometric", factor=6.0)
mesh_geom = create_diffusion_mesh(
    D=D, t_max=t_max, kind="geometric", n_points=60, gamma=1.08
)

# Or construct directly via classmethods
mesh_exp = ExponentialMesh.from_diffusion_layer(
    D=D, t_max=t_max, n_points=60, beta=3.5, mapping="sinh"
)
mesh_uni = Uniform1DMesh.from_diffusion_layer(D=D, t_max=t_max, n_points=300)

print(f"Domain bounds: {mesh_geom.bounds}")
print(f"Geometric dx_min = {mesh_geom.dx_min:.3e} m, dx_max = {mesh_geom.dx_max:.3e} m")

2. Three Parameterizations of GeometricExpandingMesh

Specify any two of n_points, gamma, and dx_min:

from softpotato.grid import GeometricExpandingMesh

# 1. Specify (n_points, gamma) -> computes dx_min
m1 = GeometricExpandingMesh(0.0, 1e-3, n_points=51, gamma=1.08)

# 2. Specify (dx_min, gamma) -> determines required n_points
m2 = GeometricExpandingMesh(0.0, 1e-3, dx_min=1e-6, gamma=1.08)

# 3. Specify (n_points, dx_min) -> solves for gamma > 1
m3 = GeometricExpandingMesh(0.0, 1e-3, n_points=51, dx_min=1e-6)

# Two-sided symmetric mesh for thin-layer / SECM cells
m_thin = GeometricExpandingMesh(0.0, 50e-6, n_points=61, gamma=1.06, two_sided=True)

3. Coupling Expanding Grids with DiffusionProblem and Solvers

All Mesh1D subclasses implement the NumPy array protocol and integrate seamlessly with DiffusionProblem and all Soft Potato solvers (scipy_ivp, crank_nicolson, implicit, explicit) across planar, spherical, and cylindrical geometries:

from softpotato.grid import create_diffusion_mesh
from softpotato.solver import DiffusionProblem, DirichletBC, get_solver

D = 1e-9
t_max = 1.0
mesh = create_diffusion_mesh(
    D=D, t_max=t_max, kind="geometric", n_points=60, gamma=1.08
)

problem = DiffusionProblem(
    grid=mesh,
    diffusivity={"Ox": D},
    boundary_conditions={"Ox": (DirichletBC(0.0), DirichletBC(1.0))},
    initial_conditions={"Ox": 1.0},
    geometry="planar",
    n_electrons=1,
)

solver = get_solver("scipy_ivp")
result = solver.solve(problem, t_span=(0.0, t_max))
print(result.grid is mesh)  # True: preserves Mesh1D reference
print(result.currents["Ox"].shape)  # (100,) Faradaic current transient in Amperes

Module Contents

grid

Spatial discretization meshes and grid generators for numerical simulation.

Mesh1D

Bases: ABC

Abstract base class representing a 1D spatial discretization mesh.

Encapsulates 1D coordinate points \(x_j\), cell intervals \(\Delta x_j = x_{j+1} - x_j\), spatial domain boundaries, and grid uniformity diagnostics.

Provides full interoperability with the NumPy array protocol (slicing, indexing, mathematical operators, and universal functions).

Attributes:

Name Type Description
nodes ndarray

1D array of spatial node positions.
dx ndarray

Interval spacings between consecutive nodes: \(\Delta x_j = x_{j+1} - x_j\).
n_points int

Total number of spatial grid nodes \(N\).
bounds tuple[float, float]

Spatial domain boundaries \((x_{\min}, x_{\max})\).
is_uniform bool

Whether the grid spacing is uniform within numerical tolerance.
Source code in src/softpotato/grid/base.py 
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
class Mesh1D(ABC):
    r"""Abstract base class representing a 1D spatial discretization mesh.

    Encapsulates 1D coordinate points :math:`x_j`, cell intervals
    :math:`\Delta x_j = x_{j+1} - x_j`, spatial domain boundaries,
    and grid uniformity diagnostics.

    Provides full interoperability with the NumPy array protocol
    (slicing, indexing, mathematical operators, and universal functions).

    Attributes
    ----------
    nodes : np.ndarray
        1D array of spatial node positions.
    dx : np.ndarray
        Interval spacings between consecutive nodes: :math:`\Delta x_j = x_{j+1} - x_j`.
    n_points : int
        Total number of spatial grid nodes :math:`N`.
    bounds : tuple[float, float]
        Spatial domain boundaries :math:`(x_{\min}, x_{\max})`.
    is_uniform : bool
        Whether the grid spacing is uniform within numerical tolerance.
    """

    @property
    @abstractmethod
    def nodes(self) -> np.ndarray:
        """1D array of spatial node positions."""
        raise NotImplementedError

    @property
    @abstractmethod
    def dx(self) -> np.ndarray:
        """Interval spacings between consecutive nodes: dx[j] = nodes[j+1] - nodes[j]."""
        raise NotImplementedError

    @property
    def n_points(self) -> int:
        """Total number of spatial nodes N."""
        return len(self.nodes)

    @property
    def shape(self) -> tuple[int, ...]:
        """Array shape of the spatial nodes `(N,)`."""
        return self.nodes.shape

    @property
    def ndim(self) -> int:
        """Number of array dimensions (`1`)."""
        return 1

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

    @property
    def x_min(self) -> float:
        """Spatial domain minimum boundary coordinate."""
        return self.bounds[0]

    @property
    def x_max(self) -> float:
        """Spatial domain maximum boundary coordinate."""
        return self.bounds[1]

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

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

    @property
    def is_uniform(self) -> bool:
        """Whether the grid spacing is uniform within numerical tolerance."""
        if len(self.dx) <= 1:
            return True
        return bool(np.allclose(self.dx, self.dx[0], rtol=1e-5, atol=0.0))

    def __array__(
        self, dtype: np.typing.DTypeLike | None = None, copy: bool | None = None
    ) -> np.ndarray:
        """NumPy array protocol support: np.asarray(mesh) yields mesh.nodes."""
        if copy:
            return np.array(self.nodes, dtype=dtype, copy=True)
        return np.asarray(self.nodes, dtype=dtype)

    def __len__(self) -> int:
        """Return the number of spatial nodes."""
        return self.n_points

    def __getitem__(self, item: Any) -> Any:
        """Index or slice the spatial nodes array."""
        return self.nodes[item]

    def __iter__(self) -> Iterator[Any]:
        """Iterate over the spatial node coordinates."""
        return iter(self.nodes)

    def __array_ufunc__(
        self, ufunc: np.ufunc, method: str, *inputs: Any, **kwargs: Any
    ) -> Any:
        """Delegate NumPy ufunc evaluations to the underlying nodes array."""
        converted_inputs = tuple(
            np.asarray(x) if isinstance(x, Mesh1D) else x for x in inputs
        )
        return getattr(ufunc, method)(*converted_inputs, **kwargs)

    def __add__(self, other: Any) -> np.ndarray:
        return self.nodes + (other.nodes if isinstance(other, Mesh1D) else other)

    def __radd__(self, other: Any) -> np.ndarray:
        return (other.nodes if isinstance(other, Mesh1D) else other) + self.nodes

    def __sub__(self, other: Any) -> np.ndarray:
        return self.nodes - (other.nodes if isinstance(other, Mesh1D) else other)

    def __rsub__(self, other: Any) -> np.ndarray:
        return (other.nodes if isinstance(other, Mesh1D) else other) - self.nodes

    def __mul__(self, other: Any) -> np.ndarray:
        return self.nodes * (other.nodes if isinstance(other, Mesh1D) else other)

    def __rmul__(self, other: Any) -> np.ndarray:
        return (other.nodes if isinstance(other, Mesh1D) else other) * self.nodes

    def __truediv__(self, other: Any) -> np.ndarray:
        return self.nodes / (other.nodes if isinstance(other, Mesh1D) else other)

    def __rtruediv__(self, other: Any) -> np.ndarray:
        return (other.nodes if isinstance(other, Mesh1D) else other) / self.nodes

    def __pow__(self, other: Any) -> np.ndarray:
        return self.nodes ** (other.nodes if isinstance(other, Mesh1D) else other)

    def __neg__(self) -> np.ndarray:
        return -self.nodes

    def __pos__(self) -> np.ndarray:
        return +self.nodes

    def __abs__(self) -> np.ndarray:
        return np.abs(self.nodes)

    def __lt__(self, other: Any) -> np.ndarray:
        return self.nodes < (other.nodes if isinstance(other, Mesh1D) else other)

    def __le__(self, other: Any) -> np.ndarray:
        return self.nodes <= (other.nodes if isinstance(other, Mesh1D) else other)

    def __gt__(self, other: Any) -> np.ndarray:
        return self.nodes > (other.nodes if isinstance(other, Mesh1D) else other)

    def __ge__(self, other: Any) -> np.ndarray:
        return self.nodes >= (other.nodes if isinstance(other, Mesh1D) else other)

    def __eq__(self, other: object) -> bool:
        if isinstance(other, Mesh1D):
            return bool(
                self.n_points == other.n_points
                and np.allclose(self.nodes, other.nodes, rtol=1e-12, atol=1e-12)
            )
        if isinstance(other, (np.ndarray, list, tuple)):
            arr = np.asarray(other)
            if arr.shape != self.nodes.shape:
                return False
            return bool(np.allclose(self.nodes, arr, rtol=1e-12, atol=1e-12))
        return False

    def __repr__(self) -> str:
        return (
            f"<{self.__class__.__name__} "
            f"n_points={self.n_points} "
            f"bounds=({self.bounds[0]:.4g}, {self.bounds[1]:.4g})>"
        )

nodes abstractmethod property 

nodes: ndarray

1D array of spatial node positions.

dx abstractmethod property 

dx: ndarray

Interval spacings between consecutive nodes: dx[j] = nodes[j+1] - nodes[j].

n_points property 

n_points: int

Total number of spatial nodes N.

shape property 

shape: tuple[int, ...]

Array shape of the spatial nodes (N,).

ndim property 

ndim: int

Number of array dimensions (1).

bounds property 

bounds: tuple[float, float]

Spatial domain boundaries (x_min, x_max).

x_min property 

x_min: float

Spatial domain minimum boundary coordinate.

x_max property 

x_max: float

Spatial domain maximum boundary coordinate.

dx_min property 

dx_min: float

Minimum interval spacing across the mesh.

dx_max property 

dx_max: float

Maximum interval spacing across the mesh.

is_uniform property 

is_uniform: bool

Whether the grid spacing is uniform within numerical tolerance.

__array__ 

__array__(dtype: DTypeLike | None = None, copy: bool | None = None) -> np.ndarray

NumPy array protocol support: np.asarray(mesh) yields mesh.nodes.

Source code in src/softpotato/grid/base.py 
149
150
151
152
153
154
155
def __array__(
    self, dtype: np.typing.DTypeLike | None = None, copy: bool | None = None
) -> np.ndarray:
    """NumPy array protocol support: np.asarray(mesh) yields mesh.nodes."""
    if copy:
        return np.array(self.nodes, dtype=dtype, copy=True)
    return np.asarray(self.nodes, dtype=dtype)

__len__ 

__len__() -> int

Return the number of spatial nodes.

Source code in src/softpotato/grid/base.py 
157
158
159
def __len__(self) -> int:
    """Return the number of spatial nodes."""
    return self.n_points

__getitem__ 

__getitem__(item: Any) -> Any

Index or slice the spatial nodes array.

Source code in src/softpotato/grid/base.py 
161
162
163
def __getitem__(self, item: Any) -> Any:
    """Index or slice the spatial nodes array."""
    return self.nodes[item]

__iter__ 

__iter__() -> Iterator[Any]

Iterate over the spatial node coordinates.

Source code in src/softpotato/grid/base.py 
165
166
167
def __iter__(self) -> Iterator[Any]:
    """Iterate over the spatial node coordinates."""
    return iter(self.nodes)

__array_ufunc__ 

__array_ufunc__(ufunc: ufunc, method: str, *inputs: Any, **kwargs: Any) -> Any

Delegate NumPy ufunc evaluations to the underlying nodes array.

Source code in src/softpotato/grid/base.py 
169
170
171
172
173
174
175
176
def __array_ufunc__(
    self, ufunc: np.ufunc, method: str, *inputs: Any, **kwargs: Any
) -> Any:
    """Delegate NumPy ufunc evaluations to the underlying nodes array."""
    converted_inputs = tuple(
        np.asarray(x) if isinstance(x, Mesh1D) else x for x in inputs
    )
    return getattr(ufunc, method)(*converted_inputs, **kwargs)

ExponentialMesh

Bases: Mesh1D

Exponential and conformal-mapping 1D spatial discretization mesh.

Generates a smooth, non-uniform 1D spatial grid across \([x_{\min}, x_{\max}]\) parameterized by a stretching factor \(\beta > 0\) applied to uniform computational coordinates \(\xi \in [0, 1]\).

For the default "exponential" mapping (one-sided):

\[ x(\xi) = x_{\min} + (x_{\max} - x_{\min}) \frac{\exp(\beta \xi) - 1}{\exp(\beta) - 1}, \quad \xi \in [0, 1] \]

Also supports conformal / boundary-layer coordinate transformations:

  • "tanh" (hyperbolic tangent coordinate clustering):

    \[ x(\xi) = x_{\min} + (x_{\max} - x_{\min}) \left(1 - \frac{\tanh(\beta (1 - \xi))}{\tanh(\beta)}\right) \]

  • "sinh" (hyperbolic sine conformal stretching):

    \[ x(\xi) = x_{\min} + (x_{\max} - x_{\min}) \frac{\sinh(\beta \xi)}{\sinh(\beta)} \]

When two_sided=True, nodes are symmetrically clustered near both boundaries (\(x_{\min}\) and \(x_{\max}\)) while coarsening toward the domain center.

Parameters:

Name Type Description Default
x_min float

Lower boundary of the spatial domain (\(x_{\min}\)).

 required
x_max float

Upper boundary of the spatial domain (\(x_{\max}\)). Must be strictly greater than x_min.

 required
n_points int

Total number of grid nodes \(N\). Must be \(\ge 2\) (or \(\ge 3\) when two_sided=True).

 required
beta float

Dimensionless stretching parameter \(\beta > 0\). Larger values concentrate more nodes near the boundary.

 3.0
two_sided bool

If True, clusters nodes symmetrically at both \(x_{\min}\) and \(x_{\max}\).

 False
mapping str

Coordinate transformation function: "exponential", "tanh", or "sinh".

 "exponential"

Attributes:

Name Type Description
nodes ndarray

1D array of spatial node positions of length \(N\).
dx ndarray

Interval spacings between consecutive nodes of length \(N-1\).
xi ndarray

Uniform computational coordinates \(\xi \in [0, 1]\) of length \(N\).
beta float

Stretching parameter \(\beta > 0\).
two_sided bool

Whether the mesh is symmetrically clustered at both boundaries.
mapping str

Name of the coordinate mapping transformation.
dx_min float

Minimum interval spacing across the mesh.
dx_max float

Maximum interval spacing across the mesh.
n_points int

Total number of spatial nodes \(N\).
bounds tuple[float, float]

Spatial domain boundaries \((x_{\min}, x_{\max})\).
is_uniform bool

Always False for exponential meshes.

Raises:

Type Description
ValueError

If x_max <= x_min, n_points is too small, beta <= 0, beta > 500, or mapping is unknown.
TypeError

If inputs have invalid types.

Examples:

>>> from softpotato.grid import ExponentialMesh
>>> mesh = ExponentialMesh(0.0, 1.0, n_points=51, beta=3.0)
>>> mesh.bounds
(0.0, 1.0)
>>> mesh.is_uniform
False
Source code in src/softpotato/grid/exponential.py 
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
class ExponentialMesh(Mesh1D):
    r"""Exponential and conformal-mapping 1D spatial discretization mesh.

    Generates a smooth, non-uniform 1D spatial grid across :math:`[x_{\min}, x_{\max}]`
    parameterized by a stretching factor :math:`\beta > 0` applied to uniform
    computational coordinates :math:`\xi \in [0, 1]`.

    For the default ``"exponential"`` mapping (one-sided):

    .. math::

        x(\xi) = x_{\min} + (x_{\max} - x_{\min}) \frac{\exp(\beta \xi) - 1}{\exp(\beta) - 1},
        \quad \xi \in [0, 1]

    Also supports conformal / boundary-layer coordinate transformations:

    - ``"tanh"`` (hyperbolic tangent coordinate clustering):

      .. math::

          x(\xi) = x_{\min} + (x_{\max} - x_{\min})
          \left(1 - \frac{\tanh(\beta (1 - \xi))}{\tanh(\beta)}\right)

    - ``"sinh"`` (hyperbolic sine conformal stretching):

      .. math::

          x(\xi) = x_{\min} + (x_{\max} - x_{\min}) \frac{\sinh(\beta \xi)}{\sinh(\beta)}

    When ``two_sided=True``, nodes are symmetrically clustered near both
    boundaries (:math:`x_{\min}` and :math:`x_{\max}`) while coarsening toward
    the domain center.

    Parameters
    ----------
    x_min : float
        Lower boundary of the spatial domain (:math:`x_{\min}`).
    x_max : float
        Upper boundary of the spatial domain (:math:`x_{\max}`).
        Must be strictly greater than ``x_min``.
    n_points : int
        Total number of grid nodes :math:`N`. Must be :math:`\ge 2` (or :math:`\ge 3`
        when ``two_sided=True``).
    beta : float, default 3.0
        Dimensionless stretching parameter :math:`\beta > 0`. Larger values
        concentrate more nodes near the boundary.
    two_sided : bool, default False
        If ``True``, clusters nodes symmetrically at both :math:`x_{\min}` and
        :math:`x_{\max}`.
    mapping : str, default "exponential"
        Coordinate transformation function: ``"exponential"``, ``"tanh"``, or ``"sinh"``.

    Attributes
    ----------
    nodes : np.ndarray
        1D array of spatial node positions of length :math:`N`.
    dx : np.ndarray
        Interval spacings between consecutive nodes of length :math:`N-1`.
    xi : np.ndarray
        Uniform computational coordinates :math:`\xi \in [0, 1]` of length :math:`N`.
    beta : float
        Stretching parameter :math:`\beta > 0`.
    two_sided : bool
        Whether the mesh is symmetrically clustered at both boundaries.
    mapping : str
        Name of the coordinate mapping transformation.
    dx_min : float
        Minimum interval spacing across the mesh.
    dx_max : float
        Maximum interval spacing across the mesh.
    n_points : int
        Total number of spatial nodes :math:`N`.
    bounds : tuple[float, float]
        Spatial domain boundaries :math:`(x_{\min}, x_{\max})`.
    is_uniform : bool
        Always ``False`` for exponential meshes.

    Raises
    ------
    ValueError
        If ``x_max <= x_min``, ``n_points`` is too small, ``beta <= 0``,
        ``beta > 500``, or ``mapping`` is unknown.
    TypeError
        If inputs have invalid types.

    Examples
    --------
    >>> from softpotato.grid import ExponentialMesh
    >>> mesh = ExponentialMesh(0.0, 1.0, n_points=51, beta=3.0)
    >>> mesh.bounds
    (0.0, 1.0)
    >>> mesh.is_uniform
    False
    """

    def __init__(
        self,
        x_min: float,
        x_max: float,
        n_points: int,
        beta: float = 3.0,
        two_sided: bool = False,
        mapping: str = "exponential",
    ) -> None:
        if isinstance(x_min, (bool, str, bytes)) or isinstance(
            x_max, (bool, str, bytes)
        ):
            raise TypeError("x_min and x_max must be real numbers.")
        try:
            x_min_f = float(x_min)
            x_max_f = float(x_max)
        except (TypeError, ValueError) as exc:
            raise TypeError("x_min and x_max must be real numbers.") from exc

        if not (math.isfinite(x_min_f) and math.isfinite(x_max_f)):
            raise ValueError(f"x_min and x_max must be finite, got ({x_min}, {x_max}).")

        if x_max_f <= x_min_f:
            raise ValueError(
                f"x_max ({x_max_f}) must be strictly greater than x_min ({x_min_f})."
            )

        if isinstance(n_points, bool) or not isinstance(n_points, (int, np.integer)):
            raise TypeError(
                f"n_points must be an integer, got {type(n_points).__name__}."
            )
        n_pts = int(n_points)

        if not isinstance(two_sided, bool):
            raise TypeError(
                f"two_sided must be a boolean, got {type(two_sided).__name__}."
            )

        min_pts = 3 if two_sided else 2
        if n_pts < min_pts:
            raise ValueError(f"n_points must be >= {min_pts}, got {n_pts}.")

        if isinstance(beta, (bool, str, bytes)):
            raise TypeError("beta must be a real number.")
        try:
            beta_f = float(beta)
        except (TypeError, ValueError) as exc:
            raise TypeError("beta must be a real number.") from exc

        if not math.isfinite(beta_f):
            raise ValueError(f"beta must be finite, got {beta}.")
        if beta_f <= 0.0:
            raise ValueError(f"beta must be strictly positive, got {beta_f}.")
        if beta_f > 500.0:
            raise ValueError(
                f"beta ({beta_f}) is too large (must be <= 500.0 to avoid overflow)."
            )

        if not isinstance(mapping, str):
            raise TypeError(f"mapping must be a string, got {type(mapping).__name__}.")
        mapping_norm = mapping.strip().lower()
        if mapping_norm not in _VALID_MAPPINGS:
            raise ValueError(
                f"Unknown mapping '{mapping}'. Must be one of {_VALID_MAPPINGS}."
            )

        xi = np.linspace(0.0, 1.0, n_pts, dtype=float)
        eta = self._compute_normalized_coords(xi, beta_f, two_sided, mapping_norm)

        domain_len = x_max_f - x_min_f
        nodes = x_min_f + domain_len * eta
        nodes[0] = x_min_f
        nodes[-1] = x_max_f
        if two_sided and n_pts % 2 == 1:
            nodes[n_pts // 2] = 0.5 * (x_min_f + x_max_f)

        dx_arr = np.diff(nodes)
        xi.flags.writeable = False
        nodes.flags.writeable = False
        dx_arr.flags.writeable = False

        self._xi = xi
        self._nodes = nodes
        self._dx = dx_arr
        self._beta = beta_f
        self._two_sided = two_sided
        self._mapping = mapping_norm

    @staticmethod
    def _compute_normalized_coords(
        xi: np.ndarray,
        beta: float,
        two_sided: bool,
        mapping: str,
    ) -> np.ndarray:
        """Compute normalized coordinates eta in [0, 1] from uniform xi in [0, 1]."""
        if mapping == "exponential":
            denom = math.expm1(beta)
            if not two_sided:
                return np.expm1(beta * xi) / denom
            eta = np.empty_like(xi)
            left = xi <= 0.5
            right = ~left
            eta[left] = 0.5 * np.expm1(2.0 * beta * xi[left]) / denom
            eta[right] = 1.0 - 0.5 * np.expm1(2.0 * beta * (1.0 - xi[right])) / denom
            return eta

        if mapping == "tanh":
            tanh_b = math.tanh(beta)
            if not two_sided:
                return 1.0 - np.tanh(beta * (1.0 - xi)) / tanh_b
            return 0.5 * (1.0 + np.tanh(beta * (2.0 * xi - 1.0)) / tanh_b)

        # mapping == "sinh"
        sinh_b = math.sinh(beta)
        if not two_sided:
            return np.sinh(beta * xi) / sinh_b
        eta = np.empty_like(xi)
        left = xi <= 0.5
        right = ~left
        eta[left] = 0.5 * np.sinh(2.0 * beta * xi[left]) / sinh_b
        eta[right] = 1.0 - 0.5 * np.sinh(2.0 * beta * (1.0 - xi[right])) / sinh_b
        return eta

    @classmethod
    def from_tanh(
        cls,
        x_min: float,
        x_max: float,
        n_points: int,
        beta: float = 2.5,
        two_sided: bool = False,
    ) -> ExponentialMesh:
        """Construct a mesh using hyperbolic tangent (``tanh``) coordinate stretching."""
        return cls(
            x_min=x_min,
            x_max=x_max,
            n_points=n_points,
            beta=beta,
            two_sided=two_sided,
            mapping="tanh",
        )

    @classmethod
    def from_sinh(
        cls,
        x_min: float,
        x_max: float,
        n_points: int,
        beta: float = 3.0,
        two_sided: bool = False,
    ) -> ExponentialMesh:
        """Construct a mesh using hyperbolic sine (``sinh``) conformal stretching."""
        return cls(
            x_min=x_min,
            x_max=x_max,
            n_points=n_points,
            beta=beta,
            two_sided=two_sided,
            mapping="sinh",
        )

    @classmethod
    def from_diffusion_layer(
        cls,
        D: float,
        t_max: float,
        *,
        x_min: float = 0.0,
        factor: float = 6.0,
        n_points: int = 100,
        beta: float = 3.0,
        two_sided: bool = False,
        mapping: str = "exponential",
    ) -> ExponentialMesh:
        r"""Create an exponential mesh sized to the diffusion layer thickness.

        Sets :math:`x_{\max} = x_{\min} + \text{factor} \cdot \sqrt{D \, t_{\max}}`,
        where ``factor`` defaults to :math:`6.0` (:math:`\delta \approx 6 \sqrt{D t_{\max}}`).

        Parameters
        ----------
        D : float
            Diffusion coefficient :math:`D > 0` (:math:`\text{m}^2/\text{s}`).
        t_max : float
            Maximum characteristic simulation time :math:`t_{\max} > 0` (:math:`\text{s}`).
        x_min : float, default 0.0
            Electrode surface coordinate :math:`x_{\min}`.
        factor : float, default 6.0
            Dimensionless multiplier for diffusion layer thickness
            :math:`\delta = \text{factor} \cdot \sqrt{D t_{\max}}`.
        n_points : int, default 100
            Total number of grid nodes :math:`N`.
        beta : float, default 3.0
            Dimensionless stretching parameter :math:`\beta > 0`.
        two_sided : bool, default False
            Whether to cluster nodes symmetrically at both boundaries.
        mapping : str, default "exponential"
            Coordinate transformation function: ``"exponential"``, ``"tanh"``, or ``"sinh"``.

        Returns
        -------
        ExponentialMesh
            Configured exponential mesh spanning :math:`[x_{\min}, x_{\min} + \delta]`.
        """
        delta = diffusion_layer_thickness(D, t_max, factor=factor)
        if isinstance(x_min, bool):
            raise TypeError("x_min must be a real number.")
        try:
            x_min_f = float(x_min)
        except (TypeError, ValueError) as exc:
            raise TypeError("x_min must be a real number.") from exc

        return cls(
            x_min=x_min_f,
            x_max=x_min_f + delta,
            n_points=n_points,
            beta=beta,
            two_sided=two_sided,
            mapping=mapping,
        )

    @property
    def nodes(self) -> np.ndarray:
        """1D array of spatial node positions."""
        return self._nodes

    @property
    def dx(self) -> np.ndarray:
        """Interval spacings between consecutive nodes."""
        return self._dx

    @property
    def xi(self) -> np.ndarray:
        """Uniform computational coordinate array xi in [0, 1]."""
        return self._xi

    @property
    def beta(self) -> float:
        """Stretching parameter beta > 0."""
        return self._beta

    @property
    def two_sided(self) -> bool:
        """Whether the mesh is symmetrically clustered at both boundaries."""
        return self._two_sided

    @property
    def mapping(self) -> str:
        """Coordinate mapping function ('exponential', 'tanh', or 'sinh')."""
        return self._mapping

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

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

    @property
    def is_uniform(self) -> bool:
        """Whether the grid spacing is uniform within numerical tolerance."""
        return False

    def __repr__(self) -> str:
        return (
            f"ExponentialMesh(x_min={self.bounds[0]}, x_max={self.bounds[1]}, "
            f"n_points={self.n_points}, beta={self.beta:.4g}, "
            f"two_sided={self.two_sided}, mapping='{self.mapping}')"
        )

    def __str__(self) -> str:
        return (
            f"ExponentialMesh(bounds=({self.bounds[0]}, {self.bounds[1]}), "
            f"n_points={self.n_points}, beta={self.beta:.4g}, "
            f"two_sided={self.two_sided}, mapping='{self.mapping}')"
        )

nodes property 

nodes: ndarray

1D array of spatial node positions.

dx property 

dx: ndarray

Interval spacings between consecutive nodes.

xi property 

xi: ndarray

Uniform computational coordinate array xi in [0, 1].

beta property 

beta: float

Stretching parameter beta > 0.

two_sided property 

two_sided: bool

Whether the mesh is symmetrically clustered at both boundaries.

mapping property 

mapping: str

Coordinate mapping function ('exponential', 'tanh', or 'sinh').

dx_min property 

dx_min: float

Minimum interval spacing across the mesh.

dx_max property 

dx_max: float

Maximum interval spacing across the mesh.

is_uniform property 

is_uniform: bool

Whether the grid spacing is uniform within numerical tolerance.

from_tanh classmethod 

from_tanh(x_min: float, x_max: float, n_points: int, beta: float = 2.5, two_sided: bool = False) -> ExponentialMesh

Construct a mesh using hyperbolic tangent (tanh) coordinate stretching.

Source code in src/softpotato/grid/exponential.py 
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
@classmethod
def from_tanh(
    cls,
    x_min: float,
    x_max: float,
    n_points: int,
    beta: float = 2.5,
    two_sided: bool = False,
) -> ExponentialMesh:
    """Construct a mesh using hyperbolic tangent (``tanh``) coordinate stretching."""
    return cls(
        x_min=x_min,
        x_max=x_max,
        n_points=n_points,
        beta=beta,
        two_sided=two_sided,
        mapping="tanh",
    )

from_sinh classmethod 

from_sinh(x_min: float, x_max: float, n_points: int, beta: float = 3.0, two_sided: bool = False) -> ExponentialMesh

Construct a mesh using hyperbolic sine (sinh) conformal stretching.

Source code in src/softpotato/grid/exponential.py 
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
@classmethod
def from_sinh(
    cls,
    x_min: float,
    x_max: float,
    n_points: int,
    beta: float = 3.0,
    two_sided: bool = False,
) -> ExponentialMesh:
    """Construct a mesh using hyperbolic sine (``sinh``) conformal stretching."""
    return cls(
        x_min=x_min,
        x_max=x_max,
        n_points=n_points,
        beta=beta,
        two_sided=two_sided,
        mapping="sinh",
    )

from_diffusion_layer classmethod 

from_diffusion_layer(D: float, t_max: float, *, x_min: float = 0.0, factor: float = 6.0, n_points: int = 100, beta: float = 3.0, two_sided: bool = False, mapping: str = 'exponential') -> ExponentialMesh

Create an exponential mesh sized to the diffusion layer thickness.

Sets \(x_{\max} = x_{\min} + \text{factor} \cdot \sqrt{D \, t_{\max}}\), where factor defaults to \(6.0\) (\(\delta \approx 6 \sqrt{D t_{\max}}\)).

Parameters:

Name Type Description Default
D float

Diffusion coefficient \(D > 0\) (\(\text{m}^2/\text{s}\)).

 required
t_max float

Maximum characteristic simulation time \(t_{\max} > 0\) (\(\text{s}\)).

 required
x_min float

Electrode surface coordinate \(x_{\min}\).

 0.0
factor float

Dimensionless multiplier for diffusion layer thickness \(\delta = \text{factor} \cdot \sqrt{D t_{\max}}\).

 6.0
n_points int

Total number of grid nodes \(N\).

 100
beta float

Dimensionless stretching parameter \(\beta > 0\).

 3.0
two_sided bool

Whether to cluster nodes symmetrically at both boundaries.

 False
mapping str

Coordinate transformation function: "exponential", "tanh", or "sinh".

 "exponential"

Returns:

Type Description
ExponentialMesh

Configured exponential mesh spanning \([x_{\min}, x_{\min} + \delta]\).
Source code in src/softpotato/grid/exponential.py 
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
@classmethod
def from_diffusion_layer(
    cls,
    D: float,
    t_max: float,
    *,
    x_min: float = 0.0,
    factor: float = 6.0,
    n_points: int = 100,
    beta: float = 3.0,
    two_sided: bool = False,
    mapping: str = "exponential",
) -> ExponentialMesh:
    r"""Create an exponential mesh sized to the diffusion layer thickness.

    Sets :math:`x_{\max} = x_{\min} + \text{factor} \cdot \sqrt{D \, t_{\max}}`,
    where ``factor`` defaults to :math:`6.0` (:math:`\delta \approx 6 \sqrt{D t_{\max}}`).

    Parameters
    ----------
    D : float
        Diffusion coefficient :math:`D > 0` (:math:`\text{m}^2/\text{s}`).
    t_max : float
        Maximum characteristic simulation time :math:`t_{\max} > 0` (:math:`\text{s}`).
    x_min : float, default 0.0
        Electrode surface coordinate :math:`x_{\min}`.
    factor : float, default 6.0
        Dimensionless multiplier for diffusion layer thickness
        :math:`\delta = \text{factor} \cdot \sqrt{D t_{\max}}`.
    n_points : int, default 100
        Total number of grid nodes :math:`N`.
    beta : float, default 3.0
        Dimensionless stretching parameter :math:`\beta > 0`.
    two_sided : bool, default False
        Whether to cluster nodes symmetrically at both boundaries.
    mapping : str, default "exponential"
        Coordinate transformation function: ``"exponential"``, ``"tanh"``, or ``"sinh"``.

    Returns
    -------
    ExponentialMesh
        Configured exponential mesh spanning :math:`[x_{\min}, x_{\min} + \delta]`.
    """
    delta = diffusion_layer_thickness(D, t_max, factor=factor)
    if isinstance(x_min, bool):
        raise TypeError("x_min must be a real number.")
    try:
        x_min_f = float(x_min)
    except (TypeError, ValueError) as exc:
        raise TypeError("x_min must be a real number.") from exc

    return cls(
        x_min=x_min_f,
        x_max=x_min_f + delta,
        n_points=n_points,
        beta=beta,
        two_sided=two_sided,
        mapping=mapping,
    )

GeometricExpandingMesh

Bases: Mesh1D

Geometrically expanding 1D spatial discretization mesh.

Generates a non-uniform spatial grid where consecutive interval steps scale by a constant expansion ratio \(\gamma > 1\):

\[ \Delta x_j = \Delta x_0 \cdot \gamma^j \quad (j = 0, \dots, N-2) \]
\[ x_{j+1} = x_j + \Delta x_j \]

Supports three construction parameterizations (specify any two of n_points, gamma, and dx_min):

  1. (x_min, x_max, n_points, gamma): Computes initial step \(\Delta x_0\) so that the geometric series sums to \(x_{\max} - x_{\min}\).
  2. (x_min, x_max, dx_min, gamma): Automatically determines the required number of nodes \(N\) so the initial step is at most dx_min.
  3. (x_min, x_max, n_points, dx_min): Numerically solves for the expansion ratio \(\gamma > 1\) that satisfies the domain length and dx_min.

Also supports two-sided symmetric expansion (two_sided=True), providing fine spatial resolution \(\Delta x_0\) at both boundaries (\(x_{\min}\) and \(x_{\max}\)) while coarsening toward the domain midpoint.

Parameters:

Name Type Description Default
x_min float

Lower boundary of the spatial domain (\(x_{\min}\)).

 required
x_max float

Upper boundary of the spatial domain (\(x_{\max}\)). Must be strictly greater than x_min.

 required
n_points int | None

Total number of grid nodes \(N\).

 None
gamma float | None

Geometric expansion ratio \(\gamma > 1\).

 None
dx_min float | None

Minimum (boundary) interval spacing \(\Delta x_0 > 0\).

 None
two_sided bool

If True, generates a symmetric mesh with fine resolution at both \(x_{\min}\) and \(x_{\max}\), expanding toward the center.

 False

Attributes:

Name Type Description
nodes ndarray

1D array of spatial node positions of length \(N\).
dx ndarray

Interval spacings between consecutive nodes of length \(N-1\).
gamma float

Geometric expansion ratio \(\gamma > 1\).
expansion_ratio float

Alias for gamma.
dx_min float

Initial (minimum) boundary step size \(\Delta x_0\).
dx_0 float

Alias for dx_min.
dx_max float

Maximum interval spacing across the mesh.
two_sided bool

Whether the mesh expands symmetrically from both boundaries.
n_points int

Total number of spatial nodes \(N\).
bounds tuple[float, float]

Spatial domain boundaries \((x_{\min}, x_{\max})\).
is_uniform bool

Always False for geometrically expanding meshes.

Raises:

Type Description
ValueError

If x_max <= x_min or bounds are non-finite. If not exactly two of n_points, gamma, and dx_min are provided. If gamma <= 1.0, dx_min <= 0.0, or n_points is below the minimum required node count. If dx_min * (n_points - 1) >= (x_max - x_min) when solving for gamma.
TypeError

If parameters have invalid types.

Examples:

Constructing from n_points and gamma:

>>> from softpotato.grid import GeometricExpandingMesh
>>> mesh = GeometricExpandingMesh(0.0, 1.0, n_points=21, gamma=1.1)
>>> mesh.n_points
21
>>> mesh.bounds
(0.0, 1.0)
>>> mesh.is_uniform
False
Source code in src/softpotato/grid/geometric.py 
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
class GeometricExpandingMesh(Mesh1D):
    r"""Geometrically expanding 1D spatial discretization mesh.

    Generates a non-uniform spatial grid where consecutive interval steps
    scale by a constant expansion ratio :math:`\gamma > 1`:

    .. math::

        \Delta x_j = \Delta x_0 \cdot \gamma^j \quad (j = 0, \dots, N-2)

    .. math::

        x_{j+1} = x_j + \Delta x_j

    Supports three construction parameterizations (specify any two of
    ``n_points``, ``gamma``, and ``dx_min``):

    1. ``(x_min, x_max, n_points, gamma)``: Computes initial step :math:`\Delta x_0`
       so that the geometric series sums to :math:`x_{\max} - x_{\min}`.
    2. ``(x_min, x_max, dx_min, gamma)``: Automatically determines the required
       number of nodes :math:`N` so the initial step is at most ``dx_min``.
    3. ``(x_min, x_max, n_points, dx_min)``: Numerically solves for the expansion
       ratio :math:`\gamma > 1` that satisfies the domain length and ``dx_min``.

    Also supports two-sided symmetric expansion (``two_sided=True``), providing
    fine spatial resolution :math:`\Delta x_0` at both boundaries (:math:`x_{\min}`
    and :math:`x_{\max}`) while coarsening toward the domain midpoint.

    Parameters
    ----------
    x_min : float
        Lower boundary of the spatial domain (:math:`x_{\min}`).
    x_max : float
        Upper boundary of the spatial domain (:math:`x_{\max}`).
        Must be strictly greater than ``x_min``.
    n_points : int | None, optional
        Total number of grid nodes :math:`N`.
    gamma : float | None, optional
        Geometric expansion ratio :math:`\gamma > 1`.
    dx_min : float | None, optional
        Minimum (boundary) interval spacing :math:`\Delta x_0 > 0`.
    two_sided : bool, default False
        If ``True``, generates a symmetric mesh with fine resolution at both
        :math:`x_{\min}` and :math:`x_{\max}`, expanding toward the center.

    Attributes
    ----------
    nodes : np.ndarray
        1D array of spatial node positions of length :math:`N`.
    dx : np.ndarray
        Interval spacings between consecutive nodes of length :math:`N-1`.
    gamma : float
        Geometric expansion ratio :math:`\gamma > 1`.
    expansion_ratio : float
        Alias for ``gamma``.
    dx_min : float
        Initial (minimum) boundary step size :math:`\Delta x_0`.
    dx_0 : float
        Alias for ``dx_min``.
    dx_max : float
        Maximum interval spacing across the mesh.
    two_sided : bool
        Whether the mesh expands symmetrically from both boundaries.
    n_points : int
        Total number of spatial nodes :math:`N`.
    bounds : tuple[float, float]
        Spatial domain boundaries :math:`(x_{\min}, x_{\max})`.
    is_uniform : bool
        Always ``False`` for geometrically expanding meshes.

    Raises
    ------
    ValueError
        If ``x_max <= x_min`` or bounds are non-finite.
        If not exactly two of ``n_points``, ``gamma``, and ``dx_min`` are provided.
        If ``gamma <= 1.0``, ``dx_min <= 0.0``, or ``n_points`` is below the minimum
        required node count.
        If ``dx_min * (n_points - 1) >= (x_max - x_min)`` when solving for ``gamma``.
    TypeError
        If parameters have invalid types.

    Examples
    --------
    Constructing from ``n_points`` and ``gamma``:

    >>> from softpotato.grid import GeometricExpandingMesh
    >>> mesh = GeometricExpandingMesh(0.0, 1.0, n_points=21, gamma=1.1)
    >>> mesh.n_points
    21
    >>> mesh.bounds
    (0.0, 1.0)
    >>> mesh.is_uniform
    False
    """

    def __init__(
        self,
        x_min: float,
        x_max: float,
        n_points: int | None = None,
        gamma: float | None = None,
        dx_min: float | None = None,
        two_sided: bool = False,
    ) -> None:
        if isinstance(x_min, (bool, str, bytes)) or isinstance(
            x_max, (bool, str, bytes)
        ):
            raise TypeError("x_min and x_max must be real numbers.")
        try:
            x_min_f = float(x_min)
            x_max_f = float(x_max)
        except (TypeError, ValueError) as exc:
            raise TypeError("x_min and x_max must be real numbers.") from exc

        if not (math.isfinite(x_min_f) and math.isfinite(x_max_f)):
            raise ValueError(f"x_min and x_max must be finite, got ({x_min}, {x_max}).")

        if x_max_f <= x_min_f:
            raise ValueError(
                f"x_max ({x_max_f}) must be strictly greater than x_min ({x_min_f})."
            )

        if not isinstance(two_sided, bool):
            raise TypeError(
                f"two_sided must be a boolean, got {type(two_sided).__name__}."
            )

        provided_count = sum(p is not None for p in (n_points, gamma, dx_min))
        if provided_count != 2:
            raise ValueError(
                "Must specify exactly two of 'n_points', 'gamma', and 'dx_min'."
            )

        domain_len = x_max_f - x_min_f

        # Validate n_points if provided
        n_pts: int | None = None
        if n_points is not None:
            if isinstance(n_points, bool) or not isinstance(
                n_points, (int, np.integer)
            ):
                raise TypeError(
                    f"n_points must be an integer, got {type(n_points).__name__}."
                )
            n_pts = int(n_points)
            if two_sided:
                min_pts = 4 if dx_min is not None else 3
                if n_pts < min_pts:
                    raise ValueError(
                        f"n_points must be >= {min_pts} for a two-sided expanding mesh, "
                        f"got {n_pts}."
                    )
            else:
                min_pts = 3 if dx_min is not None else 2
                if n_pts < min_pts:
                    raise ValueError(f"n_points must be >= {min_pts}, got {n_pts}.")

        # Validate gamma if provided
        gamma_f: float | None = None
        if gamma is not None:
            if isinstance(gamma, (bool, str, bytes)):
                raise TypeError("gamma must be a real number.")
            try:
                gamma_f = float(gamma)
            except (TypeError, ValueError) as exc:
                raise TypeError("gamma must be a real number.") from exc
            if not math.isfinite(gamma_f):
                raise ValueError(f"gamma must be finite, got {gamma}.")
            if gamma_f <= 1.0:
                raise ValueError(
                    f"gamma must be strictly greater than 1.0, got {gamma_f}."
                )

        # Validate dx_min if provided
        dx_min_f: float | None = None
        if dx_min is not None:
            if isinstance(dx_min, (bool, str, bytes)):
                raise TypeError("dx_min must be a real number.")
            try:
                dx_min_f = float(dx_min)
            except (TypeError, ValueError) as exc:
                raise TypeError("dx_min must be a real number.") from exc
            if not math.isfinite(dx_min_f):
                raise ValueError(f"dx_min must be finite, got {dx_min}.")
            if dx_min_f <= 0.0:
                raise ValueError(f"dx_min must be strictly positive, got {dx_min_f}.")

        # Case 1: (n_points, gamma) -> compute dx_0
        if n_pts is not None and gamma_f is not None:
            resolved_n = n_pts
            resolved_gamma = gamma_f
            exponents = self._interval_exponents(resolved_n - 1, two_sided)
            weights = resolved_gamma**exponents
            resolved_dx0 = float(domain_len / np.sum(weights))

        # Case 2: (dx_min, gamma) -> determine required n_points
        elif dx_min_f is not None and gamma_f is not None:
            resolved_gamma = gamma_f
            if two_sided:
                half_len = 0.5 * domain_len
                if dx_min_f >= half_len:
                    raise ValueError(
                        f"dx_min ({dx_min_f}) must be strictly less than half the "
                        f"domain length ({half_len}) for a two-sided mesh."
                    )
                k_exact = math.log(
                    1.0 + half_len * (resolved_gamma - 1.0) / dx_min_f
                ) / math.log(resolved_gamma)
                k_half = max(2, math.ceil(k_exact - 1e-12))
                resolved_n = 2 * k_half + 1
            else:
                if dx_min_f >= domain_len:
                    raise ValueError(
                        f"dx_min ({dx_min_f}) must be strictly less than the "
                        f"domain length ({domain_len})."
                    )
                m_exact = math.log(
                    1.0 + domain_len * (resolved_gamma - 1.0) / dx_min_f
                ) / math.log(resolved_gamma)
                m_intervals = max(2, math.ceil(m_exact - 1e-12))
                resolved_n = m_intervals + 1

            exponents = self._interval_exponents(resolved_n - 1, two_sided)
            weights = resolved_gamma**exponents
            resolved_dx0 = float(domain_len / np.sum(weights))

        # Case 3: (n_points, dx_min) -> solve for gamma > 1
        else:
            assert n_pts is not None and dx_min_f is not None
            resolved_n = n_pts
            resolved_dx0 = dx_min_f
            m_intervals = resolved_n - 1

            if resolved_dx0 * m_intervals >= domain_len:
                raise ValueError(
                    f"dx_min * (n_points - 1) ({resolved_dx0 * m_intervals}) must be "
                    f"strictly less than the domain length ({domain_len}) for an "
                    f"expanding mesh (gamma > 1)."
                )

            resolved_gamma = self._solve_gamma(
                domain_len=domain_len,
                dx_min=resolved_dx0,
                m_intervals=m_intervals,
                two_sided=two_sided,
            )
            exponents = self._interval_exponents(m_intervals, two_sided)
            weights = resolved_gamma**exponents
            # Re-normalize dx0 for exact boundary landing
            resolved_dx0 = float(domain_len / np.sum(weights))

        dx_arr = resolved_dx0 * weights
        nodes = np.empty(resolved_n, dtype=float)
        nodes[0] = x_min_f
        nodes[1:] = x_min_f + np.cumsum(dx_arr)
        nodes[-1] = x_max_f
        if two_sided and resolved_n % 2 == 1:
            nodes[resolved_n // 2] = 0.5 * (x_min_f + x_max_f)

        dx_final = np.diff(nodes)
        nodes.flags.writeable = False
        dx_final.flags.writeable = False

        self._nodes = nodes
        self._dx = dx_final
        self._gamma = float(resolved_gamma)
        self._dx_min = float(dx_final[0])
        self._two_sided = two_sided

    @staticmethod
    def _interval_exponents(m_intervals: int, two_sided: bool) -> np.ndarray:
        """Compute geometric power exponents for each of the M intervals."""
        idx = np.arange(m_intervals, dtype=float)
        if two_sided:
            return np.minimum(idx, (m_intervals - 1) - idx)
        return idx

    @classmethod
    def _solve_gamma(
        cls,
        domain_len: float,
        dx_min: float,
        m_intervals: int,
        two_sided: bool,
    ) -> float:
        """Solve for expansion ratio gamma > 1 given domain_len, dx_min, and M."""
        target_ratio = domain_len / dx_min  # > m_intervals
        exponents = cls._interval_exponents(m_intervals, two_sided)
        max_exp = float(np.max(exponents))

        # Work in log space to avoid overflow for large interval counts:
        # log(sum(gamma^exponents)) == log(target_ratio)
        log_target = math.log(target_ratio)

        def objective(g: float) -> float:
            log_g = math.log(g)
            log_terms = exponents * log_g
            max_log = float(np.max(log_terms))
            log_sum = max_log + math.log(float(np.sum(np.exp(log_terms - max_log))))
            return log_sum - log_target

        # Bracket:
        # Since sum(gamma^exponents) < m_intervals * gamma^max_exp,
        # we have gamma > (target_ratio / m_intervals)^(1 / max_exp).
        # And since sum(gamma^exponents) > gamma^max_exp,
        # we have gamma < target_ratio^(1 / max_exp).
        g_low = max(1.0 + 1e-15, (target_ratio / m_intervals) ** (1.0 / max_exp))
        g_high = max(g_low * 1.01, target_ratio ** (1.0 / max_exp))

        while objective(g_low) > 0.0 and g_low > 1.0 + 1e-15:
            g_low = 1.0 + 0.5 * (g_low - 1.0)
        while objective(g_high) < 0.0:
            g_high *= 2.0

        return float(brentq(objective, g_low, g_high, xtol=1e-15, rtol=1e-15))

    @classmethod
    def from_diffusion_layer(
        cls,
        D: float,
        t_max: float,
        *,
        x_min: float = 0.0,
        factor: float = 6.0,
        n_points: int | None = None,
        gamma: float | None = None,
        dx_min: float | None = None,
        two_sided: bool = False,
    ) -> GeometricExpandingMesh:
        r"""Create a geometric expanding mesh sized to the diffusion layer thickness.

        Sets :math:`x_{\max} = x_{\min} + \text{factor} \cdot \sqrt{D \, t_{\max}}`,
        where ``factor`` defaults to :math:`6.0` (:math:`\delta \approx 6 \sqrt{D t_{\max}}`).

        Parameters
        ----------
        D : float
            Diffusion coefficient :math:`D > 0` (:math:`\text{m}^2/\text{s}`).
        t_max : float
            Maximum characteristic simulation time :math:`t_{\max} > 0` (:math:`\text{s}`).
        x_min : float, default 0.0
            Electrode surface coordinate :math:`x_{\min}`.
        factor : float, default 6.0
            Dimensionless multiplier for diffusion layer thickness
            :math:`\delta = \text{factor} \cdot \sqrt{D t_{\max}}`.
        n_points : int | None, optional
            Total number of grid nodes :math:`N`. Defaults to ``100`` if fewer than
            two of ``(n_points, gamma, dx_min)`` are specified.
        gamma : float | None, optional
            Geometric expansion ratio :math:`\gamma > 1`. Defaults to ``1.05`` if
            fewer than two of ``(n_points, gamma, dx_min)`` are specified.
        dx_min : float | None, optional
            Initial boundary step size :math:`\Delta x_0 > 0`.
        two_sided : bool, default False
            Whether to expand symmetrically from both boundaries.

        Returns
        -------
        GeometricExpandingMesh
            Configured geometric mesh spanning :math:`[x_{\min}, x_{\min} + \delta]`.
        """
        delta = diffusion_layer_thickness(D, t_max, factor=factor)
        if isinstance(x_min, bool):
            raise TypeError("x_min must be a real number.")
        try:
            x_min_f = float(x_min)
        except (TypeError, ValueError) as exc:
            raise TypeError("x_min must be a real number.") from exc

        provided = sum(p is not None for p in (n_points, gamma, dx_min))
        if provided == 0:
            n_points = 100
            gamma = 1.05
        elif provided == 1:
            if gamma is None:
                gamma = 1.05
            else:
                n_points = 100

        return cls(
            x_min=x_min_f,
            x_max=x_min_f + delta,
            n_points=n_points,
            gamma=gamma,
            dx_min=dx_min,
            two_sided=two_sided,
        )

    @property
    def nodes(self) -> np.ndarray:
        """1D array of spatial node positions."""
        return self._nodes

    @property
    def dx(self) -> np.ndarray:
        """Interval spacings between consecutive nodes."""
        return self._dx

    @property
    def gamma(self) -> float:
        """Geometric expansion ratio gamma > 1."""
        return self._gamma

    @property
    def expansion_ratio(self) -> float:
        """Alias for geometric expansion ratio gamma."""
        return self._gamma

    @property
    def dx_min(self) -> float:
        """Minimum (boundary) interval spacing Delta x_0."""
        return self._dx_min

    @property
    def dx_0(self) -> float:
        """Alias for initial boundary step size Delta x_0."""
        return self._dx_min

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

    @property
    def two_sided(self) -> bool:
        """Whether the mesh expands symmetrically from both boundaries."""
        return self._two_sided

    @property
    def is_uniform(self) -> bool:
        """Whether the grid spacing is uniform within numerical tolerance."""
        return False

    def __repr__(self) -> str:
        return (
            f"GeometricExpandingMesh(x_min={self.bounds[0]}, x_max={self.bounds[1]}, "
            f"n_points={self.n_points}, gamma={self.gamma:.4g}, "
            f"dx_min={self.dx_min:.4g}, two_sided={self.two_sided})"
        )

    def __str__(self) -> str:
        return (
            f"GeometricExpandingMesh(bounds=({self.bounds[0]}, {self.bounds[1]}), "
            f"n_points={self.n_points}, gamma={self.gamma:.4g}, "
            f"dx_min={self.dx_min:.4g}, two_sided={self.two_sided})"
        )

nodes property 

nodes: ndarray

1D array of spatial node positions.

dx property 

dx: ndarray

Interval spacings between consecutive nodes.

gamma property 

gamma: float

Geometric expansion ratio gamma > 1.

expansion_ratio property 

expansion_ratio: float

Alias for geometric expansion ratio gamma.

dx_min property 

dx_min: float

Minimum (boundary) interval spacing Delta x_0.

dx_0 property 

dx_0: float

Alias for initial boundary step size Delta x_0.

dx_max property 

dx_max: float

Maximum interval spacing across the mesh.

two_sided property 

two_sided: bool

Whether the mesh expands symmetrically from both boundaries.

is_uniform property 

is_uniform: bool

Whether the grid spacing is uniform within numerical tolerance.

from_diffusion_layer classmethod 

from_diffusion_layer(D: float, t_max: float, *, x_min: float = 0.0, factor: float = 6.0, n_points: int | None = None, gamma: float | None = None, dx_min: float | None = None, two_sided: bool = False) -> GeometricExpandingMesh

Create a geometric expanding mesh sized to the diffusion layer thickness.

Sets \(x_{\max} = x_{\min} + \text{factor} \cdot \sqrt{D \, t_{\max}}\), where factor defaults to \(6.0\) (\(\delta \approx 6 \sqrt{D t_{\max}}\)).

Parameters:

Name Type Description Default
D float

Diffusion coefficient \(D > 0\) (\(\text{m}^2/\text{s}\)).

 required
t_max float

Maximum characteristic simulation time \(t_{\max} > 0\) (\(\text{s}\)).

 required
x_min float

Electrode surface coordinate \(x_{\min}\).

 0.0
factor float

Dimensionless multiplier for diffusion layer thickness \(\delta = \text{factor} \cdot \sqrt{D t_{\max}}\).

 6.0
n_points int | None

Total number of grid nodes \(N\). Defaults to 100 if fewer than two of (n_points, gamma, dx_min) are specified.

 None
gamma float | None

Geometric expansion ratio \(\gamma > 1\). Defaults to 1.05 if fewer than two of (n_points, gamma, dx_min) are specified.

 None
dx_min float | None

Initial boundary step size \(\Delta x_0 > 0\).

 None
two_sided bool

Whether to expand symmetrically from both boundaries.

 False

Returns:

Type Description
GeometricExpandingMesh

Configured geometric mesh spanning \([x_{\min}, x_{\min} + \delta]\).
Source code in src/softpotato/grid/geometric.py 
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
@classmethod
def from_diffusion_layer(
    cls,
    D: float,
    t_max: float,
    *,
    x_min: float = 0.0,
    factor: float = 6.0,
    n_points: int | None = None,
    gamma: float | None = None,
    dx_min: float | None = None,
    two_sided: bool = False,
) -> GeometricExpandingMesh:
    r"""Create a geometric expanding mesh sized to the diffusion layer thickness.

    Sets :math:`x_{\max} = x_{\min} + \text{factor} \cdot \sqrt{D \, t_{\max}}`,
    where ``factor`` defaults to :math:`6.0` (:math:`\delta \approx 6 \sqrt{D t_{\max}}`).

    Parameters
    ----------
    D : float
        Diffusion coefficient :math:`D > 0` (:math:`\text{m}^2/\text{s}`).
    t_max : float
        Maximum characteristic simulation time :math:`t_{\max} > 0` (:math:`\text{s}`).
    x_min : float, default 0.0
        Electrode surface coordinate :math:`x_{\min}`.
    factor : float, default 6.0
        Dimensionless multiplier for diffusion layer thickness
        :math:`\delta = \text{factor} \cdot \sqrt{D t_{\max}}`.
    n_points : int | None, optional
        Total number of grid nodes :math:`N`. Defaults to ``100`` if fewer than
        two of ``(n_points, gamma, dx_min)`` are specified.
    gamma : float | None, optional
        Geometric expansion ratio :math:`\gamma > 1`. Defaults to ``1.05`` if
        fewer than two of ``(n_points, gamma, dx_min)`` are specified.
    dx_min : float | None, optional
        Initial boundary step size :math:`\Delta x_0 > 0`.
    two_sided : bool, default False
        Whether to expand symmetrically from both boundaries.

    Returns
    -------
    GeometricExpandingMesh
        Configured geometric mesh spanning :math:`[x_{\min}, x_{\min} + \delta]`.
    """
    delta = diffusion_layer_thickness(D, t_max, factor=factor)
    if isinstance(x_min, bool):
        raise TypeError("x_min must be a real number.")
    try:
        x_min_f = float(x_min)
    except (TypeError, ValueError) as exc:
        raise TypeError("x_min must be a real number.") from exc

    provided = sum(p is not None for p in (n_points, gamma, dx_min))
    if provided == 0:
        n_points = 100
        gamma = 1.05
    elif provided == 1:
        if gamma is None:
            gamma = 1.05
        else:
            n_points = 100

    return cls(
        x_min=x_min_f,
        x_max=x_min_f + delta,
        n_points=n_points,
        gamma=gamma,
        dx_min=dx_min,
        two_sided=two_sided,
    )

Uniform1DMesh

Bases: Mesh1D

Uniform 1D spatial discretization mesh with equidistant nodes.

Generates spatial coordinate nodes \(x_j\) with constant spacing:

\[ x_j = x_{\min} + j \cdot \Delta x, \quad j = 0, 1, \dots, N-1 \]

where \(\Delta x = (x_{\max} - x_{\min}) / (N - 1)\).

Can be constructed either by specifying the total number of points \(N\) (n_points) or the desired uniform spacing \(\Delta x\) (dx).

Parameters:

Name Type Description Default
x_min float

Lower boundary of the spatial domain (\(x_{\min}\)).

 required
x_max float

Upper boundary of the spatial domain (\(x_{\max}\)). Must be strictly greater than x_min.

 required
n_points int | None

Total number of grid nodes \(N\). Must be \(\ge 2\). Exactly one of n_points or dx must be specified.

 None
dx float | None

Target uniform node spacing \(\Delta x\). Must be strictly positive and evenly divide the domain length (\(x_{\max} - x_{\min}\)). Exactly one of n_points or dx must be specified.

 None

Attributes:

Name Type Description
nodes ndarray

1D array of spatial node positions of length \(N\).
dx ndarray

Array of consecutive node intervals of length \(N-1\), where every element equals \(\Delta x\).
spacing float

Constant scalar grid spacing \(\Delta x\).
step float

Alias for spacing.
n_points int

Total number of spatial nodes \(N\).
bounds tuple[float, float]

Spatial domain boundaries \((x_{\min}, x_{\max})\).
is_uniform bool

Always True for uniform meshes.

Raises:

Type Description
ValueError

If x_max <= x_min. If neither n_points nor dx is specified, or if both are specified. If n_points < 2. If dx <= 0 or dx > (x_max - x_min). If domain length is not an integer multiple of dx.
TypeError

If n_points is not an integer, or if inputs cannot be converted to float.

Examples:

Constructing with point count:

>>> from softpotato.grid import Uniform1DMesh
>>> mesh = Uniform1DMesh(0.0, 1.0, n_points=5)
>>> mesh.nodes
array([0.  , 0.25, 0.5 , 0.75, 1.  ])
>>> mesh.spacing
0.25
>>> mesh.is_uniform
True

Constructing with spacing step:

>>> mesh = Uniform1DMesh(0.0, 1.0, dx=0.2)
>>> mesh.n_points
6
>>> mesh.bounds
(0.0, 1.0)
Source code in src/softpotato/grid/uniform.py 
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
class Uniform1DMesh(Mesh1D):
    r"""Uniform 1D spatial discretization mesh with equidistant nodes.

    Generates spatial coordinate nodes :math:`x_j` with constant spacing:

    .. math::

        x_j = x_{\min} + j \cdot \Delta x, \quad j = 0, 1, \dots, N-1

    where :math:`\Delta x = (x_{\max} - x_{\min}) / (N - 1)`.

    Can be constructed either by specifying the total number of points :math:`N`
    (``n_points``) or the desired uniform spacing :math:`\Delta x` (``dx``).

    Parameters
    ----------
    x_min : float
        Lower boundary of the spatial domain (:math:`x_{\min}`).
    x_max : float
        Upper boundary of the spatial domain (:math:`x_{\max}`).
        Must be strictly greater than ``x_min``.
    n_points : int | None, optional
        Total number of grid nodes :math:`N`. Must be :math:`\ge 2`.
        Exactly one of ``n_points`` or ``dx`` must be specified.
    dx : float | None, optional
        Target uniform node spacing :math:`\Delta x`. Must be strictly positive
        and evenly divide the domain length (:math:`x_{\max} - x_{\min}`).
        Exactly one of ``n_points`` or ``dx`` must be specified.

    Attributes
    ----------
    nodes : np.ndarray
        1D array of spatial node positions of length :math:`N`.
    dx : np.ndarray
        Array of consecutive node intervals of length :math:`N-1`, where
        every element equals :math:`\Delta x`.
    spacing : float
        Constant scalar grid spacing :math:`\Delta x`.
    step : float
        Alias for ``spacing``.
    n_points : int
        Total number of spatial nodes :math:`N`.
    bounds : tuple[float, float]
        Spatial domain boundaries :math:`(x_{\min}, x_{\max})`.
    is_uniform : bool
        Always ``True`` for uniform meshes.

    Raises
    ------
    ValueError
        If ``x_max <= x_min``.
        If neither ``n_points`` nor ``dx`` is specified, or if both are specified.
        If ``n_points < 2``.
        If ``dx <= 0`` or ``dx > (x_max - x_min)``.
        If domain length is not an integer multiple of ``dx``.
    TypeError
        If ``n_points`` is not an integer, or if inputs cannot be converted to float.

    Examples
    --------
    Constructing with point count:

    >>> from softpotato.grid import Uniform1DMesh
    >>> mesh = Uniform1DMesh(0.0, 1.0, n_points=5)
    >>> mesh.nodes
    array([0.  , 0.25, 0.5 , 0.75, 1.  ])
    >>> mesh.spacing
    0.25
    >>> mesh.is_uniform
    True

    Constructing with spacing step:

    >>> mesh = Uniform1DMesh(0.0, 1.0, dx=0.2)
    >>> mesh.n_points
    6
    >>> mesh.bounds
    (0.0, 1.0)
    """

    def __init__(
        self,
        x_min: float,
        x_max: float,
        n_points: int | None = None,
        dx: float | None = None,
    ) -> None:
        try:
            x_min_f = float(x_min)
            x_max_f = float(x_max)
        except (TypeError, ValueError) as exc:
            raise TypeError("x_min and x_max must be real numbers.") from exc

        if not (math.isfinite(x_min_f) and math.isfinite(x_max_f)):
            raise ValueError(f"x_min and x_max must be finite, got ({x_min}, {x_max}).")

        if x_max_f <= x_min_f:
            raise ValueError(
                f"x_max ({x_max_f}) must be strictly greater than x_min ({x_min_f})."
            )

        if n_points is None and dx is None:
            raise ValueError("Must specify either 'n_points' or 'dx'.")
        if n_points is not None and dx is not None:
            raise ValueError("Cannot specify both 'n_points' and 'dx'.")

        domain_len = x_max_f - x_min_f

        if n_points is not None:
            if isinstance(n_points, bool) or not isinstance(
                n_points, (int, np.integer)
            ):
                raise TypeError(
                    f"n_points must be an integer, got {type(n_points).__name__}."
                )
            if n_points < 2:
                raise ValueError(f"n_points must be >= 2, got {n_points}.")
            n_pts = int(n_points)
            spacing = domain_len / (n_pts - 1)
        else:
            assert dx is not None
            try:
                dx_f = float(dx)
            except (TypeError, ValueError) as exc:
                raise TypeError("dx must be a real number.") from exc

            if not math.isfinite(dx_f):
                raise ValueError(f"dx must be finite, got {dx}.")
            if dx_f <= 0.0:
                raise ValueError(f"dx must be strictly positive, got {dx_f}.")
            if dx_f > domain_len:
                raise ValueError(
                    f"dx ({dx_f}) cannot be greater than the domain length ({domain_len})."
                )

            n_intervals_float = domain_len / dx_f
            n_intervals = round(n_intervals_float)
            if n_intervals < 1 or not np.isclose(
                n_intervals * dx_f, domain_len, rtol=1e-5, atol=1e-8
            ):
                raise ValueError(
                    f"Domain length ({domain_len}) is not an integer multiple of dx ({dx_f})."
                )
            n_pts = int(n_intervals) + 1
            spacing = dx_f

        nodes = np.linspace(x_min_f, x_max_f, n_pts, dtype=float)
        nodes.flags.writeable = False

        self._nodes = nodes
        self._spacing = float(spacing)
        self._dx = np.diff(self._nodes)
        self._dx.flags.writeable = False

    @classmethod
    def from_diffusion_layer(
        cls,
        D: float,
        t_max: float,
        *,
        x_min: float = 0.0,
        factor: float = 6.0,
        n_points: int | None = None,
        dx: float | None = None,
    ) -> Uniform1DMesh:
        r"""Create a uniform mesh sized to the diffusion layer thickness.

        Sets :math:`x_{\max} = x_{\min} + \text{factor} \cdot \sqrt{D \, t_{\max}}`,
        where ``factor`` defaults to :math:`6.0` (:math:`\delta \approx 6 \sqrt{D t_{\max}}`).

        Parameters
        ----------
        D : float
            Diffusion coefficient :math:`D > 0` (:math:`\text{m}^2/\text{s}`).
        t_max : float
            Maximum characteristic simulation time :math:`t_{\max} > 0` (:math:`\text{s}`).
        x_min : float, default 0.0
            Electrode surface coordinate :math:`x_{\min}`.
        factor : float, default 6.0
            Dimensionless multiplier for diffusion layer thickness
            :math:`\delta = \text{factor} \cdot \sqrt{D t_{\max}}`.
        n_points : int | None, optional
            Total number of grid nodes :math:`N`. Defaults to ``100`` if neither
            ``n_points`` nor ``dx`` is specified.
        dx : float | None, optional
            Uniform grid spacing :math:`\Delta x`.

        Returns
        -------
        Uniform1DMesh
            Configured uniform mesh spanning :math:`[x_{\min}, x_{\min} + \delta]`.
        """
        delta = diffusion_layer_thickness(D, t_max, factor=factor)
        if isinstance(x_min, bool):
            raise TypeError("x_min must be a real number.")
        try:
            x_min_f = float(x_min)
        except (TypeError, ValueError) as exc:
            raise TypeError("x_min must be a real number.") from exc

        if n_points is None and dx is None:
            n_points = 100

        return cls(
            x_min=x_min_f,
            x_max=x_min_f + delta,
            n_points=n_points,
            dx=dx,
        )

    @property
    def nodes(self) -> np.ndarray:
        """1D array of spatial node positions."""
        return self._nodes

    @property
    def dx(self) -> np.ndarray:
        """Interval spacings between consecutive nodes."""
        return self._dx

    @property
    def spacing(self) -> float:
        """Constant grid spacing Delta x."""
        return self._spacing

    @property
    def step(self) -> float:
        """Alias for constant grid spacing Delta x."""
        return self._spacing

    @property
    def is_uniform(self) -> bool:
        """Whether the grid spacing is uniform within numerical tolerance."""
        return True

    def __repr__(self) -> str:
        return (
            f"Uniform1DMesh(x_min={self.bounds[0]}, x_max={self.bounds[1]}, "
            f"n_points={self.n_points}, dx={self.spacing:.4g})"
        )

    def __str__(self) -> str:
        return (
            f"Uniform1DMesh(bounds=({self.bounds[0]}, {self.bounds[1]}), "
            f"n_points={self.n_points}, dx={self.spacing:.4g})"
        )

nodes property 

nodes: ndarray

1D array of spatial node positions.

dx property 

dx: ndarray

Interval spacings between consecutive nodes.

spacing property 

spacing: float

Constant grid spacing Delta x.

step property 

step: float

Alias for constant grid spacing Delta x.

is_uniform property 

is_uniform: bool

Whether the grid spacing is uniform within numerical tolerance.

from_diffusion_layer classmethod 

from_diffusion_layer(D: float, t_max: float, *, x_min: float = 0.0, factor: float = 6.0, n_points: int | None = None, dx: float | None = None) -> Uniform1DMesh

Create a uniform mesh sized to the diffusion layer thickness.

Sets \(x_{\max} = x_{\min} + \text{factor} \cdot \sqrt{D \, t_{\max}}\), where factor defaults to \(6.0\) (\(\delta \approx 6 \sqrt{D t_{\max}}\)).

Parameters:

Name Type Description Default
D float

Diffusion coefficient \(D > 0\) (\(\text{m}^2/\text{s}\)).

 required
t_max float

Maximum characteristic simulation time \(t_{\max} > 0\) (\(\text{s}\)).

 required
x_min float

Electrode surface coordinate \(x_{\min}\).

 0.0
factor float

Dimensionless multiplier for diffusion layer thickness \(\delta = \text{factor} \cdot \sqrt{D t_{\max}}\).

 6.0
n_points int | None

Total number of grid nodes \(N\). Defaults to 100 if neither n_points nor dx is specified.

 None
dx float | None

Uniform grid spacing \(\Delta x\).

 None

Returns:

Type Description
Uniform1DMesh

Configured uniform mesh spanning \([x_{\min}, x_{\min} + \delta]\).
Source code in src/softpotato/grid/uniform.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
@classmethod
def from_diffusion_layer(
    cls,
    D: float,
    t_max: float,
    *,
    x_min: float = 0.0,
    factor: float = 6.0,
    n_points: int | None = None,
    dx: float | None = None,
) -> Uniform1DMesh:
    r"""Create a uniform mesh sized to the diffusion layer thickness.

    Sets :math:`x_{\max} = x_{\min} + \text{factor} \cdot \sqrt{D \, t_{\max}}`,
    where ``factor`` defaults to :math:`6.0` (:math:`\delta \approx 6 \sqrt{D t_{\max}}`).

    Parameters
    ----------
    D : float
        Diffusion coefficient :math:`D > 0` (:math:`\text{m}^2/\text{s}`).
    t_max : float
        Maximum characteristic simulation time :math:`t_{\max} > 0` (:math:`\text{s}`).
    x_min : float, default 0.0
        Electrode surface coordinate :math:`x_{\min}`.
    factor : float, default 6.0
        Dimensionless multiplier for diffusion layer thickness
        :math:`\delta = \text{factor} \cdot \sqrt{D t_{\max}}`.
    n_points : int | None, optional
        Total number of grid nodes :math:`N`. Defaults to ``100`` if neither
        ``n_points`` nor ``dx`` is specified.
    dx : float | None, optional
        Uniform grid spacing :math:`\Delta x`.

    Returns
    -------
    Uniform1DMesh
        Configured uniform mesh spanning :math:`[x_{\min}, x_{\min} + \delta]`.
    """
    delta = diffusion_layer_thickness(D, t_max, factor=factor)
    if isinstance(x_min, bool):
        raise TypeError("x_min must be a real number.")
    try:
        x_min_f = float(x_min)
    except (TypeError, ValueError) as exc:
        raise TypeError("x_min must be a real number.") from exc

    if n_points is None and dx is None:
        n_points = 100

    return cls(
        x_min=x_min_f,
        x_max=x_min_f + delta,
        n_points=n_points,
        dx=dx,
    )

diffusion_layer_thickness 

diffusion_layer_thickness(D: float, t_max: float, factor: float = 6.0) -> float

Estimate the semi-infinite diffusion layer thickness \(\delta\).

Computes the characteristic domain extent required to simulate semi-infinite diffusion over duration \(t_{\max}\) without outer boundary interference:

\[ \delta = \text{factor} \cdot \sqrt{D \, t_{\max}} \]

At the default factor=6.0, \(\operatorname{erfc}(6 / 2) \approx 2.2 \times 10^{-5}\), ensuring concentration perturbations at \(x = \delta\) remain negligible.

Parameters:

Name Type Description Default
D float

Species diffusion coefficient \(D > 0\) (\(\text{m}^2/\text{s}\)).

 required
t_max float

Maximum characteristic simulation time \(t_{\max} > 0\) (\(\text{s}\)).

 required
factor float

Dimensionless multiplier \(> 0\) applied to \(\sqrt{D t_{\max}}\).

 6.0

Returns:

Type Description
float

Estimated diffusion layer thickness \(\delta > 0\).

Raises:

Type Description
ValueError

If D <= 0, t_max <= 0, factor <= 0, or any input is non-finite.
TypeError

If inputs are not real numbers.
Source code in src/softpotato/grid/base.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def diffusion_layer_thickness(
    D: float,
    t_max: float,
    factor: float = 6.0,
) -> float:
    r"""Estimate the semi-infinite diffusion layer thickness :math:`\delta`.

    Computes the characteristic domain extent required to simulate semi-infinite
    diffusion over duration :math:`t_{\max}` without outer boundary interference:

    .. math::

        \delta = \text{factor} \cdot \sqrt{D \, t_{\max}}

    At the default ``factor=6.0``, :math:`\operatorname{erfc}(6 / 2) \approx 2.2 \times 10^{-5}`,
    ensuring concentration perturbations at :math:`x = \delta` remain negligible.

    Parameters
    ----------
    D : float
        Species diffusion coefficient :math:`D > 0` (:math:`\text{m}^2/\text{s}`).
    t_max : float
        Maximum characteristic simulation time :math:`t_{\max} > 0` (:math:`\text{s}`).
    factor : float, default 6.0
        Dimensionless multiplier :math:`> 0` applied to :math:`\sqrt{D t_{\max}}`.

    Returns
    -------
    float
        Estimated diffusion layer thickness :math:`\delta > 0`.

    Raises
    ------
    ValueError
        If ``D <= 0``, ``t_max <= 0``, ``factor <= 0``, or any input is non-finite.
    TypeError
        If inputs are not real numbers.
    """
    for name, val in (("D", D), ("t_max", t_max), ("factor", factor)):
        if isinstance(val, (bool, str, bytes)):
            raise TypeError(f"{name} must be a real number.")
        try:
            val_f = float(val)
        except (TypeError, ValueError) as exc:
            raise TypeError(f"{name} must be a real number.") from exc
        if not math.isfinite(val_f):
            raise ValueError(f"{name} must be finite, got {val}.")
        if val_f <= 0.0:
            raise ValueError(f"{name} must be strictly positive, got {val_f}.")

    return float(factor) * math.sqrt(float(D) * float(t_max))

create_diffusion_mesh 

create_diffusion_mesh(D: float, t_max: float, *, kind: str = 'geometric', x_min: float = 0.0, factor: float = 6.0, n_points: int | None = None, gamma: float | None = None, dx_min: float | None = None, beta: float = 3.0, dx: float | None = None, two_sided: bool = False, mapping: str = 'exponential') -> Mesh1D

Create a 1D spatial mesh sized to the diffusion layer thickness \(\delta\).

Automatically sets the outer domain boundary to \(x_{\max} = x_{\min} + \text{factor} \cdot \sqrt{D \, t_{\max}}\) (with default \(\delta = 6 \sqrt{D t_{\max}}\)) and constructs a GeometricExpandingMesh, ExponentialMesh, or Uniform1DMesh.

Parameters:

Name Type Description Default
D float

Species diffusion coefficient \(D > 0\) (\(\text{m}^2/\text{s}\)).

 required
t_max float

Maximum characteristic simulation time \(t_{\max} > 0\) (\(\text{s}\)).

 required
kind str

Mesh generator type: "geometric", "exponential", or "uniform".

 "geometric"
x_min float

Lower boundary coordinate \(x_{\min}\).

 0.0
factor float

Dimensionless multiplier for diffusion layer thickness \(\delta = \text{factor} \cdot \sqrt{D t_{\max}}\).

 6.0
n_points int | None

Total number of spatial nodes \(N\).

 None
gamma float | None

Expansion ratio \(\gamma > 1\) (for kind="geometric").

 None
dx_min float | None

Minimum boundary step size \(\Delta x_0 > 0\) (for kind="geometric").

 None
beta float

Stretching parameter \(\beta > 0\) (for kind="exponential").

 3.0
dx float | None

Uniform step size \(\Delta x > 0\) (for kind="uniform").

 None
two_sided bool

Whether to cluster nodes symmetrically at both boundaries (for kind="geometric" or kind="exponential").

 False
mapping str

Coordinate transformation ("exponential", "tanh", or "sinh") when kind="exponential".

 "exponential"

Returns:

Type Description
Mesh1D

Constructed 1D mesh instance spanning \([x_{\min}, x_{\min} + \delta]\).

Raises:

Type Description
ValueError

If kind is unknown or mesh parameters are invalid.
TypeError

If kind is not a string or parameters have invalid types.
Source code in src/softpotato/grid/__init__.py 
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def create_diffusion_mesh(
    D: float,
    t_max: float,
    *,
    kind: str = "geometric",
    x_min: float = 0.0,
    factor: float = 6.0,
    n_points: int | None = None,
    gamma: float | None = None,
    dx_min: float | None = None,
    beta: float = 3.0,
    dx: float | None = None,
    two_sided: bool = False,
    mapping: str = "exponential",
) -> Mesh1D:
    r"""Create a 1D spatial mesh sized to the diffusion layer thickness :math:`\delta`.

    Automatically sets the outer domain boundary to
    :math:`x_{\max} = x_{\min} + \text{factor} \cdot \sqrt{D \, t_{\max}}`
    (with default :math:`\delta = 6 \sqrt{D t_{\max}}`) and constructs a
    ``GeometricExpandingMesh``, ``ExponentialMesh``, or ``Uniform1DMesh``.

    Parameters
    ----------
    D : float
        Species diffusion coefficient :math:`D > 0` (:math:`\text{m}^2/\text{s}`).
    t_max : float
        Maximum characteristic simulation time :math:`t_{\max} > 0` (:math:`\text{s}`).
    kind : str, default "geometric"
        Mesh generator type: ``"geometric"``, ``"exponential"``, or ``"uniform"``.
    x_min : float, default 0.0
        Lower boundary coordinate :math:`x_{\min}`.
    factor : float, default 6.0
        Dimensionless multiplier for diffusion layer thickness
        :math:`\delta = \text{factor} \cdot \sqrt{D t_{\max}}`.
    n_points : int | None, optional
        Total number of spatial nodes :math:`N`.
    gamma : float | None, optional
        Expansion ratio :math:`\gamma > 1` (for ``kind="geometric"``).
    dx_min : float | None, optional
        Minimum boundary step size :math:`\Delta x_0 > 0` (for ``kind="geometric"``).
    beta : float, default 3.0
        Stretching parameter :math:`\beta > 0` (for ``kind="exponential"``).
    dx : float | None, optional
        Uniform step size :math:`\Delta x > 0` (for ``kind="uniform"``).
    two_sided : bool, default False
        Whether to cluster nodes symmetrically at both boundaries
        (for ``kind="geometric"`` or ``kind="exponential"``).
    mapping : str, default "exponential"
        Coordinate transformation (``"exponential"``, ``"tanh"``, or ``"sinh"``)
        when ``kind="exponential"``.

    Returns
    -------
    Mesh1D
        Constructed 1D mesh instance spanning :math:`[x_{\min}, x_{\min} + \delta]`.

    Raises
    ------
    ValueError
        If ``kind`` is unknown or mesh parameters are invalid.
    TypeError
        If ``kind`` is not a string or parameters have invalid types.
    """
    if not isinstance(kind, str):
        raise TypeError(f"kind must be a string, got {type(kind).__name__}.")
    kind_norm = kind.strip().lower()

    if kind_norm == "geometric":
        return GeometricExpandingMesh.from_diffusion_layer(
            D=D,
            t_max=t_max,
            x_min=x_min,
            factor=factor,
            n_points=n_points,
            gamma=gamma,
            dx_min=dx_min,
            two_sided=two_sided,
        )
    if kind_norm == "exponential":
        return ExponentialMesh.from_diffusion_layer(
            D=D,
            t_max=t_max,
            x_min=x_min,
            factor=factor,
            n_points=100 if n_points is None else n_points,
            beta=beta,
            two_sided=two_sided,
            mapping=mapping,
        )
    if kind_norm == "uniform":
        return Uniform1DMesh.from_diffusion_layer(
            D=D,
            t_max=t_max,
            x_min=x_min,
            factor=factor,
            n_points=n_points,
            dx=dx,
        )

    raise ValueError(
        f"Unknown mesh kind '{kind}'. Must be one of ('geometric', 'exponential', 'uniform')."
    )