From e5ed12d01192292536f82099ae086c47254b29a2 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Mon, 25 May 2026 22:17:29 -0700 Subject: [PATCH 01/66] Initial implementation by Angel Boada. --- examples/python/simple_operators_check.py | 86 +++++++++++++++++++++++ src/python/pymole/__init__.py | 15 ++++ src/python/pymole/dim1d/__init__.py | 13 ++++ src/python/pymole/dim1d/divergence.py | 36 ++++++++++ src/python/pymole/dim1d/gradient.py | 38 ++++++++++ src/python/pymole/dim1d/grid.py | 58 +++++++++++++++ src/python/pymole/dim1d/laplacian.py | 20 ++++++ src/python/pymole/dim1d/operator.py | 19 +++++ src/python/pymole/dim2d/__init__.py | 11 +++ src/python/pymole/dim2d/divergence.py | 26 +++++++ src/python/pymole/dim2d/gradient.py | 26 +++++++ src/python/pymole/dim2d/laplacian.py | 20 ++++++ src/python/pymole/dim2d/operator.py | 23 ++++++ src/python/pymole/dim3d/__init__.py | 7 ++ src/python/pymole/dim3d/gradient.py | 30 ++++++++ src/python/pymole/dim3d/operator.py | 33 +++++++++ src/python/pymole/operator.py | 56 +++++++++++++++ src/python/pymole/requirements.txt | 2 + src/python/pyproject.toml | 9 +++ src/python/test/__init__.py | 0 src/python/test/dim1d/__init__.py | 0 src/python/test/dim1d/test_divergence.py | 73 +++++++++++++++++++ src/python/test/dim1d/test_gradient.py | 80 +++++++++++++++++++++ src/python/test/dim1d/test_grid.py | 85 ++++++++++++++++++++++ src/python/test/dim1d/test_laplacian.py | 54 ++++++++++++++ src/python/test/dim1d/test_operator.py | 47 +++++++++++++ src/python/test/dim2d/__init__.py | 0 src/python/test/dim2d/test_divergence.py | 57 +++++++++++++++ src/python/test/dim2d/test_gradient.py | 63 +++++++++++++++++ src/python/test/dim2d/test_laplacian.py | 68 ++++++++++++++++++ src/python/test/dim2d/test_operator.py | 59 ++++++++++++++++ src/python/test/dim3d/__init__.py | 0 src/python/test/dim3d/test_gradient.py | 76 ++++++++++++++++++++ src/python/test/dim3d/test_operator.py | 70 ++++++++++++++++++ 34 files changed, 1260 insertions(+) create mode 100644 examples/python/simple_operators_check.py create mode 100644 src/python/pymole/__init__.py create mode 100644 src/python/pymole/dim1d/__init__.py create mode 100644 src/python/pymole/dim1d/divergence.py create mode 100644 src/python/pymole/dim1d/gradient.py create mode 100644 src/python/pymole/dim1d/grid.py create mode 100644 src/python/pymole/dim1d/laplacian.py create mode 100644 src/python/pymole/dim1d/operator.py create mode 100644 src/python/pymole/dim2d/__init__.py create mode 100644 src/python/pymole/dim2d/divergence.py create mode 100644 src/python/pymole/dim2d/gradient.py create mode 100644 src/python/pymole/dim2d/laplacian.py create mode 100644 src/python/pymole/dim2d/operator.py create mode 100644 src/python/pymole/dim3d/__init__.py create mode 100644 src/python/pymole/dim3d/gradient.py create mode 100644 src/python/pymole/dim3d/operator.py create mode 100644 src/python/pymole/operator.py create mode 100644 src/python/pymole/requirements.txt create mode 100644 src/python/pyproject.toml create mode 100644 src/python/test/__init__.py create mode 100644 src/python/test/dim1d/__init__.py create mode 100644 src/python/test/dim1d/test_divergence.py create mode 100644 src/python/test/dim1d/test_gradient.py create mode 100644 src/python/test/dim1d/test_grid.py create mode 100644 src/python/test/dim1d/test_laplacian.py create mode 100644 src/python/test/dim1d/test_operator.py create mode 100644 src/python/test/dim2d/__init__.py create mode 100644 src/python/test/dim2d/test_divergence.py create mode 100644 src/python/test/dim2d/test_gradient.py create mode 100644 src/python/test/dim2d/test_laplacian.py create mode 100644 src/python/test/dim2d/test_operator.py create mode 100644 src/python/test/dim3d/__init__.py create mode 100644 src/python/test/dim3d/test_gradient.py create mode 100644 src/python/test/dim3d/test_operator.py diff --git a/examples/python/simple_operators_check.py b/examples/python/simple_operators_check.py new file mode 100644 index 000000000..eceae4ac8 --- /dev/null +++ b/examples/python/simple_operators_check.py @@ -0,0 +1,86 @@ +import sys +sys.path.insert(0, "src/python") + +import numpy as np +from numpy.typing import NDArray +from scipy import sparse + +from pymole import ( + Divergence1D, + Divergence2D, + Gradient1D, + Gradient2D, + Grid1D, + Laplacian1D, + Laplacian2D, +) + + +def report_grid(name: str, grid: Grid1D) -> None: + print(f"{name}: {grid}; cell size (resolution): {grid.cell_size:.4f}") + + +def report_nullity_check( + name: str, + op: sparse.csr_matrix, + expected_shape: tuple[int, int], + constant_field: NDArray[np.float64], +) -> None: + print(f"{name} shape: {op.shape}; expected {expected_shape}") + print( + f"||{name} @ constant||_inf = {np.linalg.norm(op @ constant_field, ord=np.inf)}" + ) + + +def check_1d_operators(accuracy_order: int = 2) -> None: + num_cells = 2 * accuracy_order + 1 + + grid = Grid1D(0.0, 1.0, num_cells) + + G = Gradient1D(grid, accuracy_order).matrix + D = Divergence1D(grid, accuracy_order).matrix + L = Laplacian1D(grid, accuracy_order).matrix + + f_constant = np.ones(num_cells + 2) + v_constant = np.ones(num_cells + 1) + + print("=== 1D operators ===") + report_grid("x grid", grid) + report_nullity_check("G", G, (num_cells + 1, num_cells + 2), f_constant) + report_nullity_check("D", D, (num_cells + 2, num_cells + 1), v_constant) + report_nullity_check("L", L, (num_cells + 2, num_cells + 2), f_constant) + print(f"D @ G == L? -> {abs(D @ G - L).max() < 1e-12}") + print() + + +def check_2d_operators(accuracy_order: int = 2) -> None: + x_grid = Grid1D(-1.0, 1.0, 2 * accuracy_order + 1) + y_grid = Grid1D(-0.5, 2.0, 2 * accuracy_order + 2) + + m = x_grid.num_cells + n = y_grid.num_cells + + G = Gradient2D(x_grid, y_grid, accuracy_order).matrix + D = Divergence2D(x_grid, y_grid, accuracy_order).matrix + L = Laplacian2D(x_grid, y_grid, accuracy_order).matrix + + f_constant = np.ones((m + 2) * (n + 2)) + v_constant = np.ones(2 * m * n + m + n) + + print("=== 2D operators ===") + report_grid("x grid", x_grid) + report_grid("y grid", y_grid) + report_nullity_check("G", G, (2 * m * n + m + n, (m + 2) * (n + 2)), f_constant) + report_nullity_check("D", D, ((m + 2) * (n + 2), 2 * m * n + m + n), v_constant) + report_nullity_check("L", L, ((m + 2) * (n + 2), (m + 2) * (n + 2)), f_constant) + print(f"D @ G == L? -> {abs(D @ G - L).max() < 1e-12}") + print() + + +def main() -> None: + check_1d_operators() + check_2d_operators() + + +if __name__ == "__main__": + main() diff --git a/src/python/pymole/__init__.py b/src/python/pymole/__init__.py new file mode 100644 index 000000000..a8be89aa0 --- /dev/null +++ b/src/python/pymole/__init__.py @@ -0,0 +1,15 @@ +from .dim1d import Divergence1D, Gradient1D, Grid1D, Laplacian1D +from .dim2d import Divergence2D, Gradient2D, Laplacian2D +from .dim3d import Gradient3D + +__version__ = "0.1.0" +__all__ = [ + "Divergence1D", + "Divergence2D", + "Gradient1D", + "Gradient2D", + "Gradient3D", + "Grid1D", + "Laplacian1D", + "Laplacian2D", +] diff --git a/src/python/pymole/dim1d/__init__.py b/src/python/pymole/dim1d/__init__.py new file mode 100644 index 000000000..39651117d --- /dev/null +++ b/src/python/pymole/dim1d/__init__.py @@ -0,0 +1,13 @@ +from .divergence import Divergence1D +from .gradient import Gradient1D +from .grid import Grid1D +from .laplacian import Laplacian1D +from .operator import MimeticOperator1D + +__all__ = [ + "Divergence1D", + "Gradient1D", + "Grid1D", + "Laplacian1D", + "MimeticOperator1D", +] diff --git a/src/python/pymole/dim1d/divergence.py b/src/python/pymole/dim1d/divergence.py new file mode 100644 index 000000000..45a37a6de --- /dev/null +++ b/src/python/pymole/dim1d/divergence.py @@ -0,0 +1,36 @@ +import numpy as np +from scipy import sparse + +from .operator import MimeticOperator1D + + +class Divergence1D(MimeticOperator1D): + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return 2 * accuracy_order + 1 + + def _construct_matrix(self) -> sparse.csr_matrix: + match self.accuracy_order: + case 2: + return self._divergence_2nd_order() + case _: + raise ValueError( + f"Unsupported {self.__class__.__name__} " + f"order of accuracy: {self.accuracy_order}." + ) + + def _divergence_2nd_order(self) -> sparse.csr_matrix: + m = self.x_grid.num_cells + delta_x = self.x_grid.cell_size + + main_diag = np.ones(m + 1) + sub_diag = np.full(m + 1, -1.0) + + main_diag[0] = 0.0 + sub_diag[-1] = 0.0 + + D = sparse.diags( + [sub_diag, main_diag], offsets=[-1, 0], shape=(m + 2, m + 1) + ).tocsr() + + return D / delta_x diff --git a/src/python/pymole/dim1d/gradient.py b/src/python/pymole/dim1d/gradient.py new file mode 100644 index 000000000..144bdf1cb --- /dev/null +++ b/src/python/pymole/dim1d/gradient.py @@ -0,0 +1,38 @@ +import numpy as np +from scipy import sparse + +from .operator import MimeticOperator1D + + +class Gradient1D(MimeticOperator1D): + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return 2 * accuracy_order + + def _construct_matrix(self) -> sparse.csr_matrix: + match self.accuracy_order: + case 2: + return self._gradient_2nd_order() + case _: + raise ValueError( + f"Unsupported {self.__class__.__name__} " + f"order of accuracy: {self.accuracy_order}." + ) + + def _gradient_2nd_order(self) -> sparse.csr_matrix: + m = self.x_grid.num_cells + delta_x = self.x_grid.cell_size + + main_diag = np.full(m + 1, -1.0) + super_diag = np.ones(m + 1) + + main_diag[-1] = 0.0 + + G = sparse.diags( + [main_diag, super_diag], offsets=[0, 1], shape=(m + 1, m + 2) + ).tolil() + + G[0, 0:3] = [-8.0 / 3.0, 3.0, -1.0 / 3.0] + G[m, m - 1 : m + 2] = [1.0 / 3.0, -3.0, 8.0 / 3.0] + + return G.tocsr() / delta_x diff --git a/src/python/pymole/dim1d/grid.py b/src/python/pymole/dim1d/grid.py new file mode 100644 index 000000000..33b585dbb --- /dev/null +++ b/src/python/pymole/dim1d/grid.py @@ -0,0 +1,58 @@ +from math import isclose + + +class Grid1D: + def __init__(self, lower_bound: float, upper_bound: float, num_cells: int): + if lower_bound >= upper_bound: + raise ValueError( + f"Lower bound ({lower_bound}) value must be less than " + f"the upper bound ({upper_bound})", + ) + if num_cells < 1: + raise ValueError(f"Number of cells ({num_cells}) must be at least 1") + + self._lower_bound = lower_bound + self._upper_bound = upper_bound + self._num_cells = num_cells + + @classmethod + def from_number_of_cells(cls, num_cells: int): + if num_cells < 1: + raise ValueError(f"Number of cells ({num_cells}) must be at least 1") + + lower_bound = 0.0 + upper_bound = 1.0 + + return cls(lower_bound, upper_bound, num_cells) + + @property + def lower_bound(self) -> float: + return self._lower_bound + + @property + def upper_bound(self) -> float: + return self._upper_bound + + @property + def num_cells(self) -> int: + return self._num_cells + + @property + def cell_size(self) -> float: + return (self.upper_bound - self.lower_bound) / self.num_cells + + def __repr__(self) -> str: + return ( + f"Grid1D(lower_bound={self.lower_bound}, " + f"upper_bound={self.upper_bound}, " + f"num_cells={self.num_cells})" + ) + + def __eq__(self, other: object) -> bool: + if not isinstance(other, Grid1D): + return False + return ( + isclose(self.lower_bound, other.lower_bound) + and isclose(self.upper_bound, other.upper_bound) + and self.num_cells == other.num_cells + ) diff --git a/src/python/pymole/dim1d/laplacian.py b/src/python/pymole/dim1d/laplacian.py new file mode 100644 index 000000000..95ec9a163 --- /dev/null +++ b/src/python/pymole/dim1d/laplacian.py @@ -0,0 +1,20 @@ +from scipy import sparse + +from .divergence import Divergence1D +from .gradient import Gradient1D +from .operator import MimeticOperator1D + + +class Laplacian1D(MimeticOperator1D): + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return max( + Divergence1D.minimum_cells_for_order(accuracy_order), + Gradient1D.minimum_cells_for_order(accuracy_order), + ) + + def _construct_matrix(self) -> sparse.csr_matrix: + D = Divergence1D(self.x_grid, self.accuracy_order) + G = Gradient1D(self.x_grid, self.accuracy_order) + + return D.matrix @ G.matrix diff --git a/src/python/pymole/dim1d/operator.py b/src/python/pymole/dim1d/operator.py new file mode 100644 index 000000000..255b832b1 --- /dev/null +++ b/src/python/pymole/dim1d/operator.py @@ -0,0 +1,19 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING + +from pymole.operator import MimeticOperator + +if TYPE_CHECKING: + from .grid import Grid1D + + +class MimeticOperator1D(MimeticOperator): + axis_names = ("x",) + + def __init__(self, x_grid: Grid1D, accuracy_order: int = 2) -> None: + super().__init__((x_grid,), accuracy_order) + + @property + def x_grid(self) -> Grid1D: + return self.grids[0] diff --git a/src/python/pymole/dim2d/__init__.py b/src/python/pymole/dim2d/__init__.py new file mode 100644 index 000000000..909b45cd5 --- /dev/null +++ b/src/python/pymole/dim2d/__init__.py @@ -0,0 +1,11 @@ +from .divergence import Divergence2D +from .gradient import Gradient2D +from .laplacian import Laplacian2D +from .operator import MimeticOperator2D + +__all__ = [ + "Divergence2D", + "Gradient2D", + "Laplacian2D", + "MimeticOperator2D", +] diff --git a/src/python/pymole/dim2d/divergence.py b/src/python/pymole/dim2d/divergence.py new file mode 100644 index 000000000..e3fcfb0ee --- /dev/null +++ b/src/python/pymole/dim2d/divergence.py @@ -0,0 +1,26 @@ +from scipy import sparse + +from pymole.dim1d.divergence import Divergence1D + +from .operator import MimeticOperator2D + + +class Divergence2D(MimeticOperator2D): + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return Divergence1D.minimum_cells_for_order(accuracy_order) + + def _construct_matrix(self) -> sparse.csr_matrix: + m = self.x_grid.num_cells + n = self.y_grid.num_cells + + Dx = Divergence1D(self.x_grid, self.accuracy_order).matrix + Dy = Divergence1D(self.y_grid, self.accuracy_order).matrix + + Im = sparse.eye(m + 2, format="csr")[:, 1:-1] + In = sparse.eye(n + 2, format="csr")[:, 1:-1] + + Sx = sparse.kron(In, Dx, format="csr") + Sy = sparse.kron(Dy, Im, format="csr") + + return sparse.hstack([Sx, Sy], format="csr") diff --git a/src/python/pymole/dim2d/gradient.py b/src/python/pymole/dim2d/gradient.py new file mode 100644 index 000000000..5cb5272fd --- /dev/null +++ b/src/python/pymole/dim2d/gradient.py @@ -0,0 +1,26 @@ +from scipy import sparse + +from pymole.dim1d.gradient import Gradient1D + +from .operator import MimeticOperator2D + + +class Gradient2D(MimeticOperator2D): + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return Gradient1D.minimum_cells_for_order(accuracy_order) + + def _construct_matrix(self) -> sparse.csr_matrix: + m = self.x_grid.num_cells + n = self.y_grid.num_cells + + Gx = Gradient1D(self.x_grid, self.accuracy_order).matrix + Gy = Gradient1D(self.y_grid, self.accuracy_order).matrix + + Im = sparse.eye(m + 2, format="csr")[1:-1, :] + In = sparse.eye(n + 2, format="csr")[1:-1, :] + + Sx = sparse.kron(In, Gx, format="csr") + Sy = sparse.kron(Gy, Im, format="csr") + + return sparse.vstack([Sx, Sy], format="csr") diff --git a/src/python/pymole/dim2d/laplacian.py b/src/python/pymole/dim2d/laplacian.py new file mode 100644 index 000000000..997545c8a --- /dev/null +++ b/src/python/pymole/dim2d/laplacian.py @@ -0,0 +1,20 @@ +from scipy import sparse + +from .divergence import Divergence2D +from .gradient import Gradient2D +from .operator import MimeticOperator2D + + +class Laplacian2D(MimeticOperator2D): + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return max( + Divergence2D.minimum_cells_for_order(accuracy_order), + Gradient2D.minimum_cells_for_order(accuracy_order), + ) + + def _construct_matrix(self) -> sparse.csr_matrix: + D = Divergence2D(self.x_grid, self.y_grid, self.accuracy_order) + G = Gradient2D(self.x_grid, self.y_grid, self.accuracy_order) + + return D.matrix @ G.matrix diff --git a/src/python/pymole/dim2d/operator.py b/src/python/pymole/dim2d/operator.py new file mode 100644 index 000000000..5e847d93a --- /dev/null +++ b/src/python/pymole/dim2d/operator.py @@ -0,0 +1,23 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING + +from pymole.operator import MimeticOperator + +if TYPE_CHECKING: + from pymole.dim1d.grid import Grid1D + + +class MimeticOperator2D(MimeticOperator): + axis_names = ("x", "y") + + def __init__(self, x_grid: Grid1D, y_grid: Grid1D, accuracy_order: int = 2) -> None: + super().__init__((x_grid, y_grid), accuracy_order) + + @property + def x_grid(self) -> Grid1D: + return self.grids[0] + + @property + def y_grid(self) -> Grid1D: + return self.grids[1] diff --git a/src/python/pymole/dim3d/__init__.py b/src/python/pymole/dim3d/__init__.py new file mode 100644 index 000000000..38df6ef1f --- /dev/null +++ b/src/python/pymole/dim3d/__init__.py @@ -0,0 +1,7 @@ +from .gradient import Gradient3D +from .operator import MimeticOperator3D + +__all__ = [ + "Gradient3D", + "MimeticOperator3D", +] diff --git a/src/python/pymole/dim3d/gradient.py b/src/python/pymole/dim3d/gradient.py new file mode 100644 index 000000000..0d4d12c93 --- /dev/null +++ b/src/python/pymole/dim3d/gradient.py @@ -0,0 +1,30 @@ +from scipy import sparse + +from pymole.dim1d.gradient import Gradient1D + +from .operator import MimeticOperator3D + + +class Gradient3D(MimeticOperator3D): + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return Gradient1D.minimum_cells_for_order(accuracy_order) + + def _construct_matrix(self) -> sparse.csr_matrix: + m = self.x_grid.num_cells + n = self.y_grid.num_cells + o = self.z_grid.num_cells + + Gx = Gradient1D(self.x_grid, self.accuracy_order).matrix + Gy = Gradient1D(self.y_grid, self.accuracy_order).matrix + Gz = Gradient1D(self.z_grid, self.accuracy_order).matrix + + Im = sparse.eye(m + 2, format="csr")[1:-1, :] + In = sparse.eye(n + 2, format="csr")[1:-1, :] + Io = sparse.eye(o + 2, format="csr")[1:-1, :] + + Sx = sparse.kron(sparse.kron(Io, In, format="csr"), Gx, format="csr") + Sy = sparse.kron(sparse.kron(Io, Gy, format="csr"), Im, format="csr") + Sz = sparse.kron(sparse.kron(Gz, In, format="csr"), Im, format="csr") + + return sparse.vstack([Sx, Sy, Sz], format="csr") diff --git a/src/python/pymole/dim3d/operator.py b/src/python/pymole/dim3d/operator.py new file mode 100644 index 000000000..43a5963ef --- /dev/null +++ b/src/python/pymole/dim3d/operator.py @@ -0,0 +1,33 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING + +from pymole.operator import MimeticOperator + +if TYPE_CHECKING: + from pymole.dim1d.grid import Grid1D + + +class MimeticOperator3D(MimeticOperator): + axis_names = ("x", "y", "z") + + def __init__( + self, + x_grid: Grid1D, + y_grid: Grid1D, + z_grid: Grid1D, + accuracy_order: int = 2, + ) -> None: + super().__init__((x_grid, y_grid, z_grid), accuracy_order) + + @property + def x_grid(self) -> Grid1D: + return self.grids[0] + + @property + def y_grid(self) -> Grid1D: + return self.grids[1] + + @property + def z_grid(self) -> Grid1D: + return self.grids[2] diff --git a/src/python/pymole/operator.py b/src/python/pymole/operator.py new file mode 100644 index 000000000..465b33e0e --- /dev/null +++ b/src/python/pymole/operator.py @@ -0,0 +1,56 @@ +from __future__ import annotations + +from abc import ABC, abstractmethod +from functools import cached_property +from typing import TYPE_CHECKING, ClassVar + +from scipy import sparse + +if TYPE_CHECKING: + from pymole.dim1d.grid import Grid1D + + +class MimeticOperator(ABC): + axis_names: ClassVar[tuple[str, ...]] + + def __init__(self, grids: tuple[Grid1D, ...], accuracy_order: int = 2) -> None: + if len(grids) != len(self.axis_names): + raise TypeError( + f"{type(self).__name__} expects {len(self.axis_names)} " + f"grid(s) (one per axis in {self.axis_names}), got {len(grids)}." + ) + + self.grids = grids + self.accuracy_order = accuracy_order + + min_cells = self.minimum_cells_for_order(accuracy_order) + undersized = [ + f"{name}: {g.num_cells}" + for name, g in zip(self.axis_names, grids) + if g.num_cells < min_cells + ] + if undersized: + raise ValueError( + f"Grid too small for requested accuracy: " + f"Number of cells in {', '.join(undersized)} " + f"must be at least {min_cells} for a " + f"{accuracy_order}-order accurate {len(grids)}D " + f"{type(self).__name__}.", + ) + + @property + def minimum_cells(self) -> int: + return self.minimum_cells_for_order(self.accuracy_order) + + @classmethod + @abstractmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + pass + + @cached_property + def matrix(self) -> sparse.csr_matrix: + return self._construct_matrix() + + @abstractmethod + def _construct_matrix(self) -> sparse.csr_matrix: + pass diff --git a/src/python/pymole/requirements.txt b/src/python/pymole/requirements.txt new file mode 100644 index 000000000..9d2ad45b9 --- /dev/null +++ b/src/python/pymole/requirements.txt @@ -0,0 +1,2 @@ +numpy==2.4.0 +scipy==1.16.3 diff --git a/src/python/pyproject.toml b/src/python/pyproject.toml new file mode 100644 index 000000000..ebc9e7df3 --- /dev/null +++ b/src/python/pyproject.toml @@ -0,0 +1,9 @@ +[tool.ruff] +src = ["."] +fix = true + +[tool.ruff.lint] +extend-select = ["I"] + +[tool.ruff.lint.isort] +known-first-party = ["pymole"] diff --git a/src/python/test/__init__.py b/src/python/test/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/python/test/dim1d/__init__.py b/src/python/test/dim1d/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/python/test/dim1d/test_divergence.py b/src/python/test/dim1d/test_divergence.py new file mode 100644 index 000000000..9bde1f9a5 --- /dev/null +++ b/src/python/test/dim1d/test_divergence.py @@ -0,0 +1,73 @@ +import numpy as np +import pytest + +from pymole import Divergence1D +from pymole.dim1d.grid import Grid1D + + +@pytest.fixture +def a_second_order_divergence(): + grid = Grid1D(0.0, 1.2, 6) + return Divergence1D(grid, 2) + + +class TestDivergence1DConstruction: + def test_minimum_cells_for_divergence(self): + grid_too_small = Grid1D(0.0, 0.4, 4) + with pytest.raises(ValueError): + Divergence1D(grid_too_small, 2) + + grid_ok = Grid1D(0.0, 0.5, 5) + divergence = Divergence1D(grid_ok, 2) + assert divergence.x_grid.num_cells == 5 + + def test_unsupported_accuracy_order_raises_error(self): + grid = Grid1D(0.0, 1.0, 10) + with pytest.raises(ValueError): + divergence = Divergence1D(grid, 4) + _ = divergence.matrix + + +class TestDivergence1DMatrix: + def test_matrix_has_correct_shape(self, a_second_order_divergence): + D = a_second_order_divergence.matrix + num_cells = a_second_order_divergence.x_grid.num_cells + + assert D.shape == (num_cells + 2, num_cells + 1) + + def test_constant_field_divergence_is_zero(self, a_second_order_divergence): + num_cells = a_second_order_divergence.x_grid.num_cells + v_constant = np.ones(num_cells + 1) + D = a_second_order_divergence.matrix + div_v = D @ v_constant + + assert np.allclose(div_v, 0.0, atol=1e-12) + + +class TestDivergence1DSecondOrder: + def test_second_order_boundary_cells(self): + grid = Grid1D(0.0, 1.25, 5) + divergence = Divergence1D(grid, 2) + num_cells = grid.num_cells + + D = divergence.matrix.toarray() + + expected_first_row = np.zeros(num_cells + 1) + expected_last_row = np.zeros(num_cells + 1) + + assert np.allclose(D[0, :], expected_first_row, atol=1e-12) + assert np.allclose(D[-1, :], expected_last_row, atol=1e-12) + + def test_second_order_interior_stencil(self, a_second_order_divergence): + num_cells = a_second_order_divergence.x_grid.num_cells + delta_x = a_second_order_divergence.x_grid.cell_size + + D = a_second_order_divergence.matrix.toarray() + + for i in range(1, num_cells + 1): + row = D[i, :] + non_zero_indices = np.where(np.abs(row) > 1e-12)[0] + + assert len(non_zero_indices) == 2 + assert row[i - 1] == pytest.approx(-1.0 / delta_x) + assert row[i] == pytest.approx(1.0 / delta_x) diff --git a/src/python/test/dim1d/test_gradient.py b/src/python/test/dim1d/test_gradient.py new file mode 100644 index 000000000..0c7aec0be --- /dev/null +++ b/src/python/test/dim1d/test_gradient.py @@ -0,0 +1,80 @@ +import numpy as np +import pytest + +from pymole import Gradient1D +from pymole.dim1d.grid import Grid1D + + +@pytest.fixture +def a_second_order_gradient(): + grid = Grid1D(0.0, 3.0, 15) + return Gradient1D(grid, 2) + + +class TestGradient1DConstruction: + def test_minimum_cells_for_gradient(self): + grid_too_small = Grid1D(0.0, 0.3, 3) + with pytest.raises(ValueError): + Gradient1D(grid_too_small, 2) + + grid_ok = Grid1D(0.0, 0.4, 4) + gradient = Gradient1D(grid_ok, 2) + + assert gradient.x_grid.num_cells == 4 + + def test_unsupported_accuracy_order_raises_error(self): + grid = Grid1D(0.0, 1.0, 10) + with pytest.raises(ValueError): + gradient = Gradient1D(grid, 4) + _ = gradient.matrix + + +class TestGradient1DMatrix: + def test_matrix_has_correct_shape(self, a_second_order_gradient): + G = a_second_order_gradient.matrix + num_cells = a_second_order_gradient.x_grid.num_cells + + assert G.shape == (num_cells + 1, num_cells + 2) + + def test_constant_field_gradient_is_zero(self, a_second_order_gradient): + num_cells = a_second_order_gradient.x_grid.num_cells + f_constant = np.ones(num_cells + 2) + G = a_second_order_gradient.matrix + grad_f = G @ f_constant + + assert np.allclose(grad_f, 0.0, atol=1e-12) + + def test_row_sums_are_zero(self, a_second_order_gradient): + G = a_second_order_gradient.matrix + row_sums = G.sum(axis=1).A1 + + assert np.allclose(row_sums, 0.0, atol=1e-12) + + +class TestGradient1DSecondOrder: + def test_second_order_boundary_stencil(self): + grid = Grid1D(0.0, 1.0, 4) + gradient = Gradient1D(grid, 2) + delta_x = grid.cell_size + + G = gradient.matrix.toarray() + + expected_first_row = np.array([-8.0 / 3.0, 3.0, -1.0 / 3.0, 0, 0, 0]) / delta_x + expected_last_row = np.array([0, 0, 0, 1.0 / 3.0, -3.0, 8.0 / 3.0]) / delta_x + + assert np.allclose(G[0, :], expected_first_row, atol=1e-12) + assert np.allclose(G[-1, :], expected_last_row, atol=1e-12) + + def test_second_order_interior_stencil(self, a_second_order_gradient): + num_cells = a_second_order_gradient.x_grid.num_cells + delta_x = a_second_order_gradient.x_grid.cell_size + + G = a_second_order_gradient.matrix.toarray() + + for i in range(1, num_cells): + row = G[i, :] + non_zero_indices = np.where(np.abs(row) > 1e-12)[0] + + assert len(non_zero_indices) == 2 + assert row[i] == pytest.approx(-1.0 / delta_x) + assert row[i + 1] == pytest.approx(1.0 / delta_x) diff --git a/src/python/test/dim1d/test_grid.py b/src/python/test/dim1d/test_grid.py new file mode 100644 index 000000000..d4d55ccf6 --- /dev/null +++ b/src/python/test/dim1d/test_grid.py @@ -0,0 +1,85 @@ +import pytest + +from pymole.dim1d.grid import Grid1D + + +@pytest.fixture +def basic_grid(): + return Grid1D(0.0, 1.0, 10) + + +class TestGrid1DConstruction: + def test_valid_grid_creation(self, basic_grid): + assert basic_grid.lower_bound == pytest.approx(0.0) + assert basic_grid.upper_bound == pytest.approx(1.0) + assert basic_grid.num_cells == 10 + assert basic_grid.cell_size == pytest.approx(0.1) + + def test_lower_bound_greater_than_upper_bound_raises_error(self): + with pytest.raises(ValueError): + Grid1D(1.0, 0.0, 10) + + def test_num_cells_less_than_one_raises_error(self): + with pytest.raises(ValueError): + Grid1D(0.0, 1.0, 0) + + +class TestGrid1DFromNumberOfCells: + def test_from_number_of_cells_basic(self): + expected_num_cells = 13 + grid = Grid1D.from_number_of_cells(expected_num_cells) + + assert grid.lower_bound == pytest.approx(0.0) + assert grid.upper_bound == pytest.approx(1.0) + assert grid.num_cells == expected_num_cells + assert grid.cell_size == pytest.approx(1.0 / expected_num_cells) + + def test_zero_number_of_cells_raises_error(self): + with pytest.raises(ValueError): + Grid1D.from_number_of_cells(0) + + def test_negative_number_of_cells_raises_error(self): + with pytest.raises(ValueError): + Grid1D.from_number_of_cells(-5) + + +class TestGrid1DProperties: + def test_cell_size_computed(self): + grid = Grid1D(0.0, 2.0, 4) + + assert grid.cell_size == pytest.approx(0.5) + + +class TestGrid1DEquality: + def test_equal_grids(self, basic_grid): + grid2 = Grid1D(0.0, 1.0, 10) + + assert basic_grid == grid2 + + def test_unequal_grids_different_lower_bound(self, basic_grid): + grid2 = Grid1D(0.5, 1.0, 10) + + assert basic_grid != grid2 + + def test_unequal_grids_different_upper_bound(self, basic_grid): + grid2 = Grid1D(0.0, 2.0, 10) + + assert basic_grid != grid2 + + def test_unequal_grids_different_num_cells(self, basic_grid): + grid2 = Grid1D(0.0, 1.0, 20) + + assert basic_grid != grid2 + + def test_grid_not_equal_to_non_grid(self, basic_grid): + assert basic_grid != "not a grid" + + +class TestGrid1DRepr: + def test_repr_strings(self, basic_grid): + repr_str = repr(basic_grid) + + assert "Grid1D" in repr_str + assert "lower_bound=0.0" in repr_str + assert "upper_bound=1.0" in repr_str + assert "num_cells=10" in repr_str diff --git a/src/python/test/dim1d/test_laplacian.py b/src/python/test/dim1d/test_laplacian.py new file mode 100644 index 000000000..699ebae76 --- /dev/null +++ b/src/python/test/dim1d/test_laplacian.py @@ -0,0 +1,54 @@ +import numpy as np +import pytest + +from pymole import Laplacian1D +from pymole.dim1d.divergence import Divergence1D +from pymole.dim1d.gradient import Gradient1D +from pymole.dim1d.grid import Grid1D + + +@pytest.fixture +def a_second_order_laplacian(): + grid = Grid1D(0.0, 1.2, 6) + return Laplacian1D(grid, 2) + + +class TestLaplacian1DConstruction: + def test_minimum_cells_for_laplacian(self): + grid_too_small = Grid1D(0.0, 0.4, 4) + with pytest.raises(ValueError): + Laplacian1D(grid_too_small, 2) + + grid_ok = Grid1D(0.0, 0.5, 5) + laplacian = Laplacian1D(grid_ok, 2) + assert laplacian.x_grid.num_cells == 5 + + def test_unsupported_accuracy_order_raises_error(self): + grid = Grid1D(0.0, 1.0, 10) + with pytest.raises(ValueError): + laplacian = Laplacian1D(grid, 4) + _ = laplacian.matrix + + +class TestLaplacian1DMatrix: + def test_matrix_has_correct_shape(self, a_second_order_laplacian): + L = a_second_order_laplacian.matrix + num_cells = a_second_order_laplacian.x_grid.num_cells + + assert L.shape == (num_cells + 2, num_cells + 2) + + def test_constant_field_laplacian_is_zero(self, a_second_order_laplacian): + num_cells = a_second_order_laplacian.x_grid.num_cells + f_constant = np.ones(num_cells + 2) + L = a_second_order_laplacian.matrix + lap_f = L @ f_constant + + assert np.allclose(lap_f, 0.0, atol=1e-12) + + def test_laplacian_equals_divergence_times_gradient(self, a_second_order_laplacian): + grid = a_second_order_laplacian.x_grid + D = Divergence1D(grid, 2).matrix + G = Gradient1D(grid, 2).matrix + L = a_second_order_laplacian.matrix + + assert np.allclose((L - D @ G).toarray(), 0.0, atol=1e-12) diff --git a/src/python/test/dim1d/test_operator.py b/src/python/test/dim1d/test_operator.py new file mode 100644 index 000000000..77285bf1a --- /dev/null +++ b/src/python/test/dim1d/test_operator.py @@ -0,0 +1,47 @@ +from unittest.mock import MagicMock + +import pytest +from scipy import sparse + +from pymole.dim1d.grid import Grid1D +from pymole.dim1d.operator import MimeticOperator1D + + +class MockOperator(MimeticOperator1D): + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return accuracy_order + + def _construct_matrix(self) -> sparse.csr_matrix: + return MagicMock(spec=sparse.csr_matrix) + + +class TestMimeticOperatorConstruction: + def test_minimum_cells_requirement(self): + grid_too_small = Grid1D(0.0, 0.3, 3) + with pytest.raises(ValueError): + MockOperator(grid_too_small, 4) + + grid_ok = Grid1D(0.0, 0.4, 4) + op = MockOperator(grid_ok, 2) + + assert op.x_grid == grid_ok + assert op.accuracy_order == 2 + assert op.minimum_cells == 2 + + +class TestMimeticOperatorMatrix: + def test_matrix_is_sparse_csr(self): + grid = Grid1D(0.0, 0.4, 4) + op = MockOperator(grid, 2) + matrix = op.matrix + + assert isinstance(matrix, MagicMock) + + def test_matrix_cached_property(self): + grid = Grid1D(0.0, 0.4, 4) + op = MockOperator(grid, 2) + matrix1 = op.matrix + matrix2 = op.matrix + + assert matrix1 is matrix2 diff --git a/src/python/test/dim2d/__init__.py b/src/python/test/dim2d/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/python/test/dim2d/test_divergence.py b/src/python/test/dim2d/test_divergence.py new file mode 100644 index 000000000..2d9246d3a --- /dev/null +++ b/src/python/test/dim2d/test_divergence.py @@ -0,0 +1,57 @@ +import numpy as np +import pytest + +from pymole.dim1d.grid import Grid1D +from pymole.dim2d.divergence import Divergence2D + + +@pytest.fixture +def a_second_order_divergence(): + x_grid = Grid1D(0.0, 3.0, 6) + y_grid = Grid1D(0.0, 2.0, 5) + return Divergence2D(x_grid, y_grid, 2) + + +class TestDivergence2DConstruction: + def test_minimum_cells_for_divergence(self): + x_grid_too_small = Grid1D(0.0, 0.4, 4) + y_grid_ok = Grid1D(0.0, 1.0, 5) + with pytest.raises(ValueError): + Divergence2D(x_grid_too_small, y_grid_ok, 2) + + x_grid_ok = Grid1D(0.0, 0.5, 5) + y_grid_too_small = Grid1D(0.0, 0.4, 4) + with pytest.raises(ValueError): + Divergence2D(x_grid_ok, y_grid_too_small, 2) + + x_grid = Grid1D(0.0, 0.5, 5) + y_grid = Grid1D(0.0, 0.5, 5) + divergence = Divergence2D(x_grid, y_grid, 2) + + assert divergence.x_grid.num_cells == 5 + assert divergence.y_grid.num_cells == 5 + + def test_unsupported_accuracy_order_raises_error(self): + x_grid = Grid1D(0.0, 1.0, 10) + y_grid = Grid1D(0.0, 1.0, 10) + with pytest.raises(ValueError): + divergence = Divergence2D(x_grid, y_grid, 4) + _ = divergence.matrix + + +class TestDivergence2DMatrix: + def test_matrix_has_correct_shape(self, a_second_order_divergence): + m = a_second_order_divergence.x_grid.num_cells + n = a_second_order_divergence.y_grid.num_cells + D = a_second_order_divergence.matrix + + assert D.shape == ((m + 2) * (n + 2), 2 * m * n + m + n) + + def test_constant_field_divergence_is_zero(self, a_second_order_divergence): + m = a_second_order_divergence.x_grid.num_cells + n = a_second_order_divergence.y_grid.num_cells + v_constant = np.ones(2 * m * n + m + n) + D = a_second_order_divergence.matrix + div_v = D @ v_constant + + assert np.allclose(div_v, 0.0, atol=1e-12) diff --git a/src/python/test/dim2d/test_gradient.py b/src/python/test/dim2d/test_gradient.py new file mode 100644 index 000000000..6f1ac86e0 --- /dev/null +++ b/src/python/test/dim2d/test_gradient.py @@ -0,0 +1,63 @@ +import numpy as np +import pytest + +from pymole.dim1d.grid import Grid1D +from pymole.dim2d.gradient import Gradient2D + + +@pytest.fixture +def a_second_order_gradient(): + x_grid = Grid1D(0.0, 3.0, 6) + y_grid = Grid1D(0.0, 2.0, 5) + return Gradient2D(x_grid, y_grid, 2) + + +class TestGradient2DConstruction: + def test_minimum_cells_for_gradient(self): + x_grid_too_small = Grid1D(0.0, 0.3, 3) + y_grid_ok = Grid1D(0.0, 1.0, 4) + with pytest.raises(ValueError): + Gradient2D(x_grid_too_small, y_grid_ok, 2) + + x_grid_ok = Grid1D(0.0, 0.4, 4) + y_grid_too_small = Grid1D(0.0, 0.3, 3) + with pytest.raises(ValueError): + Gradient2D(x_grid_ok, y_grid_too_small, 2) + + x_grid = Grid1D(0.0, 0.4, 4) + y_grid = Grid1D(0.0, 0.4, 4) + gradient = Gradient2D(x_grid, y_grid, 2) + + assert gradient.x_grid.num_cells == 4 + assert gradient.y_grid.num_cells == 4 + + def test_unsupported_accuracy_order_raises_error(self): + x_grid = Grid1D(0.0, 1.0, 10) + y_grid = Grid1D(0.0, 1.0, 10) + with pytest.raises(ValueError): + gradient = Gradient2D(x_grid, y_grid, 4) + _ = gradient.matrix + + +class TestGradient2DMatrix: + def test_matrix_has_correct_shape(self, a_second_order_gradient): + m = a_second_order_gradient.x_grid.num_cells + n = a_second_order_gradient.y_grid.num_cells + G = a_second_order_gradient.matrix + + assert G.shape == (2 * m * n + m + n, (m + 2) * (n + 2)) + + def test_constant_field_gradient_is_zero(self, a_second_order_gradient): + m = a_second_order_gradient.x_grid.num_cells + n = a_second_order_gradient.y_grid.num_cells + f_constant = np.ones((m + 2) * (n + 2)) + G = a_second_order_gradient.matrix + grad_f = G @ f_constant + + assert np.allclose(grad_f, 0.0, atol=1e-12) + + def test_row_sums_are_zero(self, a_second_order_gradient): + G = a_second_order_gradient.matrix + row_sums = G.sum(axis=1).A1 + + assert np.allclose(row_sums, 0.0, atol=1e-12) diff --git a/src/python/test/dim2d/test_laplacian.py b/src/python/test/dim2d/test_laplacian.py new file mode 100644 index 000000000..5c336caaa --- /dev/null +++ b/src/python/test/dim2d/test_laplacian.py @@ -0,0 +1,68 @@ +import numpy as np +import pytest + +from pymole.dim1d.grid import Grid1D +from pymole.dim2d.divergence import Divergence2D +from pymole.dim2d.gradient import Gradient2D +from pymole.dim2d.laplacian import Laplacian2D + + +@pytest.fixture +def a_second_order_laplacian(): + x_grid = Grid1D(0.0, 3.0, 6) + y_grid = Grid1D(0.0, 2.0, 5) + return Laplacian2D(x_grid, y_grid, 2) + + +class TestLaplacian2DConstruction: + def test_minimum_cells_for_laplacian(self): + x_grid_too_small = Grid1D(0.0, 0.4, 4) + y_grid_ok = Grid1D(0.0, 1.0, 5) + with pytest.raises(ValueError): + Laplacian2D(x_grid_too_small, y_grid_ok, 2) + + x_grid_ok = Grid1D(0.0, 0.5, 5) + y_grid_too_small = Grid1D(0.0, 0.4, 4) + with pytest.raises(ValueError): + Laplacian2D(x_grid_ok, y_grid_too_small, 2) + + x_grid = Grid1D(0.0, 0.5, 5) + y_grid = Grid1D(0.0, 0.5, 5) + laplacian = Laplacian2D(x_grid, y_grid, 2) + + assert laplacian.x_grid.num_cells == 5 + assert laplacian.y_grid.num_cells == 5 + + def test_unsupported_accuracy_order_raises_error(self): + x_grid = Grid1D(0.0, 1.0, 10) + y_grid = Grid1D(0.0, 1.0, 10) + with pytest.raises(ValueError): + laplacian = Laplacian2D(x_grid, y_grid, 4) + _ = laplacian.matrix + + +class TestLaplacian2DMatrix: + def test_matrix_has_correct_shape(self, a_second_order_laplacian): + m = a_second_order_laplacian.x_grid.num_cells + n = a_second_order_laplacian.y_grid.num_cells + L = a_second_order_laplacian.matrix + + assert L.shape == ((m + 2) * (n + 2), (m + 2) * (n + 2)) + + def test_constant_field_laplacian_is_zero(self, a_second_order_laplacian): + m = a_second_order_laplacian.x_grid.num_cells + n = a_second_order_laplacian.y_grid.num_cells + f_constant = np.ones((m + 2) * (n + 2)) + L = a_second_order_laplacian.matrix + lap_f = L @ f_constant + + assert np.allclose(lap_f, 0.0, atol=1e-12) + + def test_laplacian_equals_divergence_times_gradient(self, a_second_order_laplacian): + x_grid = a_second_order_laplacian.x_grid + y_grid = a_second_order_laplacian.y_grid + D = Divergence2D(x_grid, y_grid, 2).matrix + G = Gradient2D(x_grid, y_grid, 2).matrix + L = a_second_order_laplacian.matrix + + assert np.allclose((L - D @ G).toarray(), 0.0, atol=1e-12) diff --git a/src/python/test/dim2d/test_operator.py b/src/python/test/dim2d/test_operator.py new file mode 100644 index 000000000..58ef9087d --- /dev/null +++ b/src/python/test/dim2d/test_operator.py @@ -0,0 +1,59 @@ +from unittest.mock import MagicMock + +import pytest +from scipy import sparse + +from pymole.dim1d.grid import Grid1D +from pymole.dim2d.operator import MimeticOperator2D + + +class MockOperator2D(MimeticOperator2D): + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return accuracy_order + + def _construct_matrix(self) -> sparse.csr_matrix: + return MagicMock(spec=sparse.csr_matrix) + + +@pytest.fixture +def x_grid(): + return Grid1D(0.0, 1.0, 4) + + +@pytest.fixture +def y_grid(): + return Grid1D(0.0, 1.0, 4) + + +class TestMimeticOperator2DConstruction: + def test_x_grid_too_small_raises_error(self, y_grid): + grid_too_small = Grid1D(0.0, 0.3, 3) + with pytest.raises(ValueError): + MockOperator2D(grid_too_small, y_grid, 4) + + def test_y_grid_too_small_raises_error(self, x_grid): + grid_too_small = Grid1D(0.0, 0.3, 3) + with pytest.raises(ValueError): + MockOperator2D(x_grid, grid_too_small, 4) + + def test_stores_grids_and_accuracy_order(self, x_grid, y_grid): + op = MockOperator2D(x_grid, y_grid, 2) + + assert op.x_grid == x_grid + assert op.y_grid == y_grid + assert op.accuracy_order == 2 + + def test_minimum_cells_property(self, x_grid, y_grid): + op = MockOperator2D(x_grid, y_grid, 2) + + assert op.minimum_cells == 2 + + +class TestMimeticOperator2DMatrix: + def test_matrix_cached_property(self, x_grid, y_grid): + op = MockOperator2D(x_grid, y_grid, 2) + matrix1 = op.matrix + matrix2 = op.matrix + + assert matrix1 is matrix2 diff --git a/src/python/test/dim3d/__init__.py b/src/python/test/dim3d/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/python/test/dim3d/test_gradient.py b/src/python/test/dim3d/test_gradient.py new file mode 100644 index 000000000..1087a1616 --- /dev/null +++ b/src/python/test/dim3d/test_gradient.py @@ -0,0 +1,76 @@ +import numpy as np +import pytest + +from pymole.dim1d.grid import Grid1D +from pymole.dim3d.gradient import Gradient3D + + +@pytest.fixture +def a_second_order_gradient(): + x_grid = Grid1D(0.0, 3.0, 6) + y_grid = Grid1D(0.0, 2.0, 5) + z_grid = Grid1D(0.0, 1.5, 4) + return Gradient3D(x_grid, y_grid, z_grid, 2) + + +class TestGradient3DConstruction: + def test_minimum_cells_for_gradient(self): + x_grid_too_small = Grid1D(0.0, 0.3, 3) + y_grid_ok = Grid1D(0.0, 1.0, 4) + z_grid_ok = Grid1D(0.0, 1.0, 4) + with pytest.raises(ValueError): + Gradient3D(x_grid_too_small, y_grid_ok, z_grid_ok, 2) + + x_grid_ok = Grid1D(0.0, 0.4, 4) + y_grid_too_small = Grid1D(0.0, 0.3, 3) + with pytest.raises(ValueError): + Gradient3D(x_grid_ok, y_grid_too_small, z_grid_ok, 2) + + z_grid_too_small = Grid1D(0.0, 0.3, 3) + with pytest.raises(ValueError): + Gradient3D(x_grid_ok, y_grid_ok, z_grid_too_small, 2) + + x_grid = Grid1D(0.0, 0.4, 4) + y_grid = Grid1D(0.0, 0.4, 4) + z_grid = Grid1D(0.0, 0.4, 4) + gradient = Gradient3D(x_grid, y_grid, z_grid, 2) + + assert gradient.x_grid.num_cells == 4 + assert gradient.y_grid.num_cells == 4 + assert gradient.z_grid.num_cells == 4 + + def test_unsupported_accuracy_order_raises_error(self): + x_grid = Grid1D(0.0, 1.0, 10) + y_grid = Grid1D(0.0, 1.0, 10) + z_grid = Grid1D(0.0, 1.0, 10) + with pytest.raises(ValueError): + gradient = Gradient3D(x_grid, y_grid, z_grid, 4) + _ = gradient.matrix + + +class TestGradient3DMatrix: + def test_matrix_has_correct_shape(self, a_second_order_gradient): + m = a_second_order_gradient.x_grid.num_cells + n = a_second_order_gradient.y_grid.num_cells + o = a_second_order_gradient.z_grid.num_cells + G = a_second_order_gradient.matrix + + expected_rows = 3 * m * n * o + m * n + m * o + n * o + expected_cols = (m + 2) * (n + 2) * (o + 2) + assert G.shape == (expected_rows, expected_cols) + + def test_constant_field_gradient_is_zero(self, a_second_order_gradient): + m = a_second_order_gradient.x_grid.num_cells + n = a_second_order_gradient.y_grid.num_cells + o = a_second_order_gradient.z_grid.num_cells + f_constant = np.ones((m + 2) * (n + 2) * (o + 2)) + G = a_second_order_gradient.matrix + grad_f = G @ f_constant + + assert np.allclose(grad_f, 0.0, atol=1e-12) + + def test_row_sums_are_zero(self, a_second_order_gradient): + G = a_second_order_gradient.matrix + row_sums = G.sum(axis=1).A1 + + assert np.allclose(row_sums, 0.0, atol=1e-12) diff --git a/src/python/test/dim3d/test_operator.py b/src/python/test/dim3d/test_operator.py new file mode 100644 index 000000000..91e10e5b1 --- /dev/null +++ b/src/python/test/dim3d/test_operator.py @@ -0,0 +1,70 @@ +from unittest.mock import MagicMock + +import pytest +from scipy import sparse + +from pymole.dim1d.grid import Grid1D +from pymole.dim3d.operator import MimeticOperator3D + + +class MockOperator3D(MimeticOperator3D): + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return accuracy_order + + def _construct_matrix(self) -> sparse.csr_matrix: + return MagicMock(spec=sparse.csr_matrix) + + +@pytest.fixture +def x_grid(): + return Grid1D(-9.4, 2.6, 12) + + +@pytest.fixture +def y_grid(): + return Grid1D(0.0, 1.0, 7) + + +@pytest.fixture +def z_grid(): + return Grid1D(0.0, 1.0, 11) + + +class TestMimeticOperator3DConstruction: + def test_x_grid_too_small_raises_error(self, y_grid, z_grid): + grid_too_small = Grid1D(0.0, 0.3, 3) + with pytest.raises(ValueError): + MockOperator3D(grid_too_small, y_grid, z_grid, 4) + + def test_y_grid_too_small_raises_error(self, x_grid, z_grid): + grid_too_small = Grid1D(0.0, 0.3, 3) + with pytest.raises(ValueError): + MockOperator3D(x_grid, grid_too_small, z_grid, 4) + + def test_z_grid_too_small_raises_error(self, x_grid, y_grid): + grid_too_small = Grid1D(0.0, 0.3, 3) + with pytest.raises(ValueError): + MockOperator3D(x_grid, y_grid, grid_too_small, 4) + + def test_stores_grids_and_accuracy_order(self, x_grid, y_grid, z_grid): + op = MockOperator3D(x_grid, y_grid, z_grid, 2) + + assert op.x_grid == x_grid + assert op.y_grid == y_grid + assert op.z_grid == z_grid + assert op.accuracy_order == 2 + + def test_minimum_cells_property(self, x_grid, y_grid, z_grid): + op = MockOperator3D(x_grid, y_grid, z_grid, 2) + + assert op.minimum_cells == 2 + + +class TestMimeticOperator3DMatrix: + def test_matrix_cached_property(self, x_grid, y_grid, z_grid): + op = MockOperator3D(x_grid, y_grid, z_grid, 2) + matrix1 = op.matrix + matrix2 = op.matrix + + assert matrix1 is matrix2 From b2fe6c5ad9db1d3eab47c6dccd6ebd6b8ccae27c Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 31 May 2026 10:21:31 -0700 Subject: [PATCH 02/66] feat: implement mimetic operator classes and grid generation utilities for 1D, 2D, and 3D dimensions --- examples/python/simple_operators_check.py | 4 +- src/python/pymole/dim1d/operator.py | 2 +- src/python/pymole/dim2d/operator.py | 2 +- src/python/pymole/dim3d/operator.py | 2 +- src/python/pymole/grid.py | 80 +++++++++++++++++++ .../{operator.py => mimetic_operator.py} | 0 6 files changed, 86 insertions(+), 4 deletions(-) create mode 100644 src/python/pymole/grid.py rename src/python/pymole/{operator.py => mimetic_operator.py} (100%) diff --git a/examples/python/simple_operators_check.py b/examples/python/simple_operators_check.py index eceae4ac8..caa8ec5f9 100644 --- a/examples/python/simple_operators_check.py +++ b/examples/python/simple_operators_check.py @@ -1,5 +1,7 @@ +from pathlib import Path import sys -sys.path.insert(0, "src/python") + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent.parent / "src" / "python")) import numpy as np from numpy.typing import NDArray diff --git a/src/python/pymole/dim1d/operator.py b/src/python/pymole/dim1d/operator.py index 255b832b1..6809967df 100644 --- a/src/python/pymole/dim1d/operator.py +++ b/src/python/pymole/dim1d/operator.py @@ -2,7 +2,7 @@ from typing import TYPE_CHECKING -from pymole.operator import MimeticOperator +from pymole.mimetic_operator import MimeticOperator if TYPE_CHECKING: from .grid import Grid1D diff --git a/src/python/pymole/dim2d/operator.py b/src/python/pymole/dim2d/operator.py index 5e847d93a..de1161f91 100644 --- a/src/python/pymole/dim2d/operator.py +++ b/src/python/pymole/dim2d/operator.py @@ -2,7 +2,7 @@ from typing import TYPE_CHECKING -from pymole.operator import MimeticOperator +from pymole.mimetic_operator import MimeticOperator if TYPE_CHECKING: from pymole.dim1d.grid import Grid1D diff --git a/src/python/pymole/dim3d/operator.py b/src/python/pymole/dim3d/operator.py index 43a5963ef..a58003305 100644 --- a/src/python/pymole/dim3d/operator.py +++ b/src/python/pymole/dim3d/operator.py @@ -2,7 +2,7 @@ from typing import TYPE_CHECKING -from pymole.operator import MimeticOperator +from pymole.mimetic_operator import MimeticOperator if TYPE_CHECKING: from pymole.dim1d.grid import Grid1D diff --git a/src/python/pymole/grid.py b/src/python/pymole/grid.py new file mode 100644 index 000000000..38f160e57 --- /dev/null +++ b/src/python/pymole/grid.py @@ -0,0 +1,80 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' +# pyrefly: ignore [missing-import] +import numpy as np +from typing import Tuple, Union, List + +class Grid: + """ + A utility class to generate 1D, 2D, or 3D grids using NumPy's meshgrid. + """ + + @staticmethod + def generate( + p1: Union[float, int, np.ndarray, List[float], Tuple[float, ...]], + p2: Union[float, int, np.ndarray, List[float], Tuple[float, ...]], + n: Union[int, np.ndarray, List[int], Tuple[int, ...]], + indexing: str = 'xy' + ) -> Union[np.ndarray, Tuple[np.ndarray, ...]]: + """ + Generates a 1D, 2D, or 3D grid with n evenly spaced points between p1 and p2. + + Args: + p1: Starting grid point (scalar or 1D array/sequence of length 1, 2, or 3). + p2: Ending grid point (scalar or 1D array/sequence of length 1, 2, or 3). + n: Number of evenly spaced points (scalar integer or sequence of length 1, 2, or 3). + indexing: 'xy' (default, MATLAB-style) or 'ij' (Matrix-style). + + Returns: + A single array (for 1D) or a tuple of arrays (for 2D/3D). + """ + p1_arr = np.atleast_1d(p1) + p2_arr = np.atleast_1d(p2) + + if p1_arr.ndim > 1 or p2_arr.ndim > 1: + raise ValueError("Starting and ending grid points must be scalars or 1D sequences.") + if p1_arr.shape != p2_arr.shape: + raise ValueError("p1 and p2 must have the same shape/dimension.") + + dims = len(p1_arr) + if not (1 <= dims <= 3): + raise ValueError("Grid supports 1 to 3 dimensions only.") + + if isinstance(n, (int, np.integer)): + n_arr = np.full(dims, n, dtype=int) + else: + n_arr = np.atleast_1d(n).astype(int) + if n_arr.ndim > 1 or len(n_arr) != dims: + raise ValueError( + f"n must be a single integer or a sequence/array of length {dims}." + ) + + coords = [np.linspace(p1_arr[i], p2_arr[i], n_arr[i]) for i in range(dims)] + + if dims == 1: + return coords[0] + + return np.meshgrid(*coords, indexing=indexing) + +if __name__ == "__main__": + # Example usage: + # 1D: p1=0, p2=10, n=5 + grid1d = Grid.generate(0, 10, 5) + print(f"1D Grid shape: {grid1d.shape}") + print(f"1D Grid points: {grid1d}") + + # 2D: p1=[0, 0], p2=[10, 20], n=[5, 3] + X, Y = Grid.generate([0, 0], [10, 20], [5, 3]) + print(f"2D Grid X shape: {X.shape}, Y shape: {Y.shape}") + print(f"2D Grid points for X:\n {X}") + print(f"2D Grid points for Y:\n {Y}") + + # 3D: p1=[0, 0, 0], p2=[10, 20, 30], n=[5, 3, 4] + X, Y, Z = Grid.generate([0, 0, 0], [10, 20, 30], [5, 3, 4]) + print(f"3D Grid X shape: {X.shape}, Y shape: {Y.shape}, Z shape: {Z.shape}") + print(f"2D Grid points for X:\n {X}") + print(f"2D Grid points for Y:\n {Y}") + print(f"2D Grid points for Z:\n {Z}") diff --git a/src/python/pymole/operator.py b/src/python/pymole/mimetic_operator.py similarity index 100% rename from src/python/pymole/operator.py rename to src/python/pymole/mimetic_operator.py From 5deb8d831ec47f1fc33a6522e0d191117e6547c3 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Fri, 5 Jun 2026 17:50:05 -0700 Subject: [PATCH 03/66] feat: add Grid utility class for 1D-3D mesh generation and corresponding test suite --- src/python/pymole/grid.py | 59 ++++++++++---- src/python/test/test_grid_utility.py | 112 +++++++++++++++++++++++++++ 2 files changed, 155 insertions(+), 16 deletions(-) create mode 100644 src/python/test/test_grid_utility.py diff --git a/src/python/pymole/grid.py b/src/python/pymole/grid.py index 38f160e57..61db9e9ef 100644 --- a/src/python/pymole/grid.py +++ b/src/python/pymole/grid.py @@ -16,21 +16,26 @@ class Grid: def generate( p1: Union[float, int, np.ndarray, List[float], Tuple[float, ...]], p2: Union[float, int, np.ndarray, List[float], Tuple[float, ...]], - n: Union[int, np.ndarray, List[int], Tuple[int, ...]], + shape: Union[int, np.ndarray, List[int], Tuple[int, ...], None] = None, + spacing: Union[float, np.ndarray, List[float], Tuple[float, ...], None] = None, indexing: str = 'xy' ) -> Union[np.ndarray, Tuple[np.ndarray, ...]]: """ - Generates a 1D, 2D, or 3D grid with n evenly spaced points between p1 and p2. + Generates a 1D, 2D, or 3D grid using shape (number of points/cells) or spacing (physical mesh size). Args: p1: Starting grid point (scalar or 1D array/sequence of length 1, 2, or 3). p2: Ending grid point (scalar or 1D array/sequence of length 1, 2, or 3). - n: Number of evenly spaced points (scalar integer or sequence of length 1, 2, or 3). + shape: Number of evenly spaced points (scalar integer or sequence of length 1, 2, or 3). + spacing: Physical mesh size (scalar float/int or sequence of length 1, 2, or 3). indexing: 'xy' (default, MATLAB-style) or 'ij' (Matrix-style). Returns: A single array (for 1D) or a tuple of arrays (for 2D/3D). """ + if (shape is None) == (spacing is None): + raise ValueError("Exactly one of 'shape' or 'spacing' must be specified.") + p1_arr = np.atleast_1d(p1) p2_arr = np.atleast_1d(p2) @@ -43,14 +48,36 @@ def generate( if not (1 <= dims <= 3): raise ValueError("Grid supports 1 to 3 dimensions only.") - if isinstance(n, (int, np.integer)): - n_arr = np.full(dims, n, dtype=int) + if shape is not None: + if isinstance(shape, (int, np.integer)): + shape_arr = np.full(dims, shape, dtype=int) + else: + shape_arr = np.atleast_1d(shape).astype(int) + if shape_arr.ndim > 1 or len(shape_arr) != dims: + raise ValueError( + f"shape must be a single integer or a sequence/array of length {dims}." + ) + n_arr = shape_arr else: - n_arr = np.atleast_1d(n).astype(int) - if n_arr.ndim > 1 or len(n_arr) != dims: - raise ValueError( - f"n must be a single integer or a sequence/array of length {dims}." - ) + if isinstance(spacing, (int, float, np.number)): + spacing_arr = np.full(dims, spacing, dtype=float) + else: + spacing_arr = np.atleast_1d(spacing).astype(float) + if spacing_arr.ndim > 1 or len(spacing_arr) != dims: + raise ValueError( + f"spacing must be a single float or a sequence/array of length {dims}." + ) + + if np.any(spacing_arr <= 0): + raise ValueError("spacing values must be greater than zero.") + + n_arr = np.empty(dims, dtype=int) + for i in range(dims): + length = abs(p2_arr[i] - p1_arr[i]) + num_cells = int(round(length / spacing_arr[i])) + if num_cells < 1: + num_cells = 1 + n_arr[i] = num_cells + 1 coords = [np.linspace(p1_arr[i], p2_arr[i], n_arr[i]) for i in range(dims)] @@ -61,19 +88,19 @@ def generate( if __name__ == "__main__": # Example usage: - # 1D: p1=0, p2=10, n=5 - grid1d = Grid.generate(0, 10, 5) + # 1D: p1=0, p2=10, shape=5 + grid1d = Grid.generate(0, 10, shape=5) print(f"1D Grid shape: {grid1d.shape}") print(f"1D Grid points: {grid1d}") - # 2D: p1=[0, 0], p2=[10, 20], n=[5, 3] - X, Y = Grid.generate([0, 0], [10, 20], [5, 3]) + # 2D: p1=[0, 0], p2=[10, 20], spacing=[2.5, 10.0] + X, Y = Grid.generate([0, 0], [10, 20], spacing=[2.5, 10.0]) print(f"2D Grid X shape: {X.shape}, Y shape: {Y.shape}") print(f"2D Grid points for X:\n {X}") print(f"2D Grid points for Y:\n {Y}") - # 3D: p1=[0, 0, 0], p2=[10, 20, 30], n=[5, 3, 4] - X, Y, Z = Grid.generate([0, 0, 0], [10, 20, 30], [5, 3, 4]) + # 3D: p1=[0, 0, 0], p2=[10, 20, 30], shape=[5, 3, 4] + X, Y, Z = Grid.generate([0, 0, 0], [10, 20, 30], shape=[5, 3, 4]) print(f"3D Grid X shape: {X.shape}, Y shape: {Y.shape}, Z shape: {Z.shape}") print(f"2D Grid points for X:\n {X}") print(f"2D Grid points for Y:\n {Y}") diff --git a/src/python/test/test_grid_utility.py b/src/python/test/test_grid_utility.py new file mode 100644 index 000000000..753e22c63 --- /dev/null +++ b/src/python/test/test_grid_utility.py @@ -0,0 +1,112 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' +import pytest +import numpy as np +import sys +from pathlib import Path +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) + +from pymole.grid import Grid + +def test_shape_grid_1d(): + # Test 1D grid generation using shape + grid = Grid.generate(0, 10, shape=5) + expected = np.array([0.0, 2.5, 5.0, 7.5, 10.0]) + assert np.allclose(grid, expected) + assert grid.shape == (5,) + +def test_spacing_grid_1d(): + # Test 1D grid generation using spacing + grid = Grid.generate(0, 10, spacing=2.0) + expected = np.array([0.0, 2.0, 4.0, 6.0, 8.0, 10.0]) + assert np.allclose(grid, expected) + assert grid.shape == (6,) + +def test_spacing_grid_1d_non_divisible(): + # Test 1D grid with non-divisible spacing (should round to nearest number of cells) + # Length is 10, spacing is 3. 10/3 = 3.33 -> rounds to 3 cells (4 points) + grid = Grid.generate(0, 10, spacing=3.0) + expected = np.linspace(0, 10, 4) + assert np.allclose(grid, expected) + assert grid.shape == (4,) + +def test_shape_grid_2d(): + # Test 2D grid generation using shape + X, Y = Grid.generate([0, 0], [10, 20], shape=[5, 3]) + assert X.shape == (3, 5) # for 'xy' indexing, shape of X is (Ny, Nx) + assert Y.shape == (3, 5) + + # Check boundaries + assert X[0, 0] == 0.0 + assert X[0, -1] == 10.0 + assert Y[0, 0] == 0.0 + assert Y[-1, 0] == 20.0 + +def test_spacing_grid_2d(): + # Test 2D grid generation using spacing + X, Y = Grid.generate([0, 0], [10, 20], spacing=[2.5, 10.0]) + # x spacing 2.5 -> 10/2.5 = 4 cells -> 5 points + # y spacing 10.0 -> 20/10 = 2 cells -> 3 points + assert X.shape == (3, 5) + assert Y.shape == (3, 5) + + assert np.allclose(X[0, :], [0.0, 2.5, 5.0, 7.5, 10.0]) + assert np.allclose(Y[:, 0], [0.0, 10.0, 20.0]) + +def test_shape_grid_3d(): + # Test 3D grid generation using shape + X, Y, Z = Grid.generate([0, 0, 0], [10, 20, 30], shape=[5, 3, 4]) + # indexing 'xy' -> shape is (shape[1], shape[0], shape[2]) + # which is (3, 5, 4) + assert X.shape == (3, 5, 4) + assert Y.shape == (3, 5, 4) + assert Z.shape == (3, 5, 4) + +def test_spacing_grid_3d(): + # Test 3D grid generation using spacing + X, Y, Z = Grid.generate([0, 0, 0], [10, 20, 30], spacing=[2.5, 10.0, 10.0]) + # x spacing 2.5 -> 4 cells -> 5 points + # y spacing 10.0 -> 2 cells -> 3 points + # z spacing 10.0 -> 3 cells -> 4 points + # shape is (3, 5, 4) + assert X.shape == (3, 5, 4) + assert Y.shape == (3, 5, 4) + assert Z.shape == (3, 5, 4) + +def test_both_specified_raises_error(): + with pytest.raises(ValueError, match="Exactly one of 'shape' or 'spacing' must be specified"): + Grid.generate(0, 10, shape=5, spacing=2.0) + +def test_neither_specified_raises_error(): + with pytest.raises(ValueError, match="Exactly one of 'shape' or 'spacing' must be specified"): + Grid.generate(0, 10) + +def test_invalid_dimension_raises_error(): + # 4 dimensions is not supported + with pytest.raises(ValueError, match="Grid supports 1 to 3 dimensions only"): + Grid.generate([0, 0, 0, 0], [1, 1, 1, 1], shape=5) + +def test_mismatched_p1_p2_raises_error(): + with pytest.raises(ValueError, match="p1 and p2 must have the same shape/dimension"): + Grid.generate([0, 0], [1, 1, 1], shape=5) + +def test_mismatched_shape_dimension_raises_error(): + with pytest.raises(ValueError, match="shape must be a single integer or a sequence/array of length"): + Grid.generate([0, 0], [1, 1], shape=[5, 5, 5]) + +def test_mismatched_spacing_dimension_raises_error(): + with pytest.raises(ValueError, match="spacing must be a single float or a sequence/array of length"): + Grid.generate([0, 0], [1, 1], spacing=[0.5, 0.5, 0.5]) + +def test_invalid_spacing_value_raises_error(): + with pytest.raises(ValueError, match="spacing values must be greater than zero"): + Grid.generate(0, 10, spacing=0) + + with pytest.raises(ValueError, match="spacing values must be greater than zero"): + Grid.generate(0, 10, spacing=-1.0) + + with pytest.raises(ValueError, match="spacing values must be greater than zero"): + Grid.generate([0, 0], [10, 10], spacing=[1.0, -0.5]) From 3989e7fa5b7e20b64769d9d73d31746d8d0a1ddf Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Thu, 11 Jun 2026 15:02:42 -0700 Subject: [PATCH 04/66] feat: implement Grid utility class for 1D-3D mesh generation and add unit tests --- src/python/pymole/grid.py | 101 ++++++++++++++++++++++----- src/python/test/test_grid_utility.py | 73 ++++++++++++++++--- 2 files changed, 145 insertions(+), 29 deletions(-) diff --git a/src/python/pymole/grid.py b/src/python/pymole/grid.py index 61db9e9ef..6966599a1 100644 --- a/src/python/pymole/grid.py +++ b/src/python/pymole/grid.py @@ -5,13 +5,76 @@ ''' # pyrefly: ignore [missing-import] import numpy as np -from typing import Tuple, Union, List +from typing import Tuple, Union, List, Optional class Grid: """ - A utility class to generate 1D, 2D, or 3D grids using NumPy's meshgrid. + A class that encapsulates 1D, 2D, or 3D grids generated using NumPy's meshgrid, + providing access to individual dimension arrays and unpacking support. """ + def __init__( + self, + x: np.ndarray, + y: Optional[np.ndarray] = None, + z: Optional[np.ndarray] = None + ): + self._x = x + self._y = y + self._z = z + + @property + def x(self) -> np.ndarray: + return self._x + + @property + def y(self) -> Optional[np.ndarray]: + return self._y + + @property + def z(self) -> Optional[np.ndarray]: + return self._z + + @property + def X(self) -> np.ndarray: + return self._x + + @property + def Y(self) -> Optional[np.ndarray]: + return self._y + + @property + def Z(self) -> Optional[np.ndarray]: + return self._z + + @property + def ndim(self) -> int: + if self._z is not None: + return 3 + if self._y is not None: + return 2 + return 1 + + def __iter__(self): + """Allows unpacking the grid object, e.g., X, Y = Grid.generate(...)""" + if self.ndim == 1: + yield self.x + elif self.ndim == 2: + yield self.x + yield self.y + else: + yield self.x + yield self.y + yield self.z + + def __repr__(self) -> str: + if self.ndim == 1: + return f"Grid(ndim=1, x={self.x.shape})" + elif self.ndim == 2: + return f"Grid(ndim=2, x={self.x.shape}, y={self.y.shape})" + else: + return f"Grid(ndim=3, x={self.x.shape}, y={self.y.shape}, z={self.z.shape})" + @staticmethod def generate( p1: Union[float, int, np.ndarray, List[float], Tuple[float, ...]], @@ -19,7 +82,7 @@ def generate( shape: Union[int, np.ndarray, List[int], Tuple[int, ...], None] = None, spacing: Union[float, np.ndarray, List[float], Tuple[float, ...], None] = None, indexing: str = 'xy' - ) -> Union[np.ndarray, Tuple[np.ndarray, ...]]: + ) -> 'Grid': """ Generates a 1D, 2D, or 3D grid using shape (number of points/cells) or spacing (physical mesh size). @@ -31,7 +94,7 @@ def generate( indexing: 'xy' (default, MATLAB-style) or 'ij' (Matrix-style). Returns: - A single array (for 1D) or a tuple of arrays (for 2D/3D). + A Grid object encapsulating the coordinate arrays. """ if (shape is None) == (spacing is None): raise ValueError("Exactly one of 'shape' or 'spacing' must be specified.") @@ -82,26 +145,28 @@ def generate( coords = [np.linspace(p1_arr[i], p2_arr[i], n_arr[i]) for i in range(dims)] if dims == 1: - return coords[0] + return Grid(coords[0]) - return np.meshgrid(*coords, indexing=indexing) + mesh = np.meshgrid(*coords, indexing=indexing) + return Grid(*mesh) if __name__ == "__main__": # Example usage: # 1D: p1=0, p2=10, shape=5 grid1d = Grid.generate(0, 10, shape=5) - print(f"1D Grid shape: {grid1d.shape}") - print(f"1D Grid points: {grid1d}") + print(f"1D Grid object: {grid1d}") + print(f"1D Grid x shape: {grid1d.x.shape}") + print(f"1D Grid points: {grid1d.x}") # 2D: p1=[0, 0], p2=[10, 20], spacing=[2.5, 10.0] - X, Y = Grid.generate([0, 0], [10, 20], spacing=[2.5, 10.0]) - print(f"2D Grid X shape: {X.shape}, Y shape: {Y.shape}") - print(f"2D Grid points for X:\n {X}") - print(f"2D Grid points for Y:\n {Y}") + grid2d = Grid.generate([0, 0], [10, 20], spacing=[2.5, 10.0]) + print(f"2D Grid object: {grid2d}") + print(f"2D Grid X shape: {grid2d.X.shape}, Y shape: {grid2d.Y.shape}") + print(f"2D Grid points for X:\n {grid2d.X}") + print(f"2D Grid points for Y:\n {grid2d.Y}") - # 3D: p1=[0, 0, 0], p2=[10, 20, 30], shape=[5, 3, 4] - X, Y, Z = Grid.generate([0, 0, 0], [10, 20, 30], shape=[5, 3, 4]) - print(f"3D Grid X shape: {X.shape}, Y shape: {Y.shape}, Z shape: {Z.shape}") - print(f"2D Grid points for X:\n {X}") - print(f"2D Grid points for Y:\n {Y}") - print(f"2D Grid points for Z:\n {Z}") + # 3D: p1=[0, 0, 0], p2=[10, 20, 30], shape=[5, 3, 4] (using unpacking) + grid3d = Grid.generate([0, 0, 0], [10, 20, 30], shape=[5, 3, 4]) + print(f"3D Grid object: {grid3d}") + X, Y, Z = grid3d + print(f"Unpacked 3D Grid X shape: {X.shape}, Y shape: {Y.shape}, Z shape: {Z.shape}") diff --git a/src/python/test/test_grid_utility.py b/src/python/test/test_grid_utility.py index 753e22c63..12bbf9ee8 100644 --- a/src/python/test/test_grid_utility.py +++ b/src/python/test/test_grid_utility.py @@ -3,10 +3,11 @@ © 2008-2024 San Diego State University Research Foundation (SDSURF). See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. ''' +from pathlib import Path import pytest import numpy as np import sys -from pathlib import Path + sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) from pymole.grid import Grid @@ -14,28 +15,36 @@ def test_shape_grid_1d(): # Test 1D grid generation using shape grid = Grid.generate(0, 10, shape=5) + assert grid.ndim == 1 + assert grid.y is None + assert grid.z is None expected = np.array([0.0, 2.5, 5.0, 7.5, 10.0]) - assert np.allclose(grid, expected) - assert grid.shape == (5,) + assert np.allclose(grid.x, expected) + assert grid.x.shape == (5,) def test_spacing_grid_1d(): # Test 1D grid generation using spacing grid = Grid.generate(0, 10, spacing=2.0) + assert grid.ndim == 1 expected = np.array([0.0, 2.0, 4.0, 6.0, 8.0, 10.0]) - assert np.allclose(grid, expected) - assert grid.shape == (6,) + assert np.allclose(grid.x, expected) + assert grid.x.shape == (6,) def test_spacing_grid_1d_non_divisible(): # Test 1D grid with non-divisible spacing (should round to nearest number of cells) # Length is 10, spacing is 3. 10/3 = 3.33 -> rounds to 3 cells (4 points) grid = Grid.generate(0, 10, spacing=3.0) expected = np.linspace(0, 10, 4) - assert np.allclose(grid, expected) - assert grid.shape == (4,) + assert np.allclose(grid.x, expected) + assert grid.x.shape == (4,) def test_shape_grid_2d(): # Test 2D grid generation using shape - X, Y = Grid.generate([0, 0], [10, 20], shape=[5, 3]) + grid = Grid.generate([0, 0], [10, 20], shape=[5, 3]) + assert grid.ndim == 2 + assert grid.z is None + + X, Y = grid.x, grid.y assert X.shape == (3, 5) # for 'xy' indexing, shape of X is (Ny, Nx) assert Y.shape == (3, 5) @@ -45,11 +54,16 @@ def test_shape_grid_2d(): assert Y[0, 0] == 0.0 assert Y[-1, 0] == 20.0 + # Check uppercase properties + assert grid.X is grid.x + assert grid.Y is grid.y + def test_spacing_grid_2d(): # Test 2D grid generation using spacing - X, Y = Grid.generate([0, 0], [10, 20], spacing=[2.5, 10.0]) + grid = Grid.generate([0, 0], [10, 20], spacing=[2.5, 10.0]) # x spacing 2.5 -> 10/2.5 = 4 cells -> 5 points # y spacing 10.0 -> 20/10 = 2 cells -> 3 points + X, Y = grid.X, grid.Y assert X.shape == (3, 5) assert Y.shape == (3, 5) @@ -58,24 +72,61 @@ def test_spacing_grid_2d(): def test_shape_grid_3d(): # Test 3D grid generation using shape - X, Y, Z = Grid.generate([0, 0, 0], [10, 20, 30], shape=[5, 3, 4]) + grid = Grid.generate([0, 0, 0], [10, 20, 30], shape=[5, 3, 4]) + assert grid.ndim == 3 + X, Y, Z = grid.x, grid.y, grid.z # indexing 'xy' -> shape is (shape[1], shape[0], shape[2]) # which is (3, 5, 4) assert X.shape == (3, 5, 4) assert Y.shape == (3, 5, 4) assert Z.shape == (3, 5, 4) + # Check uppercase properties + assert grid.X is grid.x + assert grid.Y is grid.y + assert grid.Z is grid.z + def test_spacing_grid_3d(): # Test 3D grid generation using spacing - X, Y, Z = Grid.generate([0, 0, 0], [10, 20, 30], spacing=[2.5, 10.0, 10.0]) + grid = Grid.generate([0, 0, 0], [10, 20, 30], spacing=[2.5, 10.0, 10.0]) # x spacing 2.5 -> 4 cells -> 5 points # y spacing 10.0 -> 2 cells -> 3 points # z spacing 10.0 -> 3 cells -> 4 points # shape is (3, 5, 4) + X, Y, Z = grid.X, grid.Y, grid.Z assert X.shape == (3, 5, 4) assert Y.shape == (3, 5, 4) assert Z.shape == (3, 5, 4) +def test_grid_unpacking(): + # 1D + grid1d = Grid.generate(0, 10, shape=5) + x, = grid1d + assert np.allclose(x, grid1d.x) + + # 2D + grid2d = Grid.generate([0, 0], [10, 20], shape=[5, 3]) + x, y = grid2d + assert np.allclose(x, grid2d.x) + assert np.allclose(y, grid2d.y) + + # 3D + grid3d = Grid.generate([0, 0, 0], [10, 20, 30], shape=[5, 3, 4]) + x, y, z = grid3d + assert np.allclose(x, grid3d.x) + assert np.allclose(y, grid3d.y) + assert np.allclose(z, grid3d.z) + +def test_grid_repr(): + grid1d = Grid.generate(0, 10, shape=5) + assert repr(grid1d) == "Grid(ndim=1, x=(5,))" + + grid2d = Grid.generate([0, 0], [10, 20], shape=[5, 3]) + assert repr(grid2d) == "Grid(ndim=2, x=(3, 5), y=(3, 5))" + + grid3d = Grid.generate([0, 0, 0], [10, 20, 30], shape=[5, 3, 4]) + assert repr(grid3d) == "Grid(ndim=3, x=(3, 5, 4), y=(3, 5, 4), z=(3, 5, 4))" + def test_both_specified_raises_error(): with pytest.raises(ValueError, match="Exactly one of 'shape' or 'spacing' must be specified"): Grid.generate(0, 10, shape=5, spacing=2.0) From 038dc52f80157f751550f14d5a19f376abbcefe0 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 14 Jun 2026 10:09:55 -0700 Subject: [PATCH 05/66] feat: implement Grid and Divergence classes with associated unit tests --- src/python/pymole/__init__.py | 6 +- src/python/pymole/divergence.py | 159 +++++++++++++++++++++++++++ src/python/pymole/grid.py | 22 ++++ src/python/test/test_divergence.py | 170 +++++++++++++++++++++++++++++ 4 files changed, 356 insertions(+), 1 deletion(-) create mode 100644 src/python/pymole/divergence.py create mode 100644 src/python/test/test_divergence.py diff --git a/src/python/pymole/__init__.py b/src/python/pymole/__init__.py index a8be89aa0..21cb17e11 100644 --- a/src/python/pymole/__init__.py +++ b/src/python/pymole/__init__.py @@ -1,14 +1,18 @@ from .dim1d import Divergence1D, Gradient1D, Grid1D, Laplacian1D from .dim2d import Divergence2D, Gradient2D, Laplacian2D from .dim3d import Gradient3D +from .grid import Grid +from .divergence import Divergence -__version__ = "0.1.0" +__version__ = "0.1.1" __all__ = [ + "Divergence", "Divergence1D", "Divergence2D", "Gradient1D", "Gradient2D", "Gradient3D", + "Grid", "Grid1D", "Laplacian1D", "Laplacian2D", diff --git a/src/python/pymole/divergence.py b/src/python/pymole/divergence.py new file mode 100644 index 000000000..e1b054cb2 --- /dev/null +++ b/src/python/pymole/divergence.py @@ -0,0 +1,159 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' + +import numpy as np +from scipy import sparse + +class Divergence: + """ + Unified divergence operator supporting 1D and 2D Grid objects. + """ + + def __init__(self, grid, accuracy_order: int = 2): + self.grid = grid + self.accuracy_order = accuracy_order + + if self.grid.ndim not in (1, 2): + raise ValueError( + f"Divergence currently supports only 1D and 2D grids. " + f"Got ndim={self.grid.ndim}" + ) + + self._matrix = self._construct_matrix() + + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return 2 * accuracy_order + 1 + + @property + def matrix(self): + return self._matrix + + def _construct_matrix(self): + + if self.accuracy_order != 2: + raise ValueError( + f"Unsupported order of accuracy: {self.accuracy_order}" + ) + + match self.grid.ndim: + + case 1: + return self._construct_1d() + + case 2: + return self._construct_2d() + + case _: + raise ValueError( + f"Unsupported grid dimension: {self.grid.ndim}" + ) + + def _construct_1d(self): + + m = self.grid.num_cells + + if m < self.minimum_cells_for_order(self.accuracy_order): + raise ValueError( + f"Grid requires at least " + f"{self.minimum_cells_for_order(self.accuracy_order)} " + f"cells for order {self.accuracy_order}" + ) + + dx = self.grid.spacing + + return self._build_1d_divergence(m, dx) + + def _construct_2d(self): + + m, n = self.grid.num_cells + dx, dy = self.grid.spacing + + minimum_cells = self.minimum_cells_for_order( + self.accuracy_order + ) + + if m < minimum_cells or n < minimum_cells: + raise ValueError( + f"Both dimensions require at least " + f"{minimum_cells} cells." + ) + + return self._build_2d_divergence( + m, + n, + dx, + dy, + ) + + def _build_1d_divergence( + self, + m: int, + dx: float, + ): + + main_diag = np.ones(m + 1) + + sub_diag = np.full( + m + 1, + -1.0, + ) + + main_diag[0] = 0.0 + sub_diag[-1] = 0.0 + + D = sparse.diags( + [sub_diag, main_diag], + offsets=[-1, 0], + shape=(m + 2, m + 1), + ).tocsr() + + return D / dx + + def _build_2d_divergence( + self, + m: int, + n: int, + dx: float, + dy: float, + ): + + Dx = self._build_1d_divergence( + m, + dx, + ) + + Dy = self._build_1d_divergence( + n, + dy, + ) + + Im = sparse.eye( + m + 2, + format="csr", + )[:, 1:-1] + + In = sparse.eye( + n + 2, + format="csr", + )[:, 1:-1] + + Sx = sparse.kron( + In, + Dx, + format="csr", + ) + + Sy = sparse.kron( + Dy, + Im, + format="csr", + ) + + return sparse.hstack( + [Sx, Sy], + format="csr", + ) diff --git a/src/python/pymole/grid.py b/src/python/pymole/grid.py index 6966599a1..bc662827f 100644 --- a/src/python/pymole/grid.py +++ b/src/python/pymole/grid.py @@ -54,6 +54,28 @@ def ndim(self) -> int: if self._y is not None: return 2 return 1 + @property + def num_cells(self): + if self.ndim == 1: + return len(self.x) - 1 + + return tuple(s - 1 for s in self.x.shape) + + @property + def spacing(self): + if self.ndim == 1: + return self.x[1] - self.x[0] + + if self.ndim == 2: + dx = self.x[0, 1] - self.x[0, 0] + dy = self.y[1, 0] - self.y[0, 0] + return dx, dy + + dx = self.x[0, 0, 1] - self.x[0, 0, 0] + dy = self.y[0, 1, 0] - self.y[0, 0, 0] + dz = self.z[1, 0, 0] - self.z[0, 0, 0] + + return dx, dy, dz def __iter__(self): """Allows unpacking the grid object, e.g., X, Y = Grid.generate(...)""" diff --git a/src/python/test/test_divergence.py b/src/python/test/test_divergence.py new file mode 100644 index 000000000..b247745cd --- /dev/null +++ b/src/python/test/test_divergence.py @@ -0,0 +1,170 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' + +import numpy as np +import pytest +import sys +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) +#from pymole.grid import Grid +#from pymole.divergence import Divergence +from pymole import Grid, Divergence + +class TestDivergence: + + def test_unsupported_accuracy_order(self): + + grid = Grid.generate( + 0, + 10, + shape=10, + ) + + with pytest.raises(ValueError): + Divergence( + grid, + accuracy_order=4, + ) + + def test_grid_too_small_1d(self): + + grid = Grid.generate( + 0, + 10, + shape=5, + ) + + with pytest.raises(ValueError): + Divergence( + grid, + 2, + ) + + def test_grid_too_small_2d(self): + + grid = Grid.generate( + [0, 0], + [10, 10], + shape=[5, 10], + ) + + with pytest.raises(ValueError): + Divergence( + grid, + 2, + ) + + def test_1d_matrix_shape(self): + + grid = Grid.generate( + 0, + 1.2, + shape=7, + ) + + D = Divergence( + grid, + 2, + ).matrix + + m = grid.num_cells + + assert D.shape == ( + m + 2, + m + 1, + ) + + def test_2d_matrix_shape(self): + + grid = Grid.generate( + [0, 0], + [3, 2], + shape=[7, 6], + ) + + D = Divergence( + grid, + 2, + ).matrix + + m, n = grid.num_cells + + expected_columns = ( + 2 * m * n + + m + + n + ) + + expected_rows = ( + (m + 2) + * (n + 2) + ) + + assert D.shape == ( + expected_rows, + expected_columns, + ) + + def test_1d_constant_field(self): + + grid = Grid.generate( + 0.0, + 1.2, + shape=7, + ) + + D = Divergence( + grid, + 2, + ).matrix + + m = grid.num_cells + + v = np.ones( + m + 1 + ) + + div_v = D @ v + + assert np.allclose( + div_v, + 0.0, + atol=1e-12, + ) + + def test_2d_constant_field(self): + + grid = Grid.generate( + [0.0, 0.0], + [3.0, 2.0], + shape=[7, 6], + ) + + D = Divergence( + grid, + 2, + ).matrix + + m, n = grid.num_cells + + vector_size = ( + 2 * m * n + + m + + n + ) + + v = np.ones( + vector_size + ) + + div_v = D @ v + + assert np.allclose( + div_v, + 0.0, + atol=1e-12, + ) From a76c4aaad80d0801ed9fd89ec681b30f0a21e660 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 14 Jun 2026 19:10:14 -0700 Subject: [PATCH 06/66] feat: implement Gradient operator for 1D, 2D, and 3D grids with associated unit tests and CI configuration. --- .github/workflows/ci.yml | 37 ++++ src/python/pymole/__init__.py | 2 + src/python/pymole/gradient.py | 182 +++++++++++++++++++ src/python/test/test_divergence.py | 2 - src/python/test/test_gradient.py | 278 +++++++++++++++++++++++++++++ 5 files changed, 499 insertions(+), 2 deletions(-) create mode 100644 src/python/pymole/gradient.py create mode 100644 src/python/test/test_gradient.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 140c760b7..7db473486 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -271,3 +271,40 @@ jobs: path: | doc/sphinx/build/html/ doc/sphinx/build/latex/*.pdf + + Python-unit-tests: + runs-on: ubuntu-latest + + strategy: + matrix: + python-version: + - "3.10" + - "3.11" + - "3.12" + + steps: + - name: Check out repository + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + + - name: Cache pip dependencies + uses: actions/cache@v4 + with: + path: ~/.cache/pip + key: ${{ runner.os }}-pip-${{ matrix.python-version }}-${{ hashFiles('**/requirements*.txt', '**/pyproject.toml') }} + restore-keys: | + ${{ runner.os }}-pip-${{ matrix.python-version }}- + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + if [ -f requirements.txt ]; then pip install -r requirements.txt; fi + pip install pytest + + - name: Run tests + run: | + pytest \ No newline at end of file diff --git a/src/python/pymole/__init__.py b/src/python/pymole/__init__.py index 21cb17e11..14e4f9df3 100644 --- a/src/python/pymole/__init__.py +++ b/src/python/pymole/__init__.py @@ -3,12 +3,14 @@ from .dim3d import Gradient3D from .grid import Grid from .divergence import Divergence +from .gradient import Gradient __version__ = "0.1.1" __all__ = [ "Divergence", "Divergence1D", "Divergence2D", + "Gradient", "Gradient1D", "Gradient2D", "Gradient3D", diff --git a/src/python/pymole/gradient.py b/src/python/pymole/gradient.py new file mode 100644 index 000000000..aff7550f8 --- /dev/null +++ b/src/python/pymole/gradient.py @@ -0,0 +1,182 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' + +import numpy as np +from scipy import sparse + + +class Gradient: + """ + Unified gradient operator supporting 1D, 2D, and 3D Grid objects. + """ + + def __init__(self, grid, accuracy_order: int = 2): + self.grid = grid + self.accuracy_order = accuracy_order + + if self.grid.ndim not in (1, 2, 3): + raise ValueError( + f"Gradient currently supports 1D, 2D, and 3D grids. " + f"Got ndim={self.grid.ndim}" + ) + + self._matrices = self._construct_matrices() + + @classmethod + def minimum_cells_for_order(cls, accuracy_order: int) -> int: + return 2 * accuracy_order + 1 + + @property + def matrices(self): + """Returns tuple of gradient matrices (Gx,) for 1D, (Gx, Gy) for 2D, or (Gx, Gy, Gz) for 3D.""" + return self._matrices + + def _construct_matrices(self): + + if self.accuracy_order != 2: + raise ValueError( + f"Unsupported order of accuracy: {self.accuracy_order}" + ) + + match self.grid.ndim: + + case 1: + return (self._construct_1d(),) + + case 2: + Gx, Gy = self._construct_2d() + return (Gx, Gy) + + case 3: + Gx, Gy, Gz = self._construct_3d() + return (Gx, Gy, Gz) + + case _: + raise ValueError( + f"Unsupported grid dimension: {self.grid.ndim}" + ) + + def _construct_1d(self): + """Constructs 1D gradient operator (du/dx).""" + + n = len(self.grid.x) + + if n < self.minimum_cells_for_order(self.accuracy_order): + raise ValueError( + f"Grid requires at least " + f"{self.minimum_cells_for_order(self.accuracy_order)} " + f"points for order {self.accuracy_order}" + ) + + dx = self.grid.spacing + + return self._build_1d_gradient(n, dx) + + def _construct_2d(self): + """Constructs 2D gradient operators (du/dx, du/dy).""" + + m, n = self.grid.X.shape + dx, dy = self.grid.spacing + + minimum_points = self.minimum_cells_for_order( + self.accuracy_order + ) + + if m < minimum_points or n < minimum_points: + raise ValueError( + f"Both dimensions require at least " + f"{minimum_points} points." + ) + + Gx = self._build_2d_gradient_x(m, n, dx) + Gy = self._build_2d_gradient_y(m, n, dy) + + return Gx, Gy + + def _construct_3d(self): + """Constructs 3D gradient operators (du/dx, du/dy, du/dz).""" + + m, n, p = self.grid.X.shape + dx, dy, dz = self.grid.spacing + + minimum_points = self.minimum_cells_for_order( + self.accuracy_order + ) + + if m < minimum_points or n < minimum_points or p < minimum_points: + raise ValueError( + f"All dimensions require at least " + f"{minimum_points} points." + ) + + Gx = self._build_3d_gradient_x(m, n, p, dx) + Gy = self._build_3d_gradient_y(m, n, p, dy) + Gz = self._build_3d_gradient_z(m, n, p, dz) + + return Gx, Gy, Gz + + def _build_1d_gradient(self, n: int, dx: float): + """Builds 1D gradient matrix using central differences.""" + + G = sparse.lil_matrix((n, n), dtype=float) + + # Central differences for interior points: du/dx ≈ (u[i+1] - u[i-1]) / (2*dx) + for i in range(1, n - 1): + G[i, i - 1] = -1.0 / (2.0 * dx) + G[i, i + 1] = 1.0 / (2.0 * dx) + + # Forward difference at left boundary: du/dx ≈ (u[1] - u[0]) / dx + G[0, 0] = -1.0 / dx + G[0, 1] = 1.0 / dx + + # Backward difference at right boundary: du/dx ≈ (u[n-1] - u[n-2]) / dx + G[n - 1, n - 2] = -1.0 / dx + G[n - 1, n - 1] = 1.0 / dx + + return G.tocsr() + + def _build_2d_gradient_x(self, m: int, n: int, dx: float): + """Builds 2D gradient operator in x-direction.""" + + Gx_1d = self._build_1d_gradient(n, dx) + Im = sparse.eye(m, format="csr") + + return sparse.kron(Im, Gx_1d, format="csr") + + def _build_2d_gradient_y(self, m: int, n: int, dy: float): + """Builds 2D gradient operator in y-direction.""" + + Gy_1d = self._build_1d_gradient(m, dy) + In = sparse.eye(n, format="csr") + + return sparse.kron(Gy_1d, In, format="csr") + + def _build_3d_gradient_x(self, m: int, n: int, p: int, dx: float): + """Builds 3D gradient operator in x-direction.""" + + Gx_1d = self._build_1d_gradient(n, dx) + Im = sparse.eye(m, format="csr") + Ip = sparse.eye(p, format="csr") + + return sparse.kron(sparse.kron(Im, Gx_1d, format="csr"), Ip, format="csr") + + def _build_3d_gradient_y(self, m: int, n: int, p: int, dy: float): + """Builds 3D gradient operator in y-direction.""" + + Gy_1d = self._build_1d_gradient(m, dy) + In = sparse.eye(n, format="csr") + Ip = sparse.eye(p, format="csr") + + return sparse.kron(sparse.kron(Gy_1d, In, format="csr"), Ip, format="csr") + + def _build_3d_gradient_z(self, m: int, n: int, p: int, dz: float): + """Builds 3D gradient operator in z-direction.""" + + Gz_1d = self._build_1d_gradient(p, dz) + Im = sparse.eye(m, format="csr") + In = sparse.eye(n, format="csr") + + return sparse.kron(sparse.kron(Im, In, format="csr"), Gz_1d, format="csr") diff --git a/src/python/test/test_divergence.py b/src/python/test/test_divergence.py index b247745cd..273cc413e 100644 --- a/src/python/test/test_divergence.py +++ b/src/python/test/test_divergence.py @@ -10,8 +10,6 @@ from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) -#from pymole.grid import Grid -#from pymole.divergence import Divergence from pymole import Grid, Divergence class TestDivergence: diff --git a/src/python/test/test_gradient.py b/src/python/test/test_gradient.py new file mode 100644 index 000000000..b278f75dc --- /dev/null +++ b/src/python/test/test_gradient.py @@ -0,0 +1,278 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' + +import numpy as np +import pytest +import sys +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) +from pymole import Grid, Gradient + +class TestGradient: + + def test_unsupported_accuracy_order(self): + + grid = Grid.generate( + 0, + 10, + shape=10, + ) + + with pytest.raises(ValueError): + Gradient( + grid, + accuracy_order=4, + ) + + def test_grid_too_small_1d(self): + + grid = Grid.generate( + 0, + 10, + shape=4, + ) + + with pytest.raises(ValueError): + Gradient( + grid, + 2, + ) + + def test_grid_too_small_2d(self): + + grid = Grid.generate( + [0, 0], + [10, 10], + shape=[4, 4], + ) + + with pytest.raises(ValueError): + Gradient( + grid, + 2, + ) + + def test_grid_too_small_3d(self): + + grid = Grid.generate( + [0, 0, 0], + [10, 10, 10], + shape=[3, 5, 6], + ) + + with pytest.raises(ValueError): + Gradient( + grid, + 2, + ) + + def test_1d_matrix_shape(self): + + grid = Grid.generate( + 0, + 1.2, + shape=7, + ) + + Gx = Gradient( + grid, + 2, + ).matrices[0] + + n = len(grid.x) + + assert Gx.shape == ( + n, + n, + ) + + def test_2d_matrix_shape(self): + + grid = Grid.generate( + [0, 0], + [3, 2], + shape=[7, 6], + ) + + Gx, Gy = Gradient( + grid, + 2, + ).matrices + + m, n = grid.X.shape + + total_nodes = m * n + + assert Gx.shape == (total_nodes, total_nodes) + assert Gy.shape == (total_nodes, total_nodes) + + def test_3d_matrix_shape(self): + + grid = Grid.generate( + [0, 0, 0], + [3, 2, 4], + shape=[5, 6, 7], + ) + + Gx, Gy, Gz = Gradient( + grid, + 2, + ).matrices + + m, n, p = grid.X.shape + + total_nodes = m * n * p + + assert Gx.shape == (total_nodes, total_nodes) + assert Gy.shape == (total_nodes, total_nodes) + assert Gz.shape == (total_nodes, total_nodes) + + def test_1d_constant_field(self): + + grid = Grid.generate( + 0.0, + 1.2, + shape=7, + ) + + Gx = Gradient( + grid, + 2, + ).matrices[0] + + n = len(grid.x) + + u = np.ones(n) + + grad_u = Gx @ u + + assert np.allclose( + grad_u, + 0.0, + atol=1e-12, + ) + + def test_2d_constant_field(self): + + grid = Grid.generate( + [0.0, 0.0], + [3.0, 2.0], + shape=[7, 6], + ) + + Gx, Gy = Gradient( + grid, + 2, + ).matrices + + m, n = grid.X.shape + + total_nodes = m * n + + u = np.ones(total_nodes) + + grad_u_x = Gx @ u + grad_u_y = Gy @ u + + assert np.allclose(grad_u_x, 0.0, atol=1e-12) + assert np.allclose(grad_u_y, 0.0, atol=1e-12) + + def test_3d_constant_field(self): + + grid = Grid.generate( + [0.0, 0.0, 0.0], + [3.0, 2.0, 4.0], + shape=[5, 6, 7], + ) + + Gx, Gy, Gz = Gradient( + grid, + 2, + ).matrices + + m, n, p = grid.X.shape + + total_nodes = m * n * p + + u = np.ones(total_nodes) + + grad_u_x = Gx @ u + grad_u_y = Gy @ u + grad_u_z = Gz @ u + + assert np.allclose(grad_u_x, 0.0, atol=1e-12) + assert np.allclose(grad_u_y, 0.0, atol=1e-12) + assert np.allclose(grad_u_z, 0.0, atol=1e-12) + + def test_1d_linear_field(self): + + grid = Grid.generate( + 0.0, + 10.0, + shape=11, + ) + + Gx = Gradient( + grid, + 2, + ).matrices[0] + + # u = 2*x, so du/dx = 2 + u = 2 * grid.x + + grad_u = Gx @ u + + # Should be approximately 2 everywhere + assert np.allclose( + grad_u, + 2.0, + atol=1e-10, + ) + + def test_2d_linear_field_x(self): + + grid = Grid.generate( + [0.0, 0.0], + [10.0, 10.0], + shape=[11, 11], + ) + + Gx, Gy = Gradient( + grid, + 2, + ).matrices + + # u = 2*x, so du/dx = 2, du/dy = 0 + u = 2 * grid.X.flatten() + + grad_u_x = Gx @ u + grad_u_y = Gy @ u + + assert np.allclose(grad_u_x, 2.0, atol=1e-10) + assert np.allclose(grad_u_y, 0.0, atol=1e-12) + + def test_2d_linear_field_y(self): + + grid = Grid.generate( + [0.0, 0.0], + [10.0, 10.0], + shape=[11, 11], + ) + + Gx, Gy = Gradient( + grid, + 2, + ).matrices + + # u = 3*y, so du/dx = 0, du/dy = 3 + u = 3 * grid.Y.flatten() + + grad_u_x = Gx @ u + grad_u_y = Gy @ u + + assert np.allclose(grad_u_x, 0.0, atol=1e-12) + assert np.allclose(grad_u_y, 3.0, atol=1e-10) From a73f32da68c1d777b2e49352f1dc0f71698dc0d8 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 14 Jun 2026 19:35:47 -0700 Subject: [PATCH 07/66] feat: add CI workflow for automated builds, testing, and documentation, and include Python dependencies file --- .github/workflows/ci.yml | 5 +++-- src/python/pymole/requirements.txt | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7db473486..ddaed3359 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -271,7 +271,7 @@ jobs: path: | doc/sphinx/build/html/ doc/sphinx/build/latex/*.pdf - + Python-unit-tests: runs-on: ubuntu-latest @@ -303,7 +303,8 @@ jobs: run: | python -m pip install --upgrade pip if [ -f requirements.txt ]; then pip install -r requirements.txt; fi - pip install pytest + pip install pytest numpy scipy + pwd - name: Run tests run: | diff --git a/src/python/pymole/requirements.txt b/src/python/pymole/requirements.txt index 9d2ad45b9..77909df50 100644 --- a/src/python/pymole/requirements.txt +++ b/src/python/pymole/requirements.txt @@ -1,2 +1,3 @@ -numpy==2.4.0 -scipy==1.16.3 +numpy +scipy +pytest From 5c38db2e78768e044bc7e1f9c84f004b0d989f95 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 14 Jun 2026 20:01:06 -0700 Subject: [PATCH 08/66] feat: implement Grid class to manage 1D, 2D, and 3D coordinate meshes --- src/python/pymole/grid.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/python/pymole/grid.py b/src/python/pymole/grid.py index bc662827f..062f60598 100644 --- a/src/python/pymole/grid.py +++ b/src/python/pymole/grid.py @@ -70,11 +70,10 @@ def spacing(self): dx = self.x[0, 1] - self.x[0, 0] dy = self.y[1, 0] - self.y[0, 0] return dx, dy - - dx = self.x[0, 0, 1] - self.x[0, 0, 0] - dy = self.y[0, 1, 0] - self.y[0, 0, 0] - dz = self.z[1, 0, 0] - self.z[0, 0, 0] - + + dx = self.x[0, 1, 0] - self.x[0, 0, 0] + dy = self.y[1, 0, 0] - self.y[0, 0, 0] + dz = self.z[0, 0, 1] - self.z[0, 0, 0] return dx, dy, dz def __iter__(self): From 1576bf9612ddfcbe3d51e69a75f7cdb297d4e0fd Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Thu, 18 Jun 2026 15:22:25 -0700 Subject: [PATCH 09/66] feat: implement unified Laplacian operator for 1D and 2D grids with accompanying tests --- src/python/pymole/__init__.py | 2 + src/python/pymole/laplacian.py | 35 ++++++++++++ src/python/test/test_laplacian.py | 88 +++++++++++++++++++++++++++++++ 3 files changed, 125 insertions(+) create mode 100644 src/python/pymole/laplacian.py create mode 100644 src/python/test/test_laplacian.py diff --git a/src/python/pymole/__init__.py b/src/python/pymole/__init__.py index 14e4f9df3..c84fa20c7 100644 --- a/src/python/pymole/__init__.py +++ b/src/python/pymole/__init__.py @@ -4,6 +4,7 @@ from .grid import Grid from .divergence import Divergence from .gradient import Gradient +from .laplacian import Laplacian __version__ = "0.1.1" __all__ = [ @@ -16,6 +17,7 @@ "Gradient3D", "Grid", "Grid1D", + "Laplacian", "Laplacian1D", "Laplacian2D", ] diff --git a/src/python/pymole/laplacian.py b/src/python/pymole/laplacian.py new file mode 100644 index 000000000..f449489d4 --- /dev/null +++ b/src/python/pymole/laplacian.py @@ -0,0 +1,35 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' + +class Laplacian: + """ + Unified Laplacian operator for 1D and 2D Grid objects. + + The Laplacian is constructed as the matrix product of divergence and + gradient operators. + """ + + def __init__(self, grid, accuracy_order: int = 2): + self.grid = grid + self.accuracy_order = accuracy_order + + if self.grid.ndim not in (1, 2): + raise ValueError( + f"Laplacian currently supports only 1D and 2D grids. " + f"Got ndim={self.grid.ndim}" + ) + + self._matrix = self._construct_matrix() + + @property + def matrix(self): + return self._matrix + + def _construct_matrix(self): + D = Divergence(self.grid, self.accuracy_order) + G = Gradient(self.grid, self.accuracy_order) + + return D.matrix @ G.matrix diff --git a/src/python/test/test_laplacian.py b/src/python/test/test_laplacian.py new file mode 100644 index 000000000..4ba4abfdf --- /dev/null +++ b/src/python/test/test_laplacian.py @@ -0,0 +1,88 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' +import numpy as np +import pytest +import sys +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) +from pymole import Grid, Laplacian + +class TestLaplacian: + + def test_unsupported_dimension(self): + grid = Grid.generate([0, 0, 0], [10, 10, 10], shape=[5, 5, 5]) + + with pytest.raises(ValueError): + Laplacian(grid) + + def test_1d_matrix_shape(self): + grid = Grid.generate(0.0, 1.0, shape=7) + + L = Laplacian(grid).matrix + + n = len(grid.x) + + assert L.shape == (n + 1, n) + + def test_2d_matrix_shape(self): + grid = Grid.generate([0.0, 0.0], [3.0, 2.0], shape=[7, 6]) + + L = Laplacian(grid).matrix + + m, n = grid.num_cells + expected_rows = (m + 2) * (n + 2) + expected_cols = (m + 1) * (n + 1) + + assert L.shape == (expected_rows, expected_cols) + + def test_2d_constant_field_shape(self): + grid = Grid.generate([0.0, 0.0], [3.0, 2.0], shape=[7, 6]) + + L = Laplacian(grid).matrix + + m, n = grid.num_cells + total_nodes = (m + 1) * (n + 1) + + assert L.shape == ((m + 2) * (n + 2), total_nodes) + + def test_1d_constant_field(self): + grid = Grid.generate(0.0, 1.0, shape=7) + + L = Laplacian(grid).matrix + + u = np.ones(len(grid.x)) + Lu = L @ u + + assert np.allclose(Lu, 0.0, atol=1e-12) + + def test_2d_constant_field(self): + grid = Grid.generate([0.0, 0.0], [3.0, 2.0], shape=[7, 6]) + + L = Laplacian(grid).matrix + + m, n = grid.num_cells + total_nodes = (m + 1) * (n + 1) + + u = np.ones(total_nodes) + Lu = L @ u + + assert np.allclose(Lu, 0.0, atol=1e-12) + + def test_1d_quadratic_field(self): + grid = Grid.generate(0.0, 1.0, shape=7) + + L = Laplacian(grid).matrix + x = grid.x + + u = x**2 + Lu = L @ u + + # For u = x^2, d^2u/dx^2 = 2 for the interior nodes. + # The first and last rows correspond to boundary divergence values, + # so we verify the interior portion only. + assert np.allclose(Lu[2:-2], 2.0, atol=1e-1) + assert np.isfinite(Lu).all() From 87e2799c1444863282910c867c09fcdef582c6ab Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Thu, 18 Jun 2026 22:06:12 -0700 Subject: [PATCH 10/66] test: add unit tests for Laplacian matrix shape and constant field evaluation --- src/python/test/test_laplacian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/python/test/test_laplacian.py b/src/python/test/test_laplacian.py index 4ba4abfdf..25271ba0f 100644 --- a/src/python/test/test_laplacian.py +++ b/src/python/test/test_laplacian.py @@ -9,7 +9,7 @@ from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) -from pymole import Grid, Laplacian +from pymole import Grid, Gradient, Divergence, Laplacian class TestLaplacian: From 47232a86d83124b11c57c4ecd2940cc3705d18f5 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Thu, 18 Jun 2026 22:21:45 -0700 Subject: [PATCH 11/66] feat: implement Laplacian operator for 1D and 2D grids and add corresponding unit tests --- src/python/pymole/laplacian.py | 3 +++ src/python/test/test_laplacian.py | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/python/pymole/laplacian.py b/src/python/pymole/laplacian.py index f449489d4..398e3e76f 100644 --- a/src/python/pymole/laplacian.py +++ b/src/python/pymole/laplacian.py @@ -4,6 +4,9 @@ See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. ''' +from pymole import Gradient +from pymole import Divergence + class Laplacian: """ Unified Laplacian operator for 1D and 2D Grid objects. diff --git a/src/python/test/test_laplacian.py b/src/python/test/test_laplacian.py index 25271ba0f..4ba4abfdf 100644 --- a/src/python/test/test_laplacian.py +++ b/src/python/test/test_laplacian.py @@ -9,7 +9,7 @@ from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) -from pymole import Grid, Gradient, Divergence, Laplacian +from pymole import Grid, Laplacian class TestLaplacian: From c6df51bae080e839a8f678d3b6d118ecbc5847ad Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Thu, 18 Jun 2026 23:34:19 -0700 Subject: [PATCH 12/66] feat: implement Divergence and Gradient operators for multidimensional grids --- src/python/pymole/divergence.py | 1 + src/python/pymole/gradient.py | 81 ++++++++++++++++++++++++++------ src/python/test/test_gradient.py | 9 ++-- 3 files changed, 72 insertions(+), 19 deletions(-) diff --git a/src/python/pymole/divergence.py b/src/python/pymole/divergence.py index e1b054cb2..1a72b91ff 100644 --- a/src/python/pymole/divergence.py +++ b/src/python/pymole/divergence.py @@ -7,6 +7,7 @@ import numpy as np from scipy import sparse + class Divergence: """ Unified divergence operator supporting 1D and 2D Grid objects. diff --git a/src/python/pymole/gradient.py b/src/python/pymole/gradient.py index aff7550f8..83960b8a5 100644 --- a/src/python/pymole/gradient.py +++ b/src/python/pymole/gradient.py @@ -34,6 +34,23 @@ def matrices(self): """Returns tuple of gradient matrices (Gx,) for 1D, (Gx, Gy) for 2D, or (Gx, Gy, Gz) for 3D.""" return self._matrices + @property + def matrix(self): + """Returns the combined gradient matrix for divergence-based Laplacian.""" + match self.grid.ndim: + case 1: + return self._construct_1d() + case 2: + return self._construct_2d_matrix() + case 3: + raise NotImplementedError( + "Gradient.matrix is only available for 1D and 2D grids." + ) + case _: + raise ValueError( + f"Unsupported grid dimension: {self.grid.ndim}" + ) + def _construct_matrices(self): if self.accuracy_order != 2: @@ -60,9 +77,10 @@ def _construct_matrices(self): ) def _construct_1d(self): - """Constructs 1D gradient operator (du/dx).""" + """Constructs 1D gradient operator from nodal scalar values.""" n = len(self.grid.x) + dx = self.grid.spacing if n < self.minimum_cells_for_order(self.accuracy_order): raise ValueError( @@ -71,24 +89,22 @@ def _construct_1d(self): f"points for order {self.accuracy_order}" ) - dx = self.grid.spacing - return self._build_1d_gradient(n, dx) def _construct_2d(self): - """Constructs 2D gradient operators (du/dx, du/dy).""" + """Constructs 2D gradient operators from nodal scalar values.""" - m, n = self.grid.X.shape + m, n = self.grid.num_cells dx, dy = self.grid.spacing - minimum_points = self.minimum_cells_for_order( + minimum_cells = self.minimum_cells_for_order( self.accuracy_order ) - if m < minimum_points or n < minimum_points: + if m < minimum_cells or n < minimum_cells: raise ValueError( f"Both dimensions require at least " - f"{minimum_points} points." + f"{minimum_cells} cells." ) Gx = self._build_2d_gradient_x(m, n, dx) @@ -96,6 +112,12 @@ def _construct_2d(self): return Gx, Gy + def _construct_2d_matrix(self): + """Constructs the combined 2D gradient matrix [Gx; Gy].""" + + Gx, Gy = self._construct_2d() + return sparse.vstack([Gx, Gy], format="csr") + def _construct_3d(self): """Constructs 3D gradient operators (du/dx, du/dy, du/dz).""" @@ -141,18 +163,49 @@ def _build_1d_gradient(self, n: int, dx: float): def _build_2d_gradient_x(self, m: int, n: int, dx: float): """Builds 2D gradient operator in x-direction.""" - Gx_1d = self._build_1d_gradient(n, dx) - Im = sparse.eye(m, format="csr") + rows = n * (m + 1) + cols = (m + 1) * (n + 1) + G = sparse.lil_matrix((rows, cols), dtype=float) + + for yi in range(n): + for xj in range(m + 1): + row = yi * (m + 1) + xj + base = yi * (m + 1) + if xj == 0: + G[row, base + xj] = -1.0 / dx + G[row, base + xj + 1] = 1.0 / dx + elif xj == m: + G[row, base + xj - 1] = -1.0 / dx + G[row, base + xj] = 1.0 / dx + else: + G[row, base + xj - 1] = -1.0 / (2.0 * dx) + G[row, base + xj + 1] = 1.0 / (2.0 * dx) - return sparse.kron(Im, Gx_1d, format="csr") + return G.tocsr() def _build_2d_gradient_y(self, m: int, n: int, dy: float): """Builds 2D gradient operator in y-direction.""" - Gy_1d = self._build_1d_gradient(m, dy) - In = sparse.eye(n, format="csr") + rows = m * (n + 1) + cols = (m + 1) * (n + 1) + G = sparse.lil_matrix((rows, cols), dtype=float) + + for xj in range(m): + for yi in range(n + 1): + row = xj * (n + 1) + yi + col_bottom = yi * (m + 1) + xj + col_top = (yi + 1) * (m + 1) + xj + if yi == 0: + G[row, col_bottom] = -1.0 / dy + G[row, col_top] = 1.0 / dy + elif yi == n: + G[row, col_bottom - (m + 1)] = -1.0 / dy + G[row, col_bottom] = 1.0 / dy + else: + G[row, col_bottom - (m + 1)] = -1.0 / (2.0 * dy) + G[row, col_top] = 1.0 / (2.0 * dy) - return sparse.kron(Gy_1d, In, format="csr") + return G.tocsr() def _build_3d_gradient_x(self, m: int, n: int, p: int, dx: float): """Builds 3D gradient operator in x-direction.""" diff --git a/src/python/test/test_gradient.py b/src/python/test/test_gradient.py index b278f75dc..62a6e6925 100644 --- a/src/python/test/test_gradient.py +++ b/src/python/test/test_gradient.py @@ -103,12 +103,11 @@ def test_2d_matrix_shape(self): 2, ).matrices - m, n = grid.X.shape - - total_nodes = m * n + m, n = grid.num_cells + total_nodes = (m + 1) * (n + 1) - assert Gx.shape == (total_nodes, total_nodes) - assert Gy.shape == (total_nodes, total_nodes) + assert Gx.shape == (n * (m + 1), total_nodes) + assert Gy.shape == (m * (n + 1), total_nodes) def test_3d_matrix_shape(self): From adbacc0be216666f6850d7cf85a8021ca295762a Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 21 Jun 2026 22:16:37 -0700 Subject: [PATCH 13/66] feat: implement 2D Curl operator and add comprehensive Laplacian and Curl unit tests --- src/python/pymole/__init__.py | 2 + src/python/pymole/curl.py | 62 +++++++++++++++++++++++++++++++ src/python/test/test_curl.py | 32 ++++++++++++++++ src/python/test/test_laplacian.py | 11 +++++- 4 files changed, 106 insertions(+), 1 deletion(-) create mode 100644 src/python/pymole/curl.py create mode 100644 src/python/test/test_curl.py diff --git a/src/python/pymole/__init__.py b/src/python/pymole/__init__.py index c84fa20c7..007696bf1 100644 --- a/src/python/pymole/__init__.py +++ b/src/python/pymole/__init__.py @@ -2,12 +2,14 @@ from .dim2d import Divergence2D, Gradient2D, Laplacian2D from .dim3d import Gradient3D from .grid import Grid +from .curl import Curl from .divergence import Divergence from .gradient import Gradient from .laplacian import Laplacian __version__ = "0.1.1" __all__ = [ + "Curl", "Divergence", "Divergence1D", "Divergence2D", diff --git a/src/python/pymole/curl.py b/src/python/pymole/curl.py new file mode 100644 index 000000000..de38cdcaf --- /dev/null +++ b/src/python/pymole/curl.py @@ -0,0 +1,62 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' + +import numpy as np + +from pymole.divergence import Divergence + +class Curl: + """ + Unified 2D curl operator for a staggered vector field. + + The operator is implemented as the divergence of a staggered field + arranged as ``[v; -u]`` so that the discrete result corresponds to + ``dv/dx - du/dy``. + """ + + def __init__(self, grid, accuracy_order: int = 2): + self.grid = grid + self.accuracy_order = accuracy_order + + if self.grid.ndim != 2: + raise ValueError( + f"Curl currently supports only 2D grids. " + f"Got ndim={self.grid.ndim}" + ) + + self._matrix = self._construct_matrix() + + @property + def matrix(self): + return self._matrix + + def _construct_matrix(self): + if self.accuracy_order != 2: + raise ValueError( + f"Unsupported order of accuracy: {self.accuracy_order}" + ) + + return Divergence( + self.grid, + self.accuracy_order, + ).matrix + + def apply(self, u, v): + """ + Apply the discrete curl to a staggered vector field. + + Parameters + ---------- + u, v : array-like + The two components of the vector field, already arranged in the + staggered layout expected by the divergence operator. + """ + if len(u) != self._matrix.shape[1] // 2 or len(v) != self._matrix.shape[1] // 2: + raise ValueError( + "Expected u and v to each provide half of the staggered vector layout." + ) + + return self._matrix @ np.concatenate([v, -u]) diff --git a/src/python/test/test_curl.py b/src/python/test/test_curl.py new file mode 100644 index 000000000..5435de84e --- /dev/null +++ b/src/python/test/test_curl.py @@ -0,0 +1,32 @@ +import numpy as np +import pytest +import sys +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) + +from pymole import Curl, Divergence, Grid + +class TestCurl: + + def test_unsupported_dimension(self): + grid = Grid.generate(0.0, 1.0, shape=7) + + with pytest.raises(ValueError): + Curl(grid) + + def test_2d_matrix_shape(self): + grid = Grid.generate([0.0, 0.0], [3.0, 2.0], shape=[7, 6]) + + C = Curl(grid).matrix + D = Divergence(grid).matrix + + assert C.shape == D.shape + + def test_2d_constant_vector_field(self): + grid = Grid.generate([0.0, 0.0], [3.0, 2.0], shape=[7, 6]) + + C = Curl(grid).matrix + w = np.ones(C.shape[1]) + + assert np.allclose(C @ w, 0.0, atol=1e-12) diff --git a/src/python/test/test_laplacian.py b/src/python/test/test_laplacian.py index 4ba4abfdf..34c47c82d 100644 --- a/src/python/test/test_laplacian.py +++ b/src/python/test/test_laplacian.py @@ -9,7 +9,7 @@ from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) -from pymole import Grid, Laplacian +from pymole import Grid, Gradient, Divergence, Laplacian class TestLaplacian: @@ -86,3 +86,12 @@ def test_1d_quadratic_field(self): # so we verify the interior portion only. assert np.allclose(Lu[2:-2], 2.0, atol=1e-1) assert np.isfinite(Lu).all() + +def test_laplacian_equals_divergence_times_gradient(self): + grid = Grid.generate([0.0, 0.0], [3.0, 2.0], shape=[7, 6]) + + D = Divergence(grid).matrix + G = Gradient(grid).matrix + L = Laplacian(grid).matrix + + assert np.allclose((L - D @ G).toarray(), 0.0, atol=1e-12) From 0c9174f5540fc9a71c6efcc2455a2a7fe354133b Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 23 Jun 2026 21:06:54 -0700 Subject: [PATCH 14/66] test: add unit tests for Laplacian operator including matrix shapes and consistency with Divergence and Gradient --- src/python/test/test_laplacian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/python/test/test_laplacian.py b/src/python/test/test_laplacian.py index 34c47c82d..a4e64f35f 100644 --- a/src/python/test/test_laplacian.py +++ b/src/python/test/test_laplacian.py @@ -87,7 +87,7 @@ def test_1d_quadratic_field(self): assert np.allclose(Lu[2:-2], 2.0, atol=1e-1) assert np.isfinite(Lu).all() -def test_laplacian_equals_divergence_times_gradient(self): + def test_laplacian_equals_divergence_times_gradient(self): grid = Grid.generate([0.0, 0.0], [3.0, 2.0], shape=[7, 6]) D = Divergence(grid).matrix From 663b3d8e52f316db5a5eb68504e827cbc6b57cbe Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 23 Jun 2026 22:14:58 -0700 Subject: [PATCH 15/66] feat: implement base BoundaryCondition class and RobinBoundaryCondition subclass --- src/python/pymole/BoundaryCondition.py | 25 +++++++++++++++++++++ src/python/pymole/RobinBoundaryCondition.py | 11 +++++++++ 2 files changed, 36 insertions(+) create mode 100644 src/python/pymole/BoundaryCondition.py create mode 100644 src/python/pymole/RobinBoundaryCondition.py diff --git a/src/python/pymole/BoundaryCondition.py b/src/python/pymole/BoundaryCondition.py new file mode 100644 index 000000000..bfda05343 --- /dev/null +++ b/src/python/pymole/BoundaryCondition.py @@ -0,0 +1,25 @@ +from abc import ABC, abstractmethod + +class BoundaryCondition(ABC): + @abstractmethod + def apply(self) -> None: + """Subclasses must implement this.""" + pass + + def isImplemented(self) -> bool: + """Returns True if the boundary condition is implemented.""" + return self.apply != BoundaryCondition.apply + + def hasImplemented(cls, bc) -> bool: + """Returns True if the boundary condition is implemented.""" + return bc.apply != BoundaryCondition.apply + + def __init__(self, grid, accuracy_order: int = 2): + self.grid = grid + self.accuracy_order = accuracy_order + + if self.grid.ndim not in (1, 2, 3): + raise ValueError( + f"BoundaryCondition currently supports 1D, 2D, and 3D grids. " + f"Got ndim={self.grid.ndim}" + ) \ No newline at end of file diff --git a/src/python/pymole/RobinBoundaryCondition.py b/src/python/pymole/RobinBoundaryCondition.py new file mode 100644 index 000000000..f90989e45 --- /dev/null +++ b/src/python/pymole/RobinBoundaryCondition.py @@ -0,0 +1,11 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' + +from .BoundaryCondition import BoundaryCondition + +class RobinBoundaryCondition(BoundaryCondition): + def apply(self) -> None: + pass \ No newline at end of file From 5924a31acc97ad9bb2c5c4a8cf09ae6f1b3573c4 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Mon, 6 Jul 2026 16:22:55 -0700 Subject: [PATCH 16/66] feat: implement base grid, vector calculus operators, and boundary condition classes for finite difference modeling --- examples/matlab_octave/elliptic1D.m | 6 + examples/python/elliptic1D.py | 33 ++++ src/python/pymole/BoundaryCondition.py | 80 ++++++++++ src/python/pymole/RobinBoundaryCondition.py | 60 +++++++ src/python/pymole/__init__.py | 5 + src/python/pymole/curl.py | 63 ++++++++ src/python/pymole/divergence.py | 75 +++++++-- src/python/pymole/gradient.py | 167 ++++++++++++++++++-- src/python/pymole/laplacian.py | 3 + src/python/test/test_curl.py | 29 ++++ 10 files changed, 487 insertions(+), 34 deletions(-) create mode 100644 examples/python/elliptic1D.py create mode 100644 src/python/pymole/BoundaryCondition.py create mode 100644 src/python/pymole/RobinBoundaryCondition.py create mode 100644 src/python/pymole/curl.py create mode 100644 src/python/test/test_curl.py diff --git a/examples/matlab_octave/elliptic1D.m b/examples/matlab_octave/elliptic1D.m index b0ef05ac9..023b3750a 100644 --- a/examples/matlab_octave/elliptic1D.m +++ b/examples/matlab_octave/elliptic1D.m @@ -14,11 +14,17 @@ L = lap(k, m, dx); % 1D Mimetic laplacian operator +full(L) +size(L) + % Impose Robin BC on laplacian operator a = 1; b = 1; L = L + robinBC(k, m, dx, a, b); +full(L) +%size(L) + % 1D Staggered grid grid = [west west+dx/2 : dx : east-dx/2 east]; diff --git a/examples/python/elliptic1D.py b/examples/python/elliptic1D.py new file mode 100644 index 000000000..da34b53ac --- /dev/null +++ b/examples/python/elliptic1D.py @@ -0,0 +1,33 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' +import numpy as np +import sys +import shutil + +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent.parent)+'/src/python') +from pymole import Grid, Gradient, Divergence, Laplacian, RobinBoundaryCondition + +west = 0 # Domain's limits +east = 1 + +k = 6 # Operator's order of accuracy +m = 2*k+1 # Minimum number of cells to attain the desired accuracy + +grid = Grid.generate(west, east, shape=m+1) + +L = Laplacian(grid, accuracy_order=k).matrix + +terminal_width, _ = shutil.get_terminal_size() +np.set_printoptions(linewidth=terminal_width, suppress=True) +print(L.toarray()) + +BC = RobinBoundaryCondition(grid, k, dirichlet_coefficient=1.0, neumann_coefficient=1.0) +print(BC.matrix.toarray()) + +L_BC = L + BC.matrix +print(L_BC.toarray()) \ No newline at end of file diff --git a/src/python/pymole/BoundaryCondition.py b/src/python/pymole/BoundaryCondition.py new file mode 100644 index 000000000..3042c3d99 --- /dev/null +++ b/src/python/pymole/BoundaryCondition.py @@ -0,0 +1,80 @@ +from abc import ABC, abstractmethod + +import numpy as np +from scipy import sparse + + +class BoundaryCondition(ABC): + @abstractmethod + def apply(self): + """Subclasses must implement this.""" + raise NotImplementedError + + def isImplemented(self) -> bool: + """Returns True if the boundary condition is implemented.""" + return self.apply != BoundaryCondition.apply + + @classmethod + def hasImplemented(cls, bc) -> bool: + """Returns True if the boundary condition is implemented.""" + return bc.apply != BoundaryCondition.apply + + def __init__(self, grid, accuracy_order: int = 2): + self.grid = grid + self.accuracy_order = accuracy_order + + if self.grid.ndim not in (1, 2, 3): + raise ValueError( + f"BoundaryCondition currently supports 1D, 2D, and 3D grids. " + f"Got ndim={self.grid.ndim}" + ) + + self._configure_dimensions() + self._matrix = None + + def _configure_dimensions(self) -> None: + if self.grid.ndim == 1: + self.m = self.grid.num_cells + self.n = None + self.o = None + elif self.grid.ndim == 2: + self.m, self.n = self.grid.num_cells + self.o = None + else: + self.m, self.n, self.o = self.grid.num_cells + + @property + def matrix(self): + if self._matrix is None: + self._matrix = self.apply() + return self._matrix + + def _spacing_for_dimension(self, axis: int) -> float: + spacing = self.grid.spacing + if self.grid.ndim == 1: + return spacing + if self.grid.ndim == 2: + return spacing[axis] + return spacing[axis] + + def _build_1d_robin_matrix(self, size: int, dx: float, a: float, b: float): + A = sparse.lil_matrix((size + 2, size + 2), dtype=float) + A[0, 0] = a + A[-1, -1] = a + + B = sparse.lil_matrix((size + 2, size + 1), dtype=float) + B[0, 0] = -b + B[-1, -1] = b + + G = self._build_1d_gradient(size, dx) + return (A + B @ G).tocsr() + + def _build_1d_gradient(self, size: int, dx: float): + G = sparse.lil_matrix((size + 1, size + 2), dtype=float) + G[0, 0:3] = [-8.0 / 3.0, 3.0, -1.0 / 3.0] + G[size, size - 1 : size + 2] = [1.0 / 3.0, -3.0, 8.0 / 3.0] + + for i in range(1, size): + G[i, i : i + 2] = [-1.0, 1.0] + + return (G / dx).tocsr() diff --git a/src/python/pymole/RobinBoundaryCondition.py b/src/python/pymole/RobinBoundaryCondition.py new file mode 100644 index 000000000..c13232641 --- /dev/null +++ b/src/python/pymole/RobinBoundaryCondition.py @@ -0,0 +1,60 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' + +from scipy import sparse + +from .BoundaryCondition import BoundaryCondition + + +class RobinBoundaryCondition(BoundaryCondition): + def __init__(self, grid, accuracy_order: int = 2, dirichlet_coefficient: float = 1.0, neumann_coefficient: float = 1.0): + self.dirichlet_coefficient = dirichlet_coefficient + self.neumann_coefficient = neumann_coefficient + super().__init__(grid, accuracy_order) + + def apply(self): + if self.grid.ndim == 1: + return self._build_1d_robin_matrix(self.m, self._spacing_for_dimension(0), self.dirichlet_coefficient, self.neumann_coefficient) + + if self.grid.ndim == 2: + Bm = self._build_1d_robin_matrix(self.m, self._spacing_for_dimension(0), self.dirichlet_coefficient, self.neumann_coefficient) + Bn = self._build_1d_robin_matrix(self.n, self._spacing_for_dimension(1), self.dirichlet_coefficient, self.neumann_coefficient) + + Im = sparse.eye(self.m + 2, format="csr") + In = sparse.eye(self.n + 2, format="csr") + In = In.tolil() + In[0, 0] = 0.0 + In[-1, -1] = 0.0 + In = In.tocsr() + + BC1 = sparse.kron(In, Bm, format="csr") + BC2 = sparse.kron(Bn, Im, format="csr") + return (BC1 + BC2).tocsr() + + if self.grid.ndim == 3: + Bm = self._build_1d_robin_matrix(self.m, self._spacing_for_dimension(0), self.dirichlet_coefficient, self.neumann_coefficient) + Bn = self._build_1d_robin_matrix(self.n, self._spacing_for_dimension(1), self.dirichlet_coefficient, self.neumann_coefficient) + Bo = self._build_1d_robin_matrix(self.o, self._spacing_for_dimension(2), self.dirichlet_coefficient, self.neumann_coefficient) + + Im = sparse.eye(self.m + 2, format="csr") + In = sparse.eye(self.n + 2, format="csr") + Io = sparse.eye(self.o + 2, format="csr") + Io = Io.tolil() + Io[0, 0] = 0.0 + Io[-1, -1] = 0.0 + Io = Io.tocsr() + + In2 = In.tolil() + In2[0, 0] = 0.0 + In2[-1, -1] = 0.0 + In2 = In2.tocsr() + + BC1 = sparse.kron(sparse.kron(Io, In2, format="csr"), Bm, format="csr") + BC2 = sparse.kron(sparse.kron(Io, Bn, format="csr"), Im, format="csr") + BC3 = sparse.kron(sparse.kron(Bo, In, format="csr"), Im, format="csr") + return (BC1 + BC2 + BC3).tocsr() + + raise ValueError(f"Unsupported grid dimension: {self.grid.ndim}") \ No newline at end of file diff --git a/src/python/pymole/__init__.py b/src/python/pymole/__init__.py index c84fa20c7..16ea616ea 100644 --- a/src/python/pymole/__init__.py +++ b/src/python/pymole/__init__.py @@ -5,9 +5,12 @@ from .divergence import Divergence from .gradient import Gradient from .laplacian import Laplacian +from .BoundaryCondition import BoundaryCondition +from .RobinBoundaryCondition import RobinBoundaryCondition __version__ = "0.1.1" __all__ = [ + "BoundaryCondition", "Divergence", "Divergence1D", "Divergence2D", @@ -20,4 +23,6 @@ "Laplacian", "Laplacian1D", "Laplacian2D", + "RobinBoundaryCondition", ] + diff --git a/src/python/pymole/curl.py b/src/python/pymole/curl.py new file mode 100644 index 000000000..fff4d8c44 --- /dev/null +++ b/src/python/pymole/curl.py @@ -0,0 +1,63 @@ +''' + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +''' + +import numpy as np + +from pymole.divergence import Divergence + + +class Curl: + """ + Unified 2D curl operator for a staggered vector field. + + The operator is implemented as the divergence of a staggered field + arranged as ``[v; -u]`` so that the discrete result corresponds to + ``dv/dx - du/dy``. + """ + + def __init__(self, grid, accuracy_order: int = 2): + self.grid = grid + self.accuracy_order = accuracy_order + + if self.grid.ndim != 2: + raise ValueError( + f"Curl currently supports only 2D grids. " + f"Got ndim={self.grid.ndim}" + ) + + self._matrix = self._construct_matrix() + + @property + def matrix(self): + return self._matrix + + def _construct_matrix(self): + if self.accuracy_order != 2: + raise ValueError( + f"Unsupported order of accuracy: {self.accuracy_order}" + ) + + return Divergence( + self.grid, + self.accuracy_order, + ).matrix + + def apply(self, u, v): + """ + Apply the discrete curl to a staggered vector field. + + Parameters + ---------- + u, v : array-like + The two components of the vector field, already arranged in the + staggered layout expected by the divergence operator. + """ + if len(u) != self._matrix.shape[1] // 2 or len(v) != self._matrix.shape[1] // 2: + raise ValueError( + "Expected u and v to each provide half of the staggered vector layout." + ) + + return self._matrix @ np.concatenate([v, -u]) diff --git a/src/python/pymole/divergence.py b/src/python/pymole/divergence.py index e1b054cb2..939c8b9be 100644 --- a/src/python/pymole/divergence.py +++ b/src/python/pymole/divergence.py @@ -7,6 +7,7 @@ import numpy as np from scipy import sparse + class Divergence: """ Unified divergence operator supporting 1D and 2D Grid objects. @@ -33,8 +34,7 @@ def matrix(self): return self._matrix def _construct_matrix(self): - - if self.accuracy_order != 2: + if self.accuracy_order not in (2, 4, 6): raise ValueError( f"Unsupported order of accuracy: {self.accuracy_order}" ) @@ -45,6 +45,10 @@ def _construct_matrix(self): return self._construct_1d() case 2: + if self.accuracy_order != 2: + raise ValueError( + f"Unsupported order of accuracy for 2D Divergence: {self.accuracy_order}" + ) return self._construct_2d() case _: @@ -56,11 +60,15 @@ def _construct_1d(self): m = self.grid.num_cells - if m < self.minimum_cells_for_order(self.accuracy_order): + # In this Grid API, users often specify the number of points via shape. + # Validate against points so k=2,4,6 examples match MATLAB-style usage. + num_points = len(self.grid.x) + + if num_points < self.minimum_cells_for_order(self.accuracy_order): raise ValueError( f"Grid requires at least " f"{self.minimum_cells_for_order(self.accuracy_order)} " - f"cells for order {self.accuracy_order}" + f"points for order {self.accuracy_order}" ) dx = self.grid.spacing @@ -95,23 +103,56 @@ def _build_1d_divergence( dx: float, ): - main_diag = np.ones(m + 1) + D = sparse.lil_matrix((m + 2, m + 1), dtype=float) - sub_diag = np.full( - m + 1, - -1.0, - ) + match self.accuracy_order: + case 2: + for i in range(1, m + 1): + D[i, i - 1 : i + 1] = [-1.0, 1.0] + + case 4: + A = np.array([-11.0 / 12.0, 17.0 / 24.0, 3.0 / 8.0, -5.0 / 24.0, 1.0 / 24.0]) + D[1, 0:5] = A + D[m, m - 4 : m + 1] = -A[::-1] + + for i in range(2, m): + D[i, i - 2 : i + 2] = [1.0 / 24.0, -9.0 / 8.0, 9.0 / 8.0, -1.0 / 24.0] + + case 6: + A = np.array( + [ + [ + -1627.0 / 1920.0, + 211.0 / 640.0, + 59.0 / 48.0, + -235.0 / 192.0, + 91.0 / 128.0, + -443.0 / 1920.0, + 31.0 / 960.0, + ], + [ + 31.0 / 960.0, + -687.0 / 640.0, + 129.0 / 128.0, + 19.0 / 192.0, + -3.0 / 32.0, + 21.0 / 640.0, + -3.0 / 640.0, + ], + ] + ) + D[1:3, 0:7] = A + D[m - 1 : m + 1, m - 6 : m + 1] = -np.rot90(A, 2) - main_diag[0] = 0.0 - sub_diag[-1] = 0.0 + for i in range(3, m - 1): + D[i, i - 3 : i + 3] = [-3.0 / 640.0, 25.0 / 384.0, -75.0 / 64.0, 75.0 / 64.0, -25.0 / 384.0, 3.0 / 640.0] - D = sparse.diags( - [sub_diag, main_diag], - offsets=[-1, 0], - shape=(m + 2, m + 1), - ).tocsr() + case _: + raise ValueError( + f"Unsupported order of accuracy: {self.accuracy_order}" + ) - return D / dx + return (D / dx).tocsr() def _build_2d_divergence( self, diff --git a/src/python/pymole/gradient.py b/src/python/pymole/gradient.py index aff7550f8..5d602a3b4 100644 --- a/src/python/pymole/gradient.py +++ b/src/python/pymole/gradient.py @@ -34,9 +34,25 @@ def matrices(self): """Returns tuple of gradient matrices (Gx,) for 1D, (Gx, Gy) for 2D, or (Gx, Gy, Gz) for 3D.""" return self._matrices - def _construct_matrices(self): + @property + def matrix(self): + """Returns the combined gradient matrix for divergence-based Laplacian.""" + match self.grid.ndim: + case 1: + return self._construct_1d() + case 2: + return self._construct_2d_matrix() + case 3: + raise NotImplementedError( + "Gradient.matrix is only available for 1D and 2D grids." + ) + case _: + raise ValueError( + f"Unsupported grid dimension: {self.grid.ndim}" + ) - if self.accuracy_order != 2: + def _construct_matrices(self): + if self.accuracy_order not in (2, 4, 6): raise ValueError( f"Unsupported order of accuracy: {self.accuracy_order}" ) @@ -47,10 +63,18 @@ def _construct_matrices(self): return (self._construct_1d(),) case 2: + if self.accuracy_order != 2: + raise ValueError( + f"Unsupported order of accuracy for 2D Gradient: {self.accuracy_order}" + ) Gx, Gy = self._construct_2d() return (Gx, Gy) case 3: + if self.accuracy_order != 2: + raise ValueError( + f"Unsupported order of accuracy for 3D Gradient: {self.accuracy_order}" + ) Gx, Gy, Gz = self._construct_3d() return (Gx, Gy, Gz) @@ -60,9 +84,10 @@ def _construct_matrices(self): ) def _construct_1d(self): - """Constructs 1D gradient operator (du/dx).""" + """Constructs 1D gradient operator from nodal scalar values.""" n = len(self.grid.x) + dx = self.grid.spacing if n < self.minimum_cells_for_order(self.accuracy_order): raise ValueError( @@ -71,24 +96,95 @@ def _construct_1d(self): f"points for order {self.accuracy_order}" ) - dx = self.grid.spacing + m = self.grid.num_cells + return self._build_1d_gradient_mimetic(m, dx) - return self._build_1d_gradient(n, dx) + def _build_1d_gradient_mimetic(self, m: int, dx: float): + """Builds 1D non-periodic mimetic gradient (MATLAB gradNonPeriodic).""" + + G = sparse.lil_matrix((m + 1, m + 2), dtype=float) + + match self.accuracy_order: + case 2: + A = np.array([-8.0 / 3.0, 3.0, -1.0 / 3.0]) + G[0, 0:3] = A + G[m, m - 1 : m + 2] = -A[::-1] + + for i in range(1, m): + G[i, i : i + 2] = [-1.0, 1.0] + + case 4: + A = np.array( + [ + [-352.0 / 105.0, 35.0 / 8.0, -35.0 / 24.0, 21.0 / 40.0, -5.0 / 56.0], + [16.0 / 105.0, -31.0 / 24.0, 29.0 / 24.0, -3.0 / 40.0, 1.0 / 168.0], + ] + ) + G[0:2, 0:5] = A + G[m - 1 : m + 1, m - 3 : m + 2] = -np.rot90(A, 2) + + for i in range(2, m - 1): + G[i, i - 1 : i + 3] = [1.0 / 24.0, -9.0 / 8.0, 9.0 / 8.0, -1.0 / 24.0] + + case 6: + A = np.array( + [ + [ + -13016.0 / 3465.0, + 693.0 / 128.0, + -385.0 / 128.0, + 693.0 / 320.0, + -495.0 / 448.0, + 385.0 / 1152.0, + -63.0 / 1408.0, + ], + [ + 496.0 / 3465.0, + -811.0 / 640.0, + 449.0 / 384.0, + -29.0 / 960.0, + -11.0 / 448.0, + 13.0 / 1152.0, + -37.0 / 21120.0, + ], + [ + -8.0 / 385.0, + 179.0 / 1920.0, + -153.0 / 128.0, + 381.0 / 320.0, + -101.0 / 1344.0, + 1.0 / 128.0, + -3.0 / 7040.0, + ], + ] + ) + G[0:3, 0:7] = A + G[m - 2 : m + 1, m - 5 : m + 2] = -np.rot90(A, 2) + + for i in range(3, m - 2): + G[i, i - 2 : i + 4] = [-3.0 / 640.0, 25.0 / 384.0, -75.0 / 64.0, 75.0 / 64.0, -25.0 / 384.0, 3.0 / 640.0] + + case _: + raise ValueError( + f"Unsupported order of accuracy: {self.accuracy_order}" + ) + + return (G / dx).tocsr() def _construct_2d(self): - """Constructs 2D gradient operators (du/dx, du/dy).""" + """Constructs 2D gradient operators from nodal scalar values.""" - m, n = self.grid.X.shape + m, n = self.grid.num_cells dx, dy = self.grid.spacing - minimum_points = self.minimum_cells_for_order( + minimum_cells = self.minimum_cells_for_order( self.accuracy_order ) - if m < minimum_points or n < minimum_points: + if m < minimum_cells or n < minimum_cells: raise ValueError( f"Both dimensions require at least " - f"{minimum_points} points." + f"{minimum_cells} cells." ) Gx = self._build_2d_gradient_x(m, n, dx) @@ -96,6 +192,12 @@ def _construct_2d(self): return Gx, Gy + def _construct_2d_matrix(self): + """Constructs the combined 2D gradient matrix [Gx; Gy].""" + + Gx, Gy = self._construct_2d() + return sparse.vstack([Gx, Gy], format="csr") + def _construct_3d(self): """Constructs 3D gradient operators (du/dx, du/dy, du/dz).""" @@ -119,7 +221,7 @@ def _construct_3d(self): return Gx, Gy, Gz def _build_1d_gradient(self, n: int, dx: float): - """Builds 1D gradient matrix using central differences.""" + """Builds square 1D central-difference gradient used by 3D helpers.""" G = sparse.lil_matrix((n, n), dtype=float) @@ -141,18 +243,49 @@ def _build_1d_gradient(self, n: int, dx: float): def _build_2d_gradient_x(self, m: int, n: int, dx: float): """Builds 2D gradient operator in x-direction.""" - Gx_1d = self._build_1d_gradient(n, dx) - Im = sparse.eye(m, format="csr") + rows = n * (m + 1) + cols = (m + 1) * (n + 1) + G = sparse.lil_matrix((rows, cols), dtype=float) + + for yi in range(n): + for xj in range(m + 1): + row = yi * (m + 1) + xj + base = yi * (m + 1) + if xj == 0: + G[row, base + xj] = -1.0 / dx + G[row, base + xj + 1] = 1.0 / dx + elif xj == m: + G[row, base + xj - 1] = -1.0 / dx + G[row, base + xj] = 1.0 / dx + else: + G[row, base + xj - 1] = -1.0 / (2.0 * dx) + G[row, base + xj + 1] = 1.0 / (2.0 * dx) - return sparse.kron(Im, Gx_1d, format="csr") + return G.tocsr() def _build_2d_gradient_y(self, m: int, n: int, dy: float): """Builds 2D gradient operator in y-direction.""" - Gy_1d = self._build_1d_gradient(m, dy) - In = sparse.eye(n, format="csr") + rows = m * (n + 1) + cols = (m + 1) * (n + 1) + G = sparse.lil_matrix((rows, cols), dtype=float) + + for xj in range(m): + for yi in range(n + 1): + row = xj * (n + 1) + yi + col_bottom = yi * (m + 1) + xj + col_top = (yi + 1) * (m + 1) + xj + if yi == 0: + G[row, col_bottom] = -1.0 / dy + G[row, col_top] = 1.0 / dy + elif yi == n: + G[row, col_bottom - (m + 1)] = -1.0 / dy + G[row, col_bottom] = 1.0 / dy + else: + G[row, col_bottom - (m + 1)] = -1.0 / (2.0 * dy) + G[row, col_top] = 1.0 / (2.0 * dy) - return sparse.kron(Gy_1d, In, format="csr") + return G.tocsr() def _build_3d_gradient_x(self, m: int, n: int, p: int, dx: float): """Builds 3D gradient operator in x-direction.""" diff --git a/src/python/pymole/laplacian.py b/src/python/pymole/laplacian.py index f449489d4..c7a14aad7 100644 --- a/src/python/pymole/laplacian.py +++ b/src/python/pymole/laplacian.py @@ -3,6 +3,9 @@ © 2008-2024 San Diego State University Research Foundation (SDSURF). See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. ''' +import numpy as np +from pymole.divergence import Divergence +from pymole.gradient import Gradient class Laplacian: """ diff --git a/src/python/test/test_curl.py b/src/python/test/test_curl.py new file mode 100644 index 000000000..a11becb51 --- /dev/null +++ b/src/python/test/test_curl.py @@ -0,0 +1,29 @@ +import numpy as np +import pytest + +from pymole import Curl, Divergence, Grid + + +class TestCurl: + + def test_unsupported_dimension(self): + grid = Grid.generate(0.0, 1.0, shape=7) + + with pytest.raises(ValueError): + Curl(grid) + + def test_2d_matrix_shape(self): + grid = Grid.generate([0.0, 0.0], [3.0, 2.0], shape=[7, 6]) + + C = Curl(grid).matrix + D = Divergence(grid).matrix + + assert C.shape == D.shape + + def test_2d_constant_vector_field(self): + grid = Grid.generate([0.0, 0.0], [3.0, 2.0], shape=[7, 6]) + + C = Curl(grid).matrix + w = np.ones(C.shape[1]) + + assert np.allclose(C @ w, 0.0, atol=1e-12) From 88ab978d5cac20b8e715eb7ad9e14b2c34b344e8 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 7 Jul 2026 23:23:30 -0700 Subject: [PATCH 17/66] feat: add BoundaryCondition base class and implement 1D Robin boundary condition example --- examples/python/elliptic1D.py | 19 +++++++++++++++---- src/python/pymole/BoundaryCondition.py | 17 ++++------------- 2 files changed, 19 insertions(+), 17 deletions(-) diff --git a/examples/python/elliptic1D.py b/examples/python/elliptic1D.py index da34b53ac..8a70a9f99 100644 --- a/examples/python/elliptic1D.py +++ b/examples/python/elliptic1D.py @@ -6,6 +6,7 @@ import numpy as np import sys import shutil +import matplotlib.pyplot as plt from pathlib import Path @@ -24,10 +25,20 @@ terminal_width, _ = shutil.get_terminal_size() np.set_printoptions(linewidth=terminal_width, suppress=True) -print(L.toarray()) BC = RobinBoundaryCondition(grid, k, dirichlet_coefficient=1.0, neumann_coefficient=1.0) -print(BC.matrix.toarray()) - L_BC = L + BC.matrix -print(L_BC.toarray()) \ No newline at end of file + +x = np.r_[west, np.arange(west + grid.spacing/2, east - grid.spacing/2 + grid.spacing/2, grid.spacing), east] +U = np.exp(x)[:, None] +U[0] = 0 # West BC +U[-1] = 2*np.exp(1) # East BC +U = np.linalg.inv(L_BC.toarray()) @ U + +plt.plot(x, U, 'o', label='Approximated') +plt.plot(x, np.exp(x), label='Analytical') +plt.legend(loc='upper left') +plt.title("Poisson's equation with Robin BC") +plt.xlabel('x') +plt.ylabel('u(x)') +plt.show() \ No newline at end of file diff --git a/src/python/pymole/BoundaryCondition.py b/src/python/pymole/BoundaryCondition.py index 3042c3d99..073d94ad4 100644 --- a/src/python/pymole/BoundaryCondition.py +++ b/src/python/pymole/BoundaryCondition.py @@ -3,6 +3,7 @@ import numpy as np from scipy import sparse +from .gradient import Gradient class BoundaryCondition(ABC): @abstractmethod @@ -65,16 +66,6 @@ def _build_1d_robin_matrix(self, size: int, dx: float, a: float, b: float): B = sparse.lil_matrix((size + 2, size + 1), dtype=float) B[0, 0] = -b B[-1, -1] = b - - G = self._build_1d_gradient(size, dx) - return (A + B @ G).tocsr() - - def _build_1d_gradient(self, size: int, dx: float): - G = sparse.lil_matrix((size + 1, size + 2), dtype=float) - G[0, 0:3] = [-8.0 / 3.0, 3.0, -1.0 / 3.0] - G[size, size - 1 : size + 2] = [1.0 / 3.0, -3.0, 8.0 / 3.0] - - for i in range(1, size): - G[i, i : i + 2] = [-1.0, 1.0] - - return (G / dx).tocsr() + + G = Gradient(self.grid, self.accuracy_order).matrix + return (A + B @ G).tocsr() \ No newline at end of file From b4ef90cd4b85ac85f63f82fa7533e0010f45e389 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 8 Jul 2026 08:37:11 -0700 Subject: [PATCH 18/66] test: add unit tests for Curl operator functionality and dimension validation --- src/python/test/test_curl.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/python/test/test_curl.py b/src/python/test/test_curl.py index eb2fd80c3..b58642629 100644 --- a/src/python/test/test_curl.py +++ b/src/python/test/test_curl.py @@ -1,10 +1,11 @@ import numpy as np import pytest +import sys +from pathlib import Path +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) from pymole import Curl, Divergence, Grid - - class TestCurl: def test_unsupported_dimension(self): From 605e04d5ebf73c336a9c933d27cb842e9523e991 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 22 Jul 2026 23:20:46 -0500 Subject: [PATCH 19/66] test: add unit tests for gradient operator functionality and shape validation --- src/python/test/test_gradient.py | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/src/python/test/test_gradient.py b/src/python/test/test_gradient.py index 62a6e6925..9885cf41f 100644 --- a/src/python/test/test_gradient.py +++ b/src/python/test/test_gradient.py @@ -25,7 +25,7 @@ def test_unsupported_accuracy_order(self): with pytest.raises(ValueError): Gradient( grid, - accuracy_order=4, + accuracy_order=8, ) def test_grid_too_small_1d(self): @@ -87,7 +87,7 @@ def test_1d_matrix_shape(self): assert Gx.shape == ( n, - n, + n+1, ) def test_2d_matrix_shape(self): @@ -141,11 +141,11 @@ def test_1d_constant_field(self): Gx = Gradient( grid, 2, - ).matrices[0] + ).matrix n = len(grid.x) - u = np.ones(n) + u = np.ones(n+1) grad_u = Gx @ u @@ -207,29 +207,33 @@ def test_3d_constant_field(self): assert np.allclose(grad_u_y, 0.0, atol=1e-12) assert np.allclose(grad_u_z, 0.0, atol=1e-12) - def test_1d_linear_field(self): + def test_1d_gradient_accuracy(self): grid = Grid.generate( 0.0, 10.0, - shape=11, + shape=101, ) + + x = np.linspace(0, 10, 101) + xSq = x ** 2 + dx = x[1] - x[0] + expected_grad_xSq = np.gradient(xSq, dx) + dx = grid.x[1] - grid.x[0] Gx = Gradient( grid, 2, - ).matrices[0] + ).matrix - # u = 2*x, so du/dx = 2 - u = 2 * grid.x + xSq = np.append(xSq,xSq[-1]+dx) - grad_u = Gx @ u + grad_xSq = Gx @ xSq - # Should be approximately 2 everywhere assert np.allclose( - grad_u, - 2.0, - atol=1e-10, + grad_xSq[:-1], + expected_grad_xSq[:-1], + atol=1e-1, ) def test_2d_linear_field_x(self): From f5bfcc95ded071d551f1b5582767b12cfbfb2a30 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 4 Aug 2026 22:23:45 -0700 Subject: [PATCH 20/66] test: add unit tests for Gradient, Divergence, and Laplacian operators --- src/python/test/test_divergence.py | 4 ++-- src/python/test/test_gradient.py | 6 ++++-- src/python/test/test_laplacian.py | 7 ++----- 3 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/python/test/test_divergence.py b/src/python/test/test_divergence.py index 273cc413e..8aa5db2bf 100644 --- a/src/python/test/test_divergence.py +++ b/src/python/test/test_divergence.py @@ -25,7 +25,7 @@ def test_unsupported_accuracy_order(self): with pytest.raises(ValueError): Divergence( grid, - accuracy_order=4, + accuracy_order=6, ) def test_grid_too_small_1d(self): @@ -39,7 +39,7 @@ def test_grid_too_small_1d(self): with pytest.raises(ValueError): Divergence( grid, - 2, + 4, ) def test_grid_too_small_2d(self): diff --git a/src/python/test/test_gradient.py b/src/python/test/test_gradient.py index 9885cf41f..b81544295 100644 --- a/src/python/test/test_gradient.py +++ b/src/python/test/test_gradient.py @@ -219,7 +219,6 @@ def test_1d_gradient_accuracy(self): xSq = x ** 2 dx = x[1] - x[0] expected_grad_xSq = np.gradient(xSq, dx) - dx = grid.x[1] - grid.x[0] Gx = Gradient( grid, @@ -227,9 +226,12 @@ def test_1d_gradient_accuracy(self): ).matrix xSq = np.append(xSq,xSq[-1]+dx) - + grad_xSq = Gx @ xSq + print(grad_xSq[:-1]) + print(expected_grad_xSq[:-1]) + assert np.allclose( grad_xSq[:-1], expected_grad_xSq[:-1], diff --git a/src/python/test/test_laplacian.py b/src/python/test/test_laplacian.py index a4e64f35f..882d74a80 100644 --- a/src/python/test/test_laplacian.py +++ b/src/python/test/test_laplacian.py @@ -26,7 +26,7 @@ def test_1d_matrix_shape(self): n = len(grid.x) - assert L.shape == (n + 1, n) + assert L.shape == (n + 1, n + 1) def test_2d_matrix_shape(self): grid = Grid.generate([0.0, 0.0], [3.0, 2.0], shape=[7, 6]) @@ -51,12 +51,9 @@ def test_2d_constant_field_shape(self): def test_1d_constant_field(self): grid = Grid.generate(0.0, 1.0, shape=7) - L = Laplacian(grid).matrix - - u = np.ones(len(grid.x)) + u = np.ones(len(grid.x) + 1) Lu = L @ u - assert np.allclose(Lu, 0.0, atol=1e-12) def test_2d_constant_field(self): From 843c339c30cab49a98b6bfc347b467009d6f423e Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 4 Aug 2026 22:49:39 -0700 Subject: [PATCH 21/66] test: add unit tests for Laplacian operator matrix construction and properties --- src/python/test/test_laplacian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/python/test/test_laplacian.py b/src/python/test/test_laplacian.py index 882d74a80..d86977a52 100644 --- a/src/python/test/test_laplacian.py +++ b/src/python/test/test_laplacian.py @@ -74,7 +74,7 @@ def test_1d_quadratic_field(self): L = Laplacian(grid).matrix x = grid.x - + x = np.append(x, x[-1] + grid.spacing) u = x**2 Lu = L @ u From fd4a2273141b115608ca54ebc9bf594a4224206c Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 11 Aug 2026 19:59:56 -0700 Subject: [PATCH 22/66] Add new directory --- python/.gitkeep | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 python/.gitkeep diff --git a/python/.gitkeep b/python/.gitkeep new file mode 100644 index 000000000..e69de29bb From 7440a91d19b103329d39c5f82419e49e2a27d579 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 11 Aug 2026 20:21:26 -0700 Subject: [PATCH 23/66] Add new hierarchy --- python/LICENSE | 0 python/README.md | 0 python/docs/.gitkeep | 0 python/examples/.gitkeep | 0 python/pyproject.toml | 0 python/src/.gitkeep | 0 python/tests/.gitkeep | 0 7 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 python/LICENSE create mode 100644 python/README.md create mode 100644 python/docs/.gitkeep create mode 100644 python/examples/.gitkeep create mode 100644 python/pyproject.toml create mode 100644 python/src/.gitkeep create mode 100644 python/tests/.gitkeep diff --git a/python/LICENSE b/python/LICENSE new file mode 100644 index 000000000..e69de29bb diff --git a/python/README.md b/python/README.md new file mode 100644 index 000000000..e69de29bb diff --git a/python/docs/.gitkeep b/python/docs/.gitkeep new file mode 100644 index 000000000..e69de29bb diff --git a/python/examples/.gitkeep b/python/examples/.gitkeep new file mode 100644 index 000000000..e69de29bb diff --git a/python/pyproject.toml b/python/pyproject.toml new file mode 100644 index 000000000..e69de29bb diff --git a/python/src/.gitkeep b/python/src/.gitkeep new file mode 100644 index 000000000..e69de29bb diff --git a/python/tests/.gitkeep b/python/tests/.gitkeep new file mode 100644 index 000000000..e69de29bb From b29252b5365fa85691ed0c1a01287ef9114eaf18 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 11 Aug 2026 20:33:59 -0700 Subject: [PATCH 24/66] Move files into new hierarchy --- {examples/python => python/examples}/elliptic1D.py | 0 python/pyproject.toml | 9 +++++++++ {src/python => python/src}/pymole/BoundaryCondition.py | 0 .../src}/pymole/RobinBoundaryCondition.py | 0 {src/python => python/src}/pymole/__init__.py | 0 {src/python => python/src}/pymole/curl.py | 0 {src/python => python/src}/pymole/divergence.py | 0 {src/python => python/src}/pymole/gradient.py | 0 {src/python => python/src}/pymole/grid.py | 0 {src/python => python/src}/pymole/laplacian.py | 0 {src/python => python/src}/pymole/mimetic_operator.py | 0 {src/python => python/src}/pymole/requirements.txt | 0 {src/python/test => python/tests}/__init__.py | 0 {src/python/test => python/tests}/test_curl.py | 0 {src/python/test => python/tests}/test_divergence.py | 0 {src/python/test => python/tests}/test_gradient.py | 0 {src/python/test => python/tests}/test_grid_utility.py | 0 {src/python/test => python/tests}/test_laplacian.py | 0 src/python/pyproject.toml | 9 --------- 19 files changed, 9 insertions(+), 9 deletions(-) rename {examples/python => python/examples}/elliptic1D.py (100%) rename {src/python => python/src}/pymole/BoundaryCondition.py (100%) rename {src/python => python/src}/pymole/RobinBoundaryCondition.py (100%) rename {src/python => python/src}/pymole/__init__.py (100%) rename {src/python => python/src}/pymole/curl.py (100%) rename {src/python => python/src}/pymole/divergence.py (100%) rename {src/python => python/src}/pymole/gradient.py (100%) rename {src/python => python/src}/pymole/grid.py (100%) rename {src/python => python/src}/pymole/laplacian.py (100%) rename {src/python => python/src}/pymole/mimetic_operator.py (100%) rename {src/python => python/src}/pymole/requirements.txt (100%) rename {src/python/test => python/tests}/__init__.py (100%) rename {src/python/test => python/tests}/test_curl.py (100%) rename {src/python/test => python/tests}/test_divergence.py (100%) rename {src/python/test => python/tests}/test_gradient.py (100%) rename {src/python/test => python/tests}/test_grid_utility.py (100%) rename {src/python/test => python/tests}/test_laplacian.py (100%) delete mode 100644 src/python/pyproject.toml diff --git a/examples/python/elliptic1D.py b/python/examples/elliptic1D.py similarity index 100% rename from examples/python/elliptic1D.py rename to python/examples/elliptic1D.py diff --git a/python/pyproject.toml b/python/pyproject.toml index e69de29bb..ebc9e7df3 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -0,0 +1,9 @@ +[tool.ruff] +src = ["."] +fix = true + +[tool.ruff.lint] +extend-select = ["I"] + +[tool.ruff.lint.isort] +known-first-party = ["pymole"] diff --git a/src/python/pymole/BoundaryCondition.py b/python/src/pymole/BoundaryCondition.py similarity index 100% rename from src/python/pymole/BoundaryCondition.py rename to python/src/pymole/BoundaryCondition.py diff --git a/src/python/pymole/RobinBoundaryCondition.py b/python/src/pymole/RobinBoundaryCondition.py similarity index 100% rename from src/python/pymole/RobinBoundaryCondition.py rename to python/src/pymole/RobinBoundaryCondition.py diff --git a/src/python/pymole/__init__.py b/python/src/pymole/__init__.py similarity index 100% rename from src/python/pymole/__init__.py rename to python/src/pymole/__init__.py diff --git a/src/python/pymole/curl.py b/python/src/pymole/curl.py similarity index 100% rename from src/python/pymole/curl.py rename to python/src/pymole/curl.py diff --git a/src/python/pymole/divergence.py b/python/src/pymole/divergence.py similarity index 100% rename from src/python/pymole/divergence.py rename to python/src/pymole/divergence.py diff --git a/src/python/pymole/gradient.py b/python/src/pymole/gradient.py similarity index 100% rename from src/python/pymole/gradient.py rename to python/src/pymole/gradient.py diff --git a/src/python/pymole/grid.py b/python/src/pymole/grid.py similarity index 100% rename from src/python/pymole/grid.py rename to python/src/pymole/grid.py diff --git a/src/python/pymole/laplacian.py b/python/src/pymole/laplacian.py similarity index 100% rename from src/python/pymole/laplacian.py rename to python/src/pymole/laplacian.py diff --git a/src/python/pymole/mimetic_operator.py b/python/src/pymole/mimetic_operator.py similarity index 100% rename from src/python/pymole/mimetic_operator.py rename to python/src/pymole/mimetic_operator.py diff --git a/src/python/pymole/requirements.txt b/python/src/pymole/requirements.txt similarity index 100% rename from src/python/pymole/requirements.txt rename to python/src/pymole/requirements.txt diff --git a/src/python/test/__init__.py b/python/tests/__init__.py similarity index 100% rename from src/python/test/__init__.py rename to python/tests/__init__.py diff --git a/src/python/test/test_curl.py b/python/tests/test_curl.py similarity index 100% rename from src/python/test/test_curl.py rename to python/tests/test_curl.py diff --git a/src/python/test/test_divergence.py b/python/tests/test_divergence.py similarity index 100% rename from src/python/test/test_divergence.py rename to python/tests/test_divergence.py diff --git a/src/python/test/test_gradient.py b/python/tests/test_gradient.py similarity index 100% rename from src/python/test/test_gradient.py rename to python/tests/test_gradient.py diff --git a/src/python/test/test_grid_utility.py b/python/tests/test_grid_utility.py similarity index 100% rename from src/python/test/test_grid_utility.py rename to python/tests/test_grid_utility.py diff --git a/src/python/test/test_laplacian.py b/python/tests/test_laplacian.py similarity index 100% rename from src/python/test/test_laplacian.py rename to python/tests/test_laplacian.py diff --git a/src/python/pyproject.toml b/src/python/pyproject.toml deleted file mode 100644 index ebc9e7df3..000000000 --- a/src/python/pyproject.toml +++ /dev/null @@ -1,9 +0,0 @@ -[tool.ruff] -src = ["."] -fix = true - -[tool.ruff.lint] -extend-select = ["I"] - -[tool.ruff.lint.isort] -known-first-party = ["pymole"] From 849e25828f0d189a279e8ade3430d87e80485131 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 11 Aug 2026 21:04:42 -0700 Subject: [PATCH 25/66] Update src/pymole/__init__.py --- .github/workflows/ci.yml | 2 +- python/pyproject.toml | 13 +++++++++++++ python/src/pymole/__init__.py | 11 ----------- 3 files changed, 14 insertions(+), 12 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ddaed3359..bb237db0d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -308,4 +308,4 @@ jobs: - name: Run tests run: | - pytest \ No newline at end of file + pytest python/tests diff --git a/python/pyproject.toml b/python/pyproject.toml index ebc9e7df3..d32ecb79e 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -7,3 +7,16 @@ extend-select = ["I"] [tool.ruff.lint.isort] known-first-party = ["pymole"] + +[build-system] +requires = ["setuptools>=68"] +build-backend = "setuptools.build_meta" + +[project] +name = "pymole" +version = "0.1.0" +requires-python = ">=3.11" + +[tool.setuptools.packages.find] +where = ["."] +include = ["pymole*"] diff --git a/python/src/pymole/__init__.py b/python/src/pymole/__init__.py index 538bbd252..e41afd6d0 100644 --- a/python/src/pymole/__init__.py +++ b/python/src/pymole/__init__.py @@ -1,6 +1,3 @@ -from .dim1d import Divergence1D, Gradient1D, Grid1D, Laplacian1D -from .dim2d import Divergence2D, Gradient2D, Laplacian2D -from .dim3d import Gradient3D from .grid import Grid from .curl import Curl from .divergence import Divergence @@ -14,17 +11,9 @@ "BoundaryCondition", "Curl", "Divergence", - "Divergence1D", - "Divergence2D", "Gradient", - "Gradient1D", - "Gradient2D", - "Gradient3D", "Grid", - "Grid1D", "Laplacian", - "Laplacian1D", - "Laplacian2D", "RobinBoundaryCondition", ] From b876a0c6ee729dc79ca8f7d877ebb01d67063f75 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 11 Aug 2026 22:04:09 -0700 Subject: [PATCH 26/66] Update pyproject.toml --- python/pyproject.toml | 7 +++---- python/tests/test_curl.py | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/python/pyproject.toml b/python/pyproject.toml index d32ecb79e..b3ca50c3c 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -1,5 +1,5 @@ [tool.ruff] -src = ["."] +src = ["src"] fix = true [tool.ruff.lint] @@ -15,8 +15,7 @@ build-backend = "setuptools.build_meta" [project] name = "pymole" version = "0.1.0" -requires-python = ">=3.11" +requires-python = ">=3.10" [tool.setuptools.packages.find] -where = ["."] -include = ["pymole*"] +where = ["src"] diff --git a/python/tests/test_curl.py b/python/tests/test_curl.py index b58642629..5de405f1a 100644 --- a/python/tests/test_curl.py +++ b/python/tests/test_curl.py @@ -3,7 +3,7 @@ import sys from pathlib import Path -sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) +#sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) from pymole import Curl, Divergence, Grid class TestCurl: From c18138e6b6eac1bea2e004584a6101b29cb8954e Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 11 Aug 2026 22:24:38 -0700 Subject: [PATCH 27/66] Add editable install for Python package Signed-off-by: Christopher Paolini --- .github/workflows/ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index bb237db0d..a711f2b33 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -304,6 +304,7 @@ jobs: python -m pip install --upgrade pip if [ -f requirements.txt ]; then pip install -r requirements.txt; fi pip install pytest numpy scipy + python -m pip install -e ./python pwd - name: Run tests From 59b54480fa6bcd5da34ec92c5e2d3532e0a21d08 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 11 Aug 2026 22:35:41 -0700 Subject: [PATCH 28/66] Create README.md for pymole Python interface Added README.md for pymole Python interface with installation and usage instructions. Signed-off-by: Christopher Paolini --- python/README.md | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/python/README.md b/python/README.md index e69de29bb..7056cd473 100644 --- a/python/README.md +++ b/python/README.md @@ -0,0 +1,35 @@ +# pymole: Python interface for [MOLE](https://github.com/csrc-sdsu/mole) + +## Description + +This is an implementation written in Python for the MOLE package. The source files for the operators are located at ```src/``` and examples can be found in ```examples/```. This is a WIP implementation; currently (as of August 11, 2026), only the 1D and 2D operators are available. + +## Installation + +pymole is not yet available in the Python Package Index (PyPI). For now, this repository needs to be cloned locally in order to use the library in Python. + +## Using pymole + +In order to use the pymole library, first navigate to the location where the repository has been cloned to. Then, go to the `mole/python` sub-directory. From here, you can access the library via the command line. + +### From the command line + +In `mole/python`, use the following command to install the MOLE module: + +```sh +python -m pip install --user ./python +``` + +## Running the test suite + +To run the unit tests, first install the MOLE module as described above, then run the command `python -m pytest python/tests` from the directory `mole/python`). The results of the unit tests should be displayed to your console. + +## Examples + +The MOLE library contains examples demonstrating how to use the operators, in a broad range of partial differential equations (PDEs). More information on the mathematical content can be found in the [main MOLE documentation](https://mole-docs.readthedocs.io/en/main/examples/index.html). + +Currently, the following examples are available in the MOLE Python package. + +- Elliptic Problems + - 1D Examples + - `elliptic1D`: A script that solves the 1D Poisson's equation with Robin boundary conditions using mimetic operators. From a41ceb76eeb43530e8491010dedf71020b4e9bdf Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 07:00:59 -0700 Subject: [PATCH 29/66] Add GitHub Actions workflow for Python documentation Signed-off-by: Christopher Paolini --- .github/workflows/python-docs.yml | 48 +++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 .github/workflows/python-docs.yml diff --git a/.github/workflows/python-docs.yml b/.github/workflows/python-docs.yml new file mode 100644 index 000000000..16ba38669 --- /dev/null +++ b/.github/workflows/python-docs.yml @@ -0,0 +1,48 @@ +name: Python Documentation + +on: + pull_request: + push: + branches: + - main + tags: '*' + +permissions: + contents: write + pull-requests: read + statuses: write + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + docbuild: + runs-on: ubuntu-24.04 + + steps: + - uses: actions/checkout@v6 + + - uses: actions/setup-python@v6 + with: + python-version: '3.12' + cache: 'pip' + + - name: Install dependencies + run: | + cd python/MOLE + python -m pip install --upgrade pip + pip install -e . + pip install -r docs/requirements.txt + + - name: Build documentation + run: | + cd python/MOLE + sphinx-build -W --keep-going -b html docs docs/_build/html + + - name: Deploy documentation + if: github.event_name == 'push' && github.ref == 'refs/heads/main' + uses: peaceiris/actions-gh-pages@v4 + with: + github_token: ${{ secrets.GITHUB_TOKEN }} + publish_dir: python/MOLE/docs/_build/html From d35a670dadc54582de12a38d9f7d144eb77a5538 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 07:01:46 -0700 Subject: [PATCH 30/66] Add Sphinx and theme to requirements Signed-off-by: Christopher Paolini --- python/docs/requirements.txt | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 python/docs/requirements.txt diff --git a/python/docs/requirements.txt b/python/docs/requirements.txt new file mode 100644 index 000000000..cbf1e3658 --- /dev/null +++ b/python/docs/requirements.txt @@ -0,0 +1,2 @@ +sphinx +sphinx-rtd-theme From dee35185fe0931854caee2c299765c3570842d19 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 07:03:29 -0700 Subject: [PATCH 31/66] Update directory path for dependency installation Signed-off-by: Christopher Paolini --- .github/workflows/python-docs.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python-docs.yml b/.github/workflows/python-docs.yml index 16ba38669..8b298fc20 100644 --- a/.github/workflows/python-docs.yml +++ b/.github/workflows/python-docs.yml @@ -30,7 +30,7 @@ jobs: - name: Install dependencies run: | - cd python/MOLE + cd python/src python -m pip install --upgrade pip pip install -e . pip install -r docs/requirements.txt From 67a800cf735d0903d88dd6a2ab64aa3a949e50e5 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 07:06:44 -0700 Subject: [PATCH 32/66] Remove cd command in documentation build step Remove unnecessary directory change before building documentation. Signed-off-by: Christopher Paolini --- .github/workflows/python-docs.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/python-docs.yml b/.github/workflows/python-docs.yml index 8b298fc20..1e98032a4 100644 --- a/.github/workflows/python-docs.yml +++ b/.github/workflows/python-docs.yml @@ -37,7 +37,6 @@ jobs: - name: Build documentation run: | - cd python/MOLE sphinx-build -W --keep-going -b html docs docs/_build/html - name: Deploy documentation From 277d274da5373ada59ea73a9828d4cf225540b1b Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 07:42:38 -0700 Subject: [PATCH 33/66] Modify Python docs workflow to trigger on all branches Updated workflow triggers to include all branches. Signed-off-by: Christopher Paolini --- .github/workflows/python-docs.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/python-docs.yml b/.github/workflows/python-docs.yml index 1e98032a4..b584834d0 100644 --- a/.github/workflows/python-docs.yml +++ b/.github/workflows/python-docs.yml @@ -2,9 +2,11 @@ name: Python Documentation on: pull_request: + branches: + - '**' push: branches: - - main + - '**' tags: '*' permissions: From fa40d5db6317eae2b99baf9ae9c5806b1ac6e7c9 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 09:16:53 -0700 Subject: [PATCH 34/66] Modify Python workflow to install dependencies Updated the workflow to conditionally install dependencies from requirements.txt and added installation of pytest, numpy, and scipy. Signed-off-by: Christopher Paolini --- .github/workflows/python-docs.yml | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/.github/workflows/python-docs.yml b/.github/workflows/python-docs.yml index b584834d0..b7b3d1f98 100644 --- a/.github/workflows/python-docs.yml +++ b/.github/workflows/python-docs.yml @@ -29,12 +29,13 @@ jobs: with: python-version: '3.12' cache: 'pip' - + - name: Install dependencies run: | - cd python/src python -m pip install --upgrade pip - pip install -e . + if [ -f requirements.txt ]; then pip install -r requirements.txt; fi + pip install pytest numpy scipy + python -m pip install -e ./python pip install -r docs/requirements.txt - name: Build documentation From 6dc9e36939a81f01540e9a1d73ff4d07d98db6d0 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 09:19:15 -0700 Subject: [PATCH 35/66] Add pwd command to Python documentation workflow Signed-off-by: Christopher Paolini --- .github/workflows/python-docs.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/python-docs.yml b/.github/workflows/python-docs.yml index b7b3d1f98..f2ea98a09 100644 --- a/.github/workflows/python-docs.yml +++ b/.github/workflows/python-docs.yml @@ -36,6 +36,7 @@ jobs: if [ -f requirements.txt ]; then pip install -r requirements.txt; fi pip install pytest numpy scipy python -m pip install -e ./python + pwd pip install -r docs/requirements.txt - name: Build documentation From ac4ead0494c26af5acb31b1dbd1c73425cfd71ca Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 09:25:37 -0700 Subject: [PATCH 36/66] Update workflow to install docs requirements in python dir Change directory to 'python' before installing documentation requirements. Signed-off-by: Christopher Paolini --- .github/workflows/python-docs.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/python-docs.yml b/.github/workflows/python-docs.yml index f2ea98a09..3dab1c299 100644 --- a/.github/workflows/python-docs.yml +++ b/.github/workflows/python-docs.yml @@ -37,7 +37,9 @@ jobs: pip install pytest numpy scipy python -m pip install -e ./python pwd - pip install -r docs/requirements.txt + cd python + pwd + if [ -f docs/requirements.txt ]; then pip install -r docs/requirements.txt; fi - name: Build documentation run: | From 87f1fd60c05367c09cf60337617df6484dd63414 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 09:27:55 -0700 Subject: [PATCH 37/66] Add pwd commands to documentation build step Add debug output for current directory during documentation build Signed-off-by: Christopher Paolini --- .github/workflows/python-docs.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.github/workflows/python-docs.yml b/.github/workflows/python-docs.yml index 3dab1c299..a79cd696e 100644 --- a/.github/workflows/python-docs.yml +++ b/.github/workflows/python-docs.yml @@ -43,6 +43,9 @@ jobs: - name: Build documentation run: | + pwd + cd python + pwd sphinx-build -W --keep-going -b html docs docs/_build/html - name: Deploy documentation From 08f8bb57a752d26a4ab80652fb616c64a7712861 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 09:34:12 -0700 Subject: [PATCH 38/66] Add Sphinx configuration for documentation Initial configuration for Sphinx documentation builder with project details and settings. Signed-off-by: Christopher Paolini --- python/docs/conf.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 python/docs/conf.py diff --git a/python/docs/conf.py b/python/docs/conf.py new file mode 100644 index 000000000..79138f95f --- /dev/null +++ b/python/docs/conf.py @@ -0,0 +1,33 @@ +# Configuration file for the Sphinx documentation builder. +# +# For the full list of built-in configuration values, see: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +import os +import sys + +# Add the project source directory to sys.path so autodoc can import it. +# Adjust this path if your package lives somewhere else. +sys.path.insert(0, os.path.abspath("..")) + +# -- Project information ----------------------------------------------------- + +project = "mole" +copyright = "2026" +author = "" + +# -- General configuration -------------------------------------------------- + +extensions = [ + "sphinx.ext.autodoc", + "sphinx.ext.napoleon", + "sphinx.ext.viewcode", +] + +templates_path = ["_templates"] +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] + +# -- Options for HTML output ------------------------------------------------ + +html_theme = "alabaster" +html_static_path = ["_static"] From 8ffa8e8f69cdfb28f27bf83a8cb0685c0192446b Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 09:37:27 -0700 Subject: [PATCH 39/66] Add initial documentation structure for mole Signed-off-by: Christopher Paolini --- python/docs/index.rst | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 python/docs/index.rst diff --git a/python/docs/index.rst b/python/docs/index.rst new file mode 100644 index 000000000..037c2bb47 --- /dev/null +++ b/python/docs/index.rst @@ -0,0 +1,8 @@ +mole +==== + +Welcome to the mole documentation. + +.. toctree:: + :maxdepth: 2 + :caption: Contents: From cd1fdacbfcd2f29ba6bd2e7124944333cbe249ca Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 09:42:53 -0700 Subject: [PATCH 40/66] Create .gitkeep Signed-off-by: Christopher Paolini --- python/docs/_static/.gitkeep | 1 + 1 file changed, 1 insertion(+) create mode 100644 python/docs/_static/.gitkeep diff --git a/python/docs/_static/.gitkeep b/python/docs/_static/.gitkeep new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/python/docs/_static/.gitkeep @@ -0,0 +1 @@ + From 8451f779eb8d0879fcf4f51a2779312614dcb319 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 09:49:11 -0700 Subject: [PATCH 41/66] Update publish directory in Python docs workflow Signed-off-by: Christopher Paolini --- .github/workflows/python-docs.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python-docs.yml b/.github/workflows/python-docs.yml index a79cd696e..6937a1834 100644 --- a/.github/workflows/python-docs.yml +++ b/.github/workflows/python-docs.yml @@ -53,4 +53,4 @@ jobs: uses: peaceiris/actions-gh-pages@v4 with: github_token: ${{ secrets.GITHUB_TOKEN }} - publish_dir: python/MOLE/docs/_build/html + publish_dir: python/docs/_build/html From d654b88fee09a7439e62d30f31c68db231f70c28 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 09:50:03 -0700 Subject: [PATCH 42/66] Create .gitkeep Signed-off-by: Christopher Paolini --- python/docs/_build/html/.gitkeep | 1 + 1 file changed, 1 insertion(+) create mode 100644 python/docs/_build/html/.gitkeep diff --git a/python/docs/_build/html/.gitkeep b/python/docs/_build/html/.gitkeep new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/python/docs/_build/html/.gitkeep @@ -0,0 +1 @@ + From 5d70636643cae857d87bde65a5f6e1f1a55e6382 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 09:54:45 -0700 Subject: [PATCH 43/66] Update documentation title for pymole module Signed-off-by: Christopher Paolini --- python/docs/index.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/docs/index.rst b/python/docs/index.rst index 037c2bb47..e49e0cf57 100644 --- a/python/docs/index.rst +++ b/python/docs/index.rst @@ -1,7 +1,7 @@ -mole -==== +Documentation for pymole, the Python MOLE module +================================================ -Welcome to the mole documentation. +Welcome to the pymole documentation. .. toctree:: :maxdepth: 2 From 365e87540893d6dd3f14b2f68702d0c1d024aadf Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 12 Aug 2026 10:06:14 -0700 Subject: [PATCH 44/66] Update deployment condition for documentation Signed-off-by: Christopher Paolini --- .github/workflows/python-docs.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python-docs.yml b/.github/workflows/python-docs.yml index 6937a1834..105849712 100644 --- a/.github/workflows/python-docs.yml +++ b/.github/workflows/python-docs.yml @@ -49,7 +49,7 @@ jobs: sphinx-build -W --keep-going -b html docs docs/_build/html - name: Deploy documentation - if: github.event_name == 'push' && github.ref == 'refs/heads/main' + if: github.event_name == 'push' uses: peaceiris/actions-gh-pages@v4 with: github_token: ${{ secrets.GITHUB_TOKEN }} From 97d8234c931fb20b5a84c45a8fbe6940d55e6b5b Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sat, 12 Sep 2026 17:50:51 -0700 Subject: [PATCH 45/66] Remove obsolete directory --- src/python/pymole/dim1d/__init__.py | 13 ---- src/python/pymole/dim1d/divergence.py | 36 ---------- src/python/pymole/dim1d/gradient.py | 38 ----------- src/python/pymole/dim1d/grid.py | 58 ---------------- src/python/pymole/dim1d/laplacian.py | 20 ------ src/python/pymole/dim1d/operator.py | 19 ------ src/python/pymole/dim2d/__init__.py | 11 --- src/python/pymole/dim2d/divergence.py | 26 -------- src/python/pymole/dim2d/gradient.py | 26 -------- src/python/pymole/dim2d/laplacian.py | 20 ------ src/python/pymole/dim2d/operator.py | 23 ------- src/python/pymole/dim3d/__init__.py | 7 -- src/python/pymole/dim3d/gradient.py | 30 --------- src/python/pymole/dim3d/operator.py | 33 --------- src/python/test/dim1d/__init__.py | 0 src/python/test/dim1d/test_divergence.py | 73 -------------------- src/python/test/dim1d/test_gradient.py | 80 ---------------------- src/python/test/dim1d/test_grid.py | 85 ------------------------ src/python/test/dim1d/test_laplacian.py | 54 --------------- src/python/test/dim1d/test_operator.py | 47 ------------- src/python/test/dim2d/__init__.py | 0 src/python/test/dim2d/test_divergence.py | 57 ---------------- src/python/test/dim2d/test_gradient.py | 63 ------------------ src/python/test/dim2d/test_laplacian.py | 68 ------------------- src/python/test/dim2d/test_operator.py | 59 ---------------- src/python/test/dim3d/__init__.py | 0 src/python/test/dim3d/test_gradient.py | 76 --------------------- src/python/test/dim3d/test_operator.py | 70 ------------------- 28 files changed, 1092 deletions(-) delete mode 100644 src/python/pymole/dim1d/__init__.py delete mode 100644 src/python/pymole/dim1d/divergence.py delete mode 100644 src/python/pymole/dim1d/gradient.py delete mode 100644 src/python/pymole/dim1d/grid.py delete mode 100644 src/python/pymole/dim1d/laplacian.py delete mode 100644 src/python/pymole/dim1d/operator.py delete mode 100644 src/python/pymole/dim2d/__init__.py delete mode 100644 src/python/pymole/dim2d/divergence.py delete mode 100644 src/python/pymole/dim2d/gradient.py delete mode 100644 src/python/pymole/dim2d/laplacian.py delete mode 100644 src/python/pymole/dim2d/operator.py delete mode 100644 src/python/pymole/dim3d/__init__.py delete mode 100644 src/python/pymole/dim3d/gradient.py delete mode 100644 src/python/pymole/dim3d/operator.py delete mode 100644 src/python/test/dim1d/__init__.py delete mode 100644 src/python/test/dim1d/test_divergence.py delete mode 100644 src/python/test/dim1d/test_gradient.py delete mode 100644 src/python/test/dim1d/test_grid.py delete mode 100644 src/python/test/dim1d/test_laplacian.py delete mode 100644 src/python/test/dim1d/test_operator.py delete mode 100644 src/python/test/dim2d/__init__.py delete mode 100644 src/python/test/dim2d/test_divergence.py delete mode 100644 src/python/test/dim2d/test_gradient.py delete mode 100644 src/python/test/dim2d/test_laplacian.py delete mode 100644 src/python/test/dim2d/test_operator.py delete mode 100644 src/python/test/dim3d/__init__.py delete mode 100644 src/python/test/dim3d/test_gradient.py delete mode 100644 src/python/test/dim3d/test_operator.py diff --git a/src/python/pymole/dim1d/__init__.py b/src/python/pymole/dim1d/__init__.py deleted file mode 100644 index 39651117d..000000000 --- a/src/python/pymole/dim1d/__init__.py +++ /dev/null @@ -1,13 +0,0 @@ -from .divergence import Divergence1D -from .gradient import Gradient1D -from .grid import Grid1D -from .laplacian import Laplacian1D -from .operator import MimeticOperator1D - -__all__ = [ - "Divergence1D", - "Gradient1D", - "Grid1D", - "Laplacian1D", - "MimeticOperator1D", -] diff --git a/src/python/pymole/dim1d/divergence.py b/src/python/pymole/dim1d/divergence.py deleted file mode 100644 index 45a37a6de..000000000 --- a/src/python/pymole/dim1d/divergence.py +++ /dev/null @@ -1,36 +0,0 @@ -import numpy as np -from scipy import sparse - -from .operator import MimeticOperator1D - - -class Divergence1D(MimeticOperator1D): - @classmethod - def minimum_cells_for_order(cls, accuracy_order: int) -> int: - return 2 * accuracy_order + 1 - - def _construct_matrix(self) -> sparse.csr_matrix: - match self.accuracy_order: - case 2: - return self._divergence_2nd_order() - case _: - raise ValueError( - f"Unsupported {self.__class__.__name__} " - f"order of accuracy: {self.accuracy_order}." - ) - - def _divergence_2nd_order(self) -> sparse.csr_matrix: - m = self.x_grid.num_cells - delta_x = self.x_grid.cell_size - - main_diag = np.ones(m + 1) - sub_diag = np.full(m + 1, -1.0) - - main_diag[0] = 0.0 - sub_diag[-1] = 0.0 - - D = sparse.diags( - [sub_diag, main_diag], offsets=[-1, 0], shape=(m + 2, m + 1) - ).tocsr() - - return D / delta_x diff --git a/src/python/pymole/dim1d/gradient.py b/src/python/pymole/dim1d/gradient.py deleted file mode 100644 index 144bdf1cb..000000000 --- a/src/python/pymole/dim1d/gradient.py +++ /dev/null @@ -1,38 +0,0 @@ -import numpy as np -from scipy import sparse - -from .operator import MimeticOperator1D - - -class Gradient1D(MimeticOperator1D): - @classmethod - def minimum_cells_for_order(cls, accuracy_order: int) -> int: - return 2 * accuracy_order - - def _construct_matrix(self) -> sparse.csr_matrix: - match self.accuracy_order: - case 2: - return self._gradient_2nd_order() - case _: - raise ValueError( - f"Unsupported {self.__class__.__name__} " - f"order of accuracy: {self.accuracy_order}." - ) - - def _gradient_2nd_order(self) -> sparse.csr_matrix: - m = self.x_grid.num_cells - delta_x = self.x_grid.cell_size - - main_diag = np.full(m + 1, -1.0) - super_diag = np.ones(m + 1) - - main_diag[-1] = 0.0 - - G = sparse.diags( - [main_diag, super_diag], offsets=[0, 1], shape=(m + 1, m + 2) - ).tolil() - - G[0, 0:3] = [-8.0 / 3.0, 3.0, -1.0 / 3.0] - G[m, m - 1 : m + 2] = [1.0 / 3.0, -3.0, 8.0 / 3.0] - - return G.tocsr() / delta_x diff --git a/src/python/pymole/dim1d/grid.py b/src/python/pymole/dim1d/grid.py deleted file mode 100644 index 33b585dbb..000000000 --- a/src/python/pymole/dim1d/grid.py +++ /dev/null @@ -1,58 +0,0 @@ -from math import isclose - - -class Grid1D: - def __init__(self, lower_bound: float, upper_bound: float, num_cells: int): - if lower_bound >= upper_bound: - raise ValueError( - f"Lower bound ({lower_bound}) value must be less than " - f"the upper bound ({upper_bound})", - ) - if num_cells < 1: - raise ValueError(f"Number of cells ({num_cells}) must be at least 1") - - self._lower_bound = lower_bound - self._upper_bound = upper_bound - self._num_cells = num_cells - - @classmethod - def from_number_of_cells(cls, num_cells: int): - if num_cells < 1: - raise ValueError(f"Number of cells ({num_cells}) must be at least 1") - - lower_bound = 0.0 - upper_bound = 1.0 - - return cls(lower_bound, upper_bound, num_cells) - - @property - def lower_bound(self) -> float: - return self._lower_bound - - @property - def upper_bound(self) -> float: - return self._upper_bound - - @property - def num_cells(self) -> int: - return self._num_cells - - @property - def cell_size(self) -> float: - return (self.upper_bound - self.lower_bound) / self.num_cells - - def __repr__(self) -> str: - return ( - f"Grid1D(lower_bound={self.lower_bound}, " - f"upper_bound={self.upper_bound}, " - f"num_cells={self.num_cells})" - ) - - def __eq__(self, other: object) -> bool: - if not isinstance(other, Grid1D): - return False - return ( - isclose(self.lower_bound, other.lower_bound) - and isclose(self.upper_bound, other.upper_bound) - and self.num_cells == other.num_cells - ) diff --git a/src/python/pymole/dim1d/laplacian.py b/src/python/pymole/dim1d/laplacian.py deleted file mode 100644 index 95ec9a163..000000000 --- a/src/python/pymole/dim1d/laplacian.py +++ /dev/null @@ -1,20 +0,0 @@ -from scipy import sparse - -from .divergence import Divergence1D -from .gradient import Gradient1D -from .operator import MimeticOperator1D - - -class Laplacian1D(MimeticOperator1D): - @classmethod - def minimum_cells_for_order(cls, accuracy_order: int) -> int: - return max( - Divergence1D.minimum_cells_for_order(accuracy_order), - Gradient1D.minimum_cells_for_order(accuracy_order), - ) - - def _construct_matrix(self) -> sparse.csr_matrix: - D = Divergence1D(self.x_grid, self.accuracy_order) - G = Gradient1D(self.x_grid, self.accuracy_order) - - return D.matrix @ G.matrix diff --git a/src/python/pymole/dim1d/operator.py b/src/python/pymole/dim1d/operator.py deleted file mode 100644 index 6809967df..000000000 --- a/src/python/pymole/dim1d/operator.py +++ /dev/null @@ -1,19 +0,0 @@ -from __future__ import annotations - -from typing import TYPE_CHECKING - -from pymole.mimetic_operator import MimeticOperator - -if TYPE_CHECKING: - from .grid import Grid1D - - -class MimeticOperator1D(MimeticOperator): - axis_names = ("x",) - - def __init__(self, x_grid: Grid1D, accuracy_order: int = 2) -> None: - super().__init__((x_grid,), accuracy_order) - - @property - def x_grid(self) -> Grid1D: - return self.grids[0] diff --git a/src/python/pymole/dim2d/__init__.py b/src/python/pymole/dim2d/__init__.py deleted file mode 100644 index 909b45cd5..000000000 --- a/src/python/pymole/dim2d/__init__.py +++ /dev/null @@ -1,11 +0,0 @@ -from .divergence import Divergence2D -from .gradient import Gradient2D -from .laplacian import Laplacian2D -from .operator import MimeticOperator2D - -__all__ = [ - "Divergence2D", - "Gradient2D", - "Laplacian2D", - "MimeticOperator2D", -] diff --git a/src/python/pymole/dim2d/divergence.py b/src/python/pymole/dim2d/divergence.py deleted file mode 100644 index e3fcfb0ee..000000000 --- a/src/python/pymole/dim2d/divergence.py +++ /dev/null @@ -1,26 +0,0 @@ -from scipy import sparse - -from pymole.dim1d.divergence import Divergence1D - -from .operator import MimeticOperator2D - - -class Divergence2D(MimeticOperator2D): - @classmethod - def minimum_cells_for_order(cls, accuracy_order: int) -> int: - return Divergence1D.minimum_cells_for_order(accuracy_order) - - def _construct_matrix(self) -> sparse.csr_matrix: - m = self.x_grid.num_cells - n = self.y_grid.num_cells - - Dx = Divergence1D(self.x_grid, self.accuracy_order).matrix - Dy = Divergence1D(self.y_grid, self.accuracy_order).matrix - - Im = sparse.eye(m + 2, format="csr")[:, 1:-1] - In = sparse.eye(n + 2, format="csr")[:, 1:-1] - - Sx = sparse.kron(In, Dx, format="csr") - Sy = sparse.kron(Dy, Im, format="csr") - - return sparse.hstack([Sx, Sy], format="csr") diff --git a/src/python/pymole/dim2d/gradient.py b/src/python/pymole/dim2d/gradient.py deleted file mode 100644 index 5cb5272fd..000000000 --- a/src/python/pymole/dim2d/gradient.py +++ /dev/null @@ -1,26 +0,0 @@ -from scipy import sparse - -from pymole.dim1d.gradient import Gradient1D - -from .operator import MimeticOperator2D - - -class Gradient2D(MimeticOperator2D): - @classmethod - def minimum_cells_for_order(cls, accuracy_order: int) -> int: - return Gradient1D.minimum_cells_for_order(accuracy_order) - - def _construct_matrix(self) -> sparse.csr_matrix: - m = self.x_grid.num_cells - n = self.y_grid.num_cells - - Gx = Gradient1D(self.x_grid, self.accuracy_order).matrix - Gy = Gradient1D(self.y_grid, self.accuracy_order).matrix - - Im = sparse.eye(m + 2, format="csr")[1:-1, :] - In = sparse.eye(n + 2, format="csr")[1:-1, :] - - Sx = sparse.kron(In, Gx, format="csr") - Sy = sparse.kron(Gy, Im, format="csr") - - return sparse.vstack([Sx, Sy], format="csr") diff --git a/src/python/pymole/dim2d/laplacian.py b/src/python/pymole/dim2d/laplacian.py deleted file mode 100644 index 997545c8a..000000000 --- a/src/python/pymole/dim2d/laplacian.py +++ /dev/null @@ -1,20 +0,0 @@ -from scipy import sparse - -from .divergence import Divergence2D -from .gradient import Gradient2D -from .operator import MimeticOperator2D - - -class Laplacian2D(MimeticOperator2D): - @classmethod - def minimum_cells_for_order(cls, accuracy_order: int) -> int: - return max( - Divergence2D.minimum_cells_for_order(accuracy_order), - Gradient2D.minimum_cells_for_order(accuracy_order), - ) - - def _construct_matrix(self) -> sparse.csr_matrix: - D = Divergence2D(self.x_grid, self.y_grid, self.accuracy_order) - G = Gradient2D(self.x_grid, self.y_grid, self.accuracy_order) - - return D.matrix @ G.matrix diff --git a/src/python/pymole/dim2d/operator.py b/src/python/pymole/dim2d/operator.py deleted file mode 100644 index de1161f91..000000000 --- a/src/python/pymole/dim2d/operator.py +++ /dev/null @@ -1,23 +0,0 @@ -from __future__ import annotations - -from typing import TYPE_CHECKING - -from pymole.mimetic_operator import MimeticOperator - -if TYPE_CHECKING: - from pymole.dim1d.grid import Grid1D - - -class MimeticOperator2D(MimeticOperator): - axis_names = ("x", "y") - - def __init__(self, x_grid: Grid1D, y_grid: Grid1D, accuracy_order: int = 2) -> None: - super().__init__((x_grid, y_grid), accuracy_order) - - @property - def x_grid(self) -> Grid1D: - return self.grids[0] - - @property - def y_grid(self) -> Grid1D: - return self.grids[1] diff --git a/src/python/pymole/dim3d/__init__.py b/src/python/pymole/dim3d/__init__.py deleted file mode 100644 index 38df6ef1f..000000000 --- a/src/python/pymole/dim3d/__init__.py +++ /dev/null @@ -1,7 +0,0 @@ -from .gradient import Gradient3D -from .operator import MimeticOperator3D - -__all__ = [ - "Gradient3D", - "MimeticOperator3D", -] diff --git a/src/python/pymole/dim3d/gradient.py b/src/python/pymole/dim3d/gradient.py deleted file mode 100644 index 0d4d12c93..000000000 --- a/src/python/pymole/dim3d/gradient.py +++ /dev/null @@ -1,30 +0,0 @@ -from scipy import sparse - -from pymole.dim1d.gradient import Gradient1D - -from .operator import MimeticOperator3D - - -class Gradient3D(MimeticOperator3D): - @classmethod - def minimum_cells_for_order(cls, accuracy_order: int) -> int: - return Gradient1D.minimum_cells_for_order(accuracy_order) - - def _construct_matrix(self) -> sparse.csr_matrix: - m = self.x_grid.num_cells - n = self.y_grid.num_cells - o = self.z_grid.num_cells - - Gx = Gradient1D(self.x_grid, self.accuracy_order).matrix - Gy = Gradient1D(self.y_grid, self.accuracy_order).matrix - Gz = Gradient1D(self.z_grid, self.accuracy_order).matrix - - Im = sparse.eye(m + 2, format="csr")[1:-1, :] - In = sparse.eye(n + 2, format="csr")[1:-1, :] - Io = sparse.eye(o + 2, format="csr")[1:-1, :] - - Sx = sparse.kron(sparse.kron(Io, In, format="csr"), Gx, format="csr") - Sy = sparse.kron(sparse.kron(Io, Gy, format="csr"), Im, format="csr") - Sz = sparse.kron(sparse.kron(Gz, In, format="csr"), Im, format="csr") - - return sparse.vstack([Sx, Sy, Sz], format="csr") diff --git a/src/python/pymole/dim3d/operator.py b/src/python/pymole/dim3d/operator.py deleted file mode 100644 index a58003305..000000000 --- a/src/python/pymole/dim3d/operator.py +++ /dev/null @@ -1,33 +0,0 @@ -from __future__ import annotations - -from typing import TYPE_CHECKING - -from pymole.mimetic_operator import MimeticOperator - -if TYPE_CHECKING: - from pymole.dim1d.grid import Grid1D - - -class MimeticOperator3D(MimeticOperator): - axis_names = ("x", "y", "z") - - def __init__( - self, - x_grid: Grid1D, - y_grid: Grid1D, - z_grid: Grid1D, - accuracy_order: int = 2, - ) -> None: - super().__init__((x_grid, y_grid, z_grid), accuracy_order) - - @property - def x_grid(self) -> Grid1D: - return self.grids[0] - - @property - def y_grid(self) -> Grid1D: - return self.grids[1] - - @property - def z_grid(self) -> Grid1D: - return self.grids[2] diff --git a/src/python/test/dim1d/__init__.py b/src/python/test/dim1d/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/src/python/test/dim1d/test_divergence.py b/src/python/test/dim1d/test_divergence.py deleted file mode 100644 index 9bde1f9a5..000000000 --- a/src/python/test/dim1d/test_divergence.py +++ /dev/null @@ -1,73 +0,0 @@ -import numpy as np -import pytest - -from pymole import Divergence1D -from pymole.dim1d.grid import Grid1D - - -@pytest.fixture -def a_second_order_divergence(): - grid = Grid1D(0.0, 1.2, 6) - return Divergence1D(grid, 2) - - -class TestDivergence1DConstruction: - def test_minimum_cells_for_divergence(self): - grid_too_small = Grid1D(0.0, 0.4, 4) - with pytest.raises(ValueError): - Divergence1D(grid_too_small, 2) - - grid_ok = Grid1D(0.0, 0.5, 5) - divergence = Divergence1D(grid_ok, 2) - assert divergence.x_grid.num_cells == 5 - - def test_unsupported_accuracy_order_raises_error(self): - grid = Grid1D(0.0, 1.0, 10) - with pytest.raises(ValueError): - divergence = Divergence1D(grid, 4) - _ = divergence.matrix - - -class TestDivergence1DMatrix: - def test_matrix_has_correct_shape(self, a_second_order_divergence): - D = a_second_order_divergence.matrix - num_cells = a_second_order_divergence.x_grid.num_cells - - assert D.shape == (num_cells + 2, num_cells + 1) - - def test_constant_field_divergence_is_zero(self, a_second_order_divergence): - num_cells = a_second_order_divergence.x_grid.num_cells - v_constant = np.ones(num_cells + 1) - D = a_second_order_divergence.matrix - div_v = D @ v_constant - - assert np.allclose(div_v, 0.0, atol=1e-12) - - -class TestDivergence1DSecondOrder: - def test_second_order_boundary_cells(self): - grid = Grid1D(0.0, 1.25, 5) - divergence = Divergence1D(grid, 2) - num_cells = grid.num_cells - - D = divergence.matrix.toarray() - - expected_first_row = np.zeros(num_cells + 1) - expected_last_row = np.zeros(num_cells + 1) - - assert np.allclose(D[0, :], expected_first_row, atol=1e-12) - assert np.allclose(D[-1, :], expected_last_row, atol=1e-12) - - def test_second_order_interior_stencil(self, a_second_order_divergence): - num_cells = a_second_order_divergence.x_grid.num_cells - delta_x = a_second_order_divergence.x_grid.cell_size - - D = a_second_order_divergence.matrix.toarray() - - for i in range(1, num_cells + 1): - row = D[i, :] - non_zero_indices = np.where(np.abs(row) > 1e-12)[0] - - assert len(non_zero_indices) == 2 - assert row[i - 1] == pytest.approx(-1.0 / delta_x) - assert row[i] == pytest.approx(1.0 / delta_x) diff --git a/src/python/test/dim1d/test_gradient.py b/src/python/test/dim1d/test_gradient.py deleted file mode 100644 index 0c7aec0be..000000000 --- a/src/python/test/dim1d/test_gradient.py +++ /dev/null @@ -1,80 +0,0 @@ -import numpy as np -import pytest - -from pymole import Gradient1D -from pymole.dim1d.grid import Grid1D - - -@pytest.fixture -def a_second_order_gradient(): - grid = Grid1D(0.0, 3.0, 15) - return Gradient1D(grid, 2) - - -class TestGradient1DConstruction: - def test_minimum_cells_for_gradient(self): - grid_too_small = Grid1D(0.0, 0.3, 3) - with pytest.raises(ValueError): - Gradient1D(grid_too_small, 2) - - grid_ok = Grid1D(0.0, 0.4, 4) - gradient = Gradient1D(grid_ok, 2) - - assert gradient.x_grid.num_cells == 4 - - def test_unsupported_accuracy_order_raises_error(self): - grid = Grid1D(0.0, 1.0, 10) - with pytest.raises(ValueError): - gradient = Gradient1D(grid, 4) - _ = gradient.matrix - - -class TestGradient1DMatrix: - def test_matrix_has_correct_shape(self, a_second_order_gradient): - G = a_second_order_gradient.matrix - num_cells = a_second_order_gradient.x_grid.num_cells - - assert G.shape == (num_cells + 1, num_cells + 2) - - def test_constant_field_gradient_is_zero(self, a_second_order_gradient): - num_cells = a_second_order_gradient.x_grid.num_cells - f_constant = np.ones(num_cells + 2) - G = a_second_order_gradient.matrix - grad_f = G @ f_constant - - assert np.allclose(grad_f, 0.0, atol=1e-12) - - def test_row_sums_are_zero(self, a_second_order_gradient): - G = a_second_order_gradient.matrix - row_sums = G.sum(axis=1).A1 - - assert np.allclose(row_sums, 0.0, atol=1e-12) - - -class TestGradient1DSecondOrder: - def test_second_order_boundary_stencil(self): - grid = Grid1D(0.0, 1.0, 4) - gradient = Gradient1D(grid, 2) - delta_x = grid.cell_size - - G = gradient.matrix.toarray() - - expected_first_row = np.array([-8.0 / 3.0, 3.0, -1.0 / 3.0, 0, 0, 0]) / delta_x - expected_last_row = np.array([0, 0, 0, 1.0 / 3.0, -3.0, 8.0 / 3.0]) / delta_x - - assert np.allclose(G[0, :], expected_first_row, atol=1e-12) - assert np.allclose(G[-1, :], expected_last_row, atol=1e-12) - - def test_second_order_interior_stencil(self, a_second_order_gradient): - num_cells = a_second_order_gradient.x_grid.num_cells - delta_x = a_second_order_gradient.x_grid.cell_size - - G = a_second_order_gradient.matrix.toarray() - - for i in range(1, num_cells): - row = G[i, :] - non_zero_indices = np.where(np.abs(row) > 1e-12)[0] - - assert len(non_zero_indices) == 2 - assert row[i] == pytest.approx(-1.0 / delta_x) - assert row[i + 1] == pytest.approx(1.0 / delta_x) diff --git a/src/python/test/dim1d/test_grid.py b/src/python/test/dim1d/test_grid.py deleted file mode 100644 index d4d55ccf6..000000000 --- a/src/python/test/dim1d/test_grid.py +++ /dev/null @@ -1,85 +0,0 @@ -import pytest - -from pymole.dim1d.grid import Grid1D - - -@pytest.fixture -def basic_grid(): - return Grid1D(0.0, 1.0, 10) - - -class TestGrid1DConstruction: - def test_valid_grid_creation(self, basic_grid): - assert basic_grid.lower_bound == pytest.approx(0.0) - assert basic_grid.upper_bound == pytest.approx(1.0) - assert basic_grid.num_cells == 10 - assert basic_grid.cell_size == pytest.approx(0.1) - - def test_lower_bound_greater_than_upper_bound_raises_error(self): - with pytest.raises(ValueError): - Grid1D(1.0, 0.0, 10) - - def test_num_cells_less_than_one_raises_error(self): - with pytest.raises(ValueError): - Grid1D(0.0, 1.0, 0) - - -class TestGrid1DFromNumberOfCells: - def test_from_number_of_cells_basic(self): - expected_num_cells = 13 - grid = Grid1D.from_number_of_cells(expected_num_cells) - - assert grid.lower_bound == pytest.approx(0.0) - assert grid.upper_bound == pytest.approx(1.0) - assert grid.num_cells == expected_num_cells - assert grid.cell_size == pytest.approx(1.0 / expected_num_cells) - - def test_zero_number_of_cells_raises_error(self): - with pytest.raises(ValueError): - Grid1D.from_number_of_cells(0) - - def test_negative_number_of_cells_raises_error(self): - with pytest.raises(ValueError): - Grid1D.from_number_of_cells(-5) - - -class TestGrid1DProperties: - def test_cell_size_computed(self): - grid = Grid1D(0.0, 2.0, 4) - - assert grid.cell_size == pytest.approx(0.5) - - -class TestGrid1DEquality: - def test_equal_grids(self, basic_grid): - grid2 = Grid1D(0.0, 1.0, 10) - - assert basic_grid == grid2 - - def test_unequal_grids_different_lower_bound(self, basic_grid): - grid2 = Grid1D(0.5, 1.0, 10) - - assert basic_grid != grid2 - - def test_unequal_grids_different_upper_bound(self, basic_grid): - grid2 = Grid1D(0.0, 2.0, 10) - - assert basic_grid != grid2 - - def test_unequal_grids_different_num_cells(self, basic_grid): - grid2 = Grid1D(0.0, 1.0, 20) - - assert basic_grid != grid2 - - def test_grid_not_equal_to_non_grid(self, basic_grid): - assert basic_grid != "not a grid" - - -class TestGrid1DRepr: - def test_repr_strings(self, basic_grid): - repr_str = repr(basic_grid) - - assert "Grid1D" in repr_str - assert "lower_bound=0.0" in repr_str - assert "upper_bound=1.0" in repr_str - assert "num_cells=10" in repr_str diff --git a/src/python/test/dim1d/test_laplacian.py b/src/python/test/dim1d/test_laplacian.py deleted file mode 100644 index 699ebae76..000000000 --- a/src/python/test/dim1d/test_laplacian.py +++ /dev/null @@ -1,54 +0,0 @@ -import numpy as np -import pytest - -from pymole import Laplacian1D -from pymole.dim1d.divergence import Divergence1D -from pymole.dim1d.gradient import Gradient1D -from pymole.dim1d.grid import Grid1D - - -@pytest.fixture -def a_second_order_laplacian(): - grid = Grid1D(0.0, 1.2, 6) - return Laplacian1D(grid, 2) - - -class TestLaplacian1DConstruction: - def test_minimum_cells_for_laplacian(self): - grid_too_small = Grid1D(0.0, 0.4, 4) - with pytest.raises(ValueError): - Laplacian1D(grid_too_small, 2) - - grid_ok = Grid1D(0.0, 0.5, 5) - laplacian = Laplacian1D(grid_ok, 2) - assert laplacian.x_grid.num_cells == 5 - - def test_unsupported_accuracy_order_raises_error(self): - grid = Grid1D(0.0, 1.0, 10) - with pytest.raises(ValueError): - laplacian = Laplacian1D(grid, 4) - _ = laplacian.matrix - - -class TestLaplacian1DMatrix: - def test_matrix_has_correct_shape(self, a_second_order_laplacian): - L = a_second_order_laplacian.matrix - num_cells = a_second_order_laplacian.x_grid.num_cells - - assert L.shape == (num_cells + 2, num_cells + 2) - - def test_constant_field_laplacian_is_zero(self, a_second_order_laplacian): - num_cells = a_second_order_laplacian.x_grid.num_cells - f_constant = np.ones(num_cells + 2) - L = a_second_order_laplacian.matrix - lap_f = L @ f_constant - - assert np.allclose(lap_f, 0.0, atol=1e-12) - - def test_laplacian_equals_divergence_times_gradient(self, a_second_order_laplacian): - grid = a_second_order_laplacian.x_grid - D = Divergence1D(grid, 2).matrix - G = Gradient1D(grid, 2).matrix - L = a_second_order_laplacian.matrix - - assert np.allclose((L - D @ G).toarray(), 0.0, atol=1e-12) diff --git a/src/python/test/dim1d/test_operator.py b/src/python/test/dim1d/test_operator.py deleted file mode 100644 index 77285bf1a..000000000 --- a/src/python/test/dim1d/test_operator.py +++ /dev/null @@ -1,47 +0,0 @@ -from unittest.mock import MagicMock - -import pytest -from scipy import sparse - -from pymole.dim1d.grid import Grid1D -from pymole.dim1d.operator import MimeticOperator1D - - -class MockOperator(MimeticOperator1D): - @classmethod - def minimum_cells_for_order(cls, accuracy_order: int) -> int: - return accuracy_order - - def _construct_matrix(self) -> sparse.csr_matrix: - return MagicMock(spec=sparse.csr_matrix) - - -class TestMimeticOperatorConstruction: - def test_minimum_cells_requirement(self): - grid_too_small = Grid1D(0.0, 0.3, 3) - with pytest.raises(ValueError): - MockOperator(grid_too_small, 4) - - grid_ok = Grid1D(0.0, 0.4, 4) - op = MockOperator(grid_ok, 2) - - assert op.x_grid == grid_ok - assert op.accuracy_order == 2 - assert op.minimum_cells == 2 - - -class TestMimeticOperatorMatrix: - def test_matrix_is_sparse_csr(self): - grid = Grid1D(0.0, 0.4, 4) - op = MockOperator(grid, 2) - matrix = op.matrix - - assert isinstance(matrix, MagicMock) - - def test_matrix_cached_property(self): - grid = Grid1D(0.0, 0.4, 4) - op = MockOperator(grid, 2) - matrix1 = op.matrix - matrix2 = op.matrix - - assert matrix1 is matrix2 diff --git a/src/python/test/dim2d/__init__.py b/src/python/test/dim2d/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/src/python/test/dim2d/test_divergence.py b/src/python/test/dim2d/test_divergence.py deleted file mode 100644 index 2d9246d3a..000000000 --- a/src/python/test/dim2d/test_divergence.py +++ /dev/null @@ -1,57 +0,0 @@ -import numpy as np -import pytest - -from pymole.dim1d.grid import Grid1D -from pymole.dim2d.divergence import Divergence2D - - -@pytest.fixture -def a_second_order_divergence(): - x_grid = Grid1D(0.0, 3.0, 6) - y_grid = Grid1D(0.0, 2.0, 5) - return Divergence2D(x_grid, y_grid, 2) - - -class TestDivergence2DConstruction: - def test_minimum_cells_for_divergence(self): - x_grid_too_small = Grid1D(0.0, 0.4, 4) - y_grid_ok = Grid1D(0.0, 1.0, 5) - with pytest.raises(ValueError): - Divergence2D(x_grid_too_small, y_grid_ok, 2) - - x_grid_ok = Grid1D(0.0, 0.5, 5) - y_grid_too_small = Grid1D(0.0, 0.4, 4) - with pytest.raises(ValueError): - Divergence2D(x_grid_ok, y_grid_too_small, 2) - - x_grid = Grid1D(0.0, 0.5, 5) - y_grid = Grid1D(0.0, 0.5, 5) - divergence = Divergence2D(x_grid, y_grid, 2) - - assert divergence.x_grid.num_cells == 5 - assert divergence.y_grid.num_cells == 5 - - def test_unsupported_accuracy_order_raises_error(self): - x_grid = Grid1D(0.0, 1.0, 10) - y_grid = Grid1D(0.0, 1.0, 10) - with pytest.raises(ValueError): - divergence = Divergence2D(x_grid, y_grid, 4) - _ = divergence.matrix - - -class TestDivergence2DMatrix: - def test_matrix_has_correct_shape(self, a_second_order_divergence): - m = a_second_order_divergence.x_grid.num_cells - n = a_second_order_divergence.y_grid.num_cells - D = a_second_order_divergence.matrix - - assert D.shape == ((m + 2) * (n + 2), 2 * m * n + m + n) - - def test_constant_field_divergence_is_zero(self, a_second_order_divergence): - m = a_second_order_divergence.x_grid.num_cells - n = a_second_order_divergence.y_grid.num_cells - v_constant = np.ones(2 * m * n + m + n) - D = a_second_order_divergence.matrix - div_v = D @ v_constant - - assert np.allclose(div_v, 0.0, atol=1e-12) diff --git a/src/python/test/dim2d/test_gradient.py b/src/python/test/dim2d/test_gradient.py deleted file mode 100644 index 6f1ac86e0..000000000 --- a/src/python/test/dim2d/test_gradient.py +++ /dev/null @@ -1,63 +0,0 @@ -import numpy as np -import pytest - -from pymole.dim1d.grid import Grid1D -from pymole.dim2d.gradient import Gradient2D - - -@pytest.fixture -def a_second_order_gradient(): - x_grid = Grid1D(0.0, 3.0, 6) - y_grid = Grid1D(0.0, 2.0, 5) - return Gradient2D(x_grid, y_grid, 2) - - -class TestGradient2DConstruction: - def test_minimum_cells_for_gradient(self): - x_grid_too_small = Grid1D(0.0, 0.3, 3) - y_grid_ok = Grid1D(0.0, 1.0, 4) - with pytest.raises(ValueError): - Gradient2D(x_grid_too_small, y_grid_ok, 2) - - x_grid_ok = Grid1D(0.0, 0.4, 4) - y_grid_too_small = Grid1D(0.0, 0.3, 3) - with pytest.raises(ValueError): - Gradient2D(x_grid_ok, y_grid_too_small, 2) - - x_grid = Grid1D(0.0, 0.4, 4) - y_grid = Grid1D(0.0, 0.4, 4) - gradient = Gradient2D(x_grid, y_grid, 2) - - assert gradient.x_grid.num_cells == 4 - assert gradient.y_grid.num_cells == 4 - - def test_unsupported_accuracy_order_raises_error(self): - x_grid = Grid1D(0.0, 1.0, 10) - y_grid = Grid1D(0.0, 1.0, 10) - with pytest.raises(ValueError): - gradient = Gradient2D(x_grid, y_grid, 4) - _ = gradient.matrix - - -class TestGradient2DMatrix: - def test_matrix_has_correct_shape(self, a_second_order_gradient): - m = a_second_order_gradient.x_grid.num_cells - n = a_second_order_gradient.y_grid.num_cells - G = a_second_order_gradient.matrix - - assert G.shape == (2 * m * n + m + n, (m + 2) * (n + 2)) - - def test_constant_field_gradient_is_zero(self, a_second_order_gradient): - m = a_second_order_gradient.x_grid.num_cells - n = a_second_order_gradient.y_grid.num_cells - f_constant = np.ones((m + 2) * (n + 2)) - G = a_second_order_gradient.matrix - grad_f = G @ f_constant - - assert np.allclose(grad_f, 0.0, atol=1e-12) - - def test_row_sums_are_zero(self, a_second_order_gradient): - G = a_second_order_gradient.matrix - row_sums = G.sum(axis=1).A1 - - assert np.allclose(row_sums, 0.0, atol=1e-12) diff --git a/src/python/test/dim2d/test_laplacian.py b/src/python/test/dim2d/test_laplacian.py deleted file mode 100644 index 5c336caaa..000000000 --- a/src/python/test/dim2d/test_laplacian.py +++ /dev/null @@ -1,68 +0,0 @@ -import numpy as np -import pytest - -from pymole.dim1d.grid import Grid1D -from pymole.dim2d.divergence import Divergence2D -from pymole.dim2d.gradient import Gradient2D -from pymole.dim2d.laplacian import Laplacian2D - - -@pytest.fixture -def a_second_order_laplacian(): - x_grid = Grid1D(0.0, 3.0, 6) - y_grid = Grid1D(0.0, 2.0, 5) - return Laplacian2D(x_grid, y_grid, 2) - - -class TestLaplacian2DConstruction: - def test_minimum_cells_for_laplacian(self): - x_grid_too_small = Grid1D(0.0, 0.4, 4) - y_grid_ok = Grid1D(0.0, 1.0, 5) - with pytest.raises(ValueError): - Laplacian2D(x_grid_too_small, y_grid_ok, 2) - - x_grid_ok = Grid1D(0.0, 0.5, 5) - y_grid_too_small = Grid1D(0.0, 0.4, 4) - with pytest.raises(ValueError): - Laplacian2D(x_grid_ok, y_grid_too_small, 2) - - x_grid = Grid1D(0.0, 0.5, 5) - y_grid = Grid1D(0.0, 0.5, 5) - laplacian = Laplacian2D(x_grid, y_grid, 2) - - assert laplacian.x_grid.num_cells == 5 - assert laplacian.y_grid.num_cells == 5 - - def test_unsupported_accuracy_order_raises_error(self): - x_grid = Grid1D(0.0, 1.0, 10) - y_grid = Grid1D(0.0, 1.0, 10) - with pytest.raises(ValueError): - laplacian = Laplacian2D(x_grid, y_grid, 4) - _ = laplacian.matrix - - -class TestLaplacian2DMatrix: - def test_matrix_has_correct_shape(self, a_second_order_laplacian): - m = a_second_order_laplacian.x_grid.num_cells - n = a_second_order_laplacian.y_grid.num_cells - L = a_second_order_laplacian.matrix - - assert L.shape == ((m + 2) * (n + 2), (m + 2) * (n + 2)) - - def test_constant_field_laplacian_is_zero(self, a_second_order_laplacian): - m = a_second_order_laplacian.x_grid.num_cells - n = a_second_order_laplacian.y_grid.num_cells - f_constant = np.ones((m + 2) * (n + 2)) - L = a_second_order_laplacian.matrix - lap_f = L @ f_constant - - assert np.allclose(lap_f, 0.0, atol=1e-12) - - def test_laplacian_equals_divergence_times_gradient(self, a_second_order_laplacian): - x_grid = a_second_order_laplacian.x_grid - y_grid = a_second_order_laplacian.y_grid - D = Divergence2D(x_grid, y_grid, 2).matrix - G = Gradient2D(x_grid, y_grid, 2).matrix - L = a_second_order_laplacian.matrix - - assert np.allclose((L - D @ G).toarray(), 0.0, atol=1e-12) diff --git a/src/python/test/dim2d/test_operator.py b/src/python/test/dim2d/test_operator.py deleted file mode 100644 index 58ef9087d..000000000 --- a/src/python/test/dim2d/test_operator.py +++ /dev/null @@ -1,59 +0,0 @@ -from unittest.mock import MagicMock - -import pytest -from scipy import sparse - -from pymole.dim1d.grid import Grid1D -from pymole.dim2d.operator import MimeticOperator2D - - -class MockOperator2D(MimeticOperator2D): - @classmethod - def minimum_cells_for_order(cls, accuracy_order: int) -> int: - return accuracy_order - - def _construct_matrix(self) -> sparse.csr_matrix: - return MagicMock(spec=sparse.csr_matrix) - - -@pytest.fixture -def x_grid(): - return Grid1D(0.0, 1.0, 4) - - -@pytest.fixture -def y_grid(): - return Grid1D(0.0, 1.0, 4) - - -class TestMimeticOperator2DConstruction: - def test_x_grid_too_small_raises_error(self, y_grid): - grid_too_small = Grid1D(0.0, 0.3, 3) - with pytest.raises(ValueError): - MockOperator2D(grid_too_small, y_grid, 4) - - def test_y_grid_too_small_raises_error(self, x_grid): - grid_too_small = Grid1D(0.0, 0.3, 3) - with pytest.raises(ValueError): - MockOperator2D(x_grid, grid_too_small, 4) - - def test_stores_grids_and_accuracy_order(self, x_grid, y_grid): - op = MockOperator2D(x_grid, y_grid, 2) - - assert op.x_grid == x_grid - assert op.y_grid == y_grid - assert op.accuracy_order == 2 - - def test_minimum_cells_property(self, x_grid, y_grid): - op = MockOperator2D(x_grid, y_grid, 2) - - assert op.minimum_cells == 2 - - -class TestMimeticOperator2DMatrix: - def test_matrix_cached_property(self, x_grid, y_grid): - op = MockOperator2D(x_grid, y_grid, 2) - matrix1 = op.matrix - matrix2 = op.matrix - - assert matrix1 is matrix2 diff --git a/src/python/test/dim3d/__init__.py b/src/python/test/dim3d/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/src/python/test/dim3d/test_gradient.py b/src/python/test/dim3d/test_gradient.py deleted file mode 100644 index 1087a1616..000000000 --- a/src/python/test/dim3d/test_gradient.py +++ /dev/null @@ -1,76 +0,0 @@ -import numpy as np -import pytest - -from pymole.dim1d.grid import Grid1D -from pymole.dim3d.gradient import Gradient3D - - -@pytest.fixture -def a_second_order_gradient(): - x_grid = Grid1D(0.0, 3.0, 6) - y_grid = Grid1D(0.0, 2.0, 5) - z_grid = Grid1D(0.0, 1.5, 4) - return Gradient3D(x_grid, y_grid, z_grid, 2) - - -class TestGradient3DConstruction: - def test_minimum_cells_for_gradient(self): - x_grid_too_small = Grid1D(0.0, 0.3, 3) - y_grid_ok = Grid1D(0.0, 1.0, 4) - z_grid_ok = Grid1D(0.0, 1.0, 4) - with pytest.raises(ValueError): - Gradient3D(x_grid_too_small, y_grid_ok, z_grid_ok, 2) - - x_grid_ok = Grid1D(0.0, 0.4, 4) - y_grid_too_small = Grid1D(0.0, 0.3, 3) - with pytest.raises(ValueError): - Gradient3D(x_grid_ok, y_grid_too_small, z_grid_ok, 2) - - z_grid_too_small = Grid1D(0.0, 0.3, 3) - with pytest.raises(ValueError): - Gradient3D(x_grid_ok, y_grid_ok, z_grid_too_small, 2) - - x_grid = Grid1D(0.0, 0.4, 4) - y_grid = Grid1D(0.0, 0.4, 4) - z_grid = Grid1D(0.0, 0.4, 4) - gradient = Gradient3D(x_grid, y_grid, z_grid, 2) - - assert gradient.x_grid.num_cells == 4 - assert gradient.y_grid.num_cells == 4 - assert gradient.z_grid.num_cells == 4 - - def test_unsupported_accuracy_order_raises_error(self): - x_grid = Grid1D(0.0, 1.0, 10) - y_grid = Grid1D(0.0, 1.0, 10) - z_grid = Grid1D(0.0, 1.0, 10) - with pytest.raises(ValueError): - gradient = Gradient3D(x_grid, y_grid, z_grid, 4) - _ = gradient.matrix - - -class TestGradient3DMatrix: - def test_matrix_has_correct_shape(self, a_second_order_gradient): - m = a_second_order_gradient.x_grid.num_cells - n = a_second_order_gradient.y_grid.num_cells - o = a_second_order_gradient.z_grid.num_cells - G = a_second_order_gradient.matrix - - expected_rows = 3 * m * n * o + m * n + m * o + n * o - expected_cols = (m + 2) * (n + 2) * (o + 2) - assert G.shape == (expected_rows, expected_cols) - - def test_constant_field_gradient_is_zero(self, a_second_order_gradient): - m = a_second_order_gradient.x_grid.num_cells - n = a_second_order_gradient.y_grid.num_cells - o = a_second_order_gradient.z_grid.num_cells - f_constant = np.ones((m + 2) * (n + 2) * (o + 2)) - G = a_second_order_gradient.matrix - grad_f = G @ f_constant - - assert np.allclose(grad_f, 0.0, atol=1e-12) - - def test_row_sums_are_zero(self, a_second_order_gradient): - G = a_second_order_gradient.matrix - row_sums = G.sum(axis=1).A1 - - assert np.allclose(row_sums, 0.0, atol=1e-12) diff --git a/src/python/test/dim3d/test_operator.py b/src/python/test/dim3d/test_operator.py deleted file mode 100644 index 91e10e5b1..000000000 --- a/src/python/test/dim3d/test_operator.py +++ /dev/null @@ -1,70 +0,0 @@ -from unittest.mock import MagicMock - -import pytest -from scipy import sparse - -from pymole.dim1d.grid import Grid1D -from pymole.dim3d.operator import MimeticOperator3D - - -class MockOperator3D(MimeticOperator3D): - @classmethod - def minimum_cells_for_order(cls, accuracy_order: int) -> int: - return accuracy_order - - def _construct_matrix(self) -> sparse.csr_matrix: - return MagicMock(spec=sparse.csr_matrix) - - -@pytest.fixture -def x_grid(): - return Grid1D(-9.4, 2.6, 12) - - -@pytest.fixture -def y_grid(): - return Grid1D(0.0, 1.0, 7) - - -@pytest.fixture -def z_grid(): - return Grid1D(0.0, 1.0, 11) - - -class TestMimeticOperator3DConstruction: - def test_x_grid_too_small_raises_error(self, y_grid, z_grid): - grid_too_small = Grid1D(0.0, 0.3, 3) - with pytest.raises(ValueError): - MockOperator3D(grid_too_small, y_grid, z_grid, 4) - - def test_y_grid_too_small_raises_error(self, x_grid, z_grid): - grid_too_small = Grid1D(0.0, 0.3, 3) - with pytest.raises(ValueError): - MockOperator3D(x_grid, grid_too_small, z_grid, 4) - - def test_z_grid_too_small_raises_error(self, x_grid, y_grid): - grid_too_small = Grid1D(0.0, 0.3, 3) - with pytest.raises(ValueError): - MockOperator3D(x_grid, y_grid, grid_too_small, 4) - - def test_stores_grids_and_accuracy_order(self, x_grid, y_grid, z_grid): - op = MockOperator3D(x_grid, y_grid, z_grid, 2) - - assert op.x_grid == x_grid - assert op.y_grid == y_grid - assert op.z_grid == z_grid - assert op.accuracy_order == 2 - - def test_minimum_cells_property(self, x_grid, y_grid, z_grid): - op = MockOperator3D(x_grid, y_grid, z_grid, 2) - - assert op.minimum_cells == 2 - - -class TestMimeticOperator3DMatrix: - def test_matrix_cached_property(self, x_grid, y_grid, z_grid): - op = MockOperator3D(x_grid, y_grid, z_grid, 2) - matrix1 = op.matrix - matrix2 = op.matrix - - assert matrix1 is matrix2 From 5d67e32c08286ed15ae06c156d09e7af6e61faaf Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sat, 19 Sep 2026 23:01:42 -0700 Subject: [PATCH 46/66] Update python/README.md --- python/README.md | 38 ++++++++++++++++++++++++++++++++++---- 1 file changed, 34 insertions(+), 4 deletions(-) diff --git a/python/README.md b/python/README.md index 7056cd473..bb83ebb76 100644 --- a/python/README.md +++ b/python/README.md @@ -8,13 +8,39 @@ This is an implementation written in Python for the MOLE package. The source fil pymole is not yet available in the Python Package Index (PyPI). For now, this repository needs to be cloned locally in order to use the library in Python. +## Python 3 + +The pymole library requires Python 3. You can check your Python version by running the command `python --version` or `python3 --version` in your terminal. If you do not have Python 3 installed, you can install it by following the instructions at [python.org](https://www.python.org/downloads/). + +If your operating system comes with Python 2 pre-installed, you may need to use `python3` instead of `python` in the commands below. In this case, you can also use the `python3` command to access your Python 3 installation. + ## Using pymole -In order to use the pymole library, first navigate to the location where the repository has been cloned to. Then, go to the `mole/python` sub-directory. From here, you can access the library via the command line. +In order to use the pymole library, first navigate to the location where the repository has been cloned to. From here, you can access the library via the command line. + +### Dependencies + +MOLE requires the following packages to be installed: + +- numpy: [installation instructions](https://numpy.org/install/) + +python -m pip install --user numpy + +- pytest: [installation instructions](https://docs.pytest.org/en/stable/getting-started/installation.html) + +python -m pip install --user pytest + +- scipy: [installation instructions](https://docs.scipy.org/doc/scipy/reference/building_scipy.html) + +python -m pip install --user scipy + +- matplotlib: [installation instructions](https://matplotlib.org/stable/users/installing.html) + +python -m pip install --user matplotlib ### From the command line -In `mole/python`, use the following command to install the MOLE module: +In the repository root, `mole`, use the following command to install the Python MOLE module: ```sh python -m pip install --user ./python @@ -22,7 +48,7 @@ python -m pip install --user ./python ## Running the test suite -To run the unit tests, first install the MOLE module as described above, then run the command `python -m pytest python/tests` from the directory `mole/python`). The results of the unit tests should be displayed to your console. +To run the unit tests, first install the MOLE module as described above, then run the command `python -m pytest python/tests` from the same root directory `mole`). The results of the unit tests should be displayed to your console. ## Examples @@ -32,4 +58,8 @@ Currently, the following examples are available in the MOLE Python package. - Elliptic Problems - 1D Examples - - `elliptic1D`: A script that solves the 1D Poisson's equation with Robin boundary conditions using mimetic operators. + - `elliptic1D`: A script that solves the 1D Poisson's equation with Robin boundary conditions using mimetic operators. Invoke as + + ```sh + python python/examples/elliptic1D.py + ``` From da91615e980eed0c2deec47fd235d65779ef2e45 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sat, 19 Sep 2026 23:27:38 -0700 Subject: [PATCH 47/66] Update python/pyproject.toml --- python/pyproject.toml | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/python/pyproject.toml b/python/pyproject.toml index b3ca50c3c..8502a1157 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -16,6 +16,16 @@ build-backend = "setuptools.build_meta" name = "pymole" version = "0.1.0" requires-python = ">=3.10" +dependencies = [ + "numpy", + "scipy", + "matplotlib", +] + +[project.optional-dependencies] +test = [ + "pytest", +] [tool.setuptools.packages.find] where = ["src"] From 8dc7d0f22959aaa1e914c33a2305d42b320c6b3b Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 20 Sep 2026 09:21:29 -0700 Subject: [PATCH 48/66] Delete python/LICENSE --- python/LICENSE | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 python/LICENSE diff --git a/python/LICENSE b/python/LICENSE deleted file mode 100644 index e69de29bb..000000000 From e01cee35bfba3a4eef03d841143f60f917fdd414 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 20 Sep 2026 10:00:54 -0700 Subject: [PATCH 49/66] Undo changes to examples/matlab_octave/elliptic1D.m --- examples/matlab_octave/elliptic1D.m | 6 ------ 1 file changed, 6 deletions(-) diff --git a/examples/matlab_octave/elliptic1D.m b/examples/matlab_octave/elliptic1D.m index 023b3750a..b0ef05ac9 100644 --- a/examples/matlab_octave/elliptic1D.m +++ b/examples/matlab_octave/elliptic1D.m @@ -14,17 +14,11 @@ L = lap(k, m, dx); % 1D Mimetic laplacian operator -full(L) -size(L) - % Impose Robin BC on laplacian operator a = 1; b = 1; L = L + robinBC(k, m, dx, a, b); -full(L) -%size(L) - % 1D Staggered grid grid = [west west+dx/2 : dx : east-dx/2 east]; From 2f8c020b62a342b681126992484de0631a59c53b Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 20 Sep 2026 10:13:41 -0700 Subject: [PATCH 50/66] Delete examples/python/simple_operators_check.py --- examples/python/simple_operators_check.py | 88 ----------------------- 1 file changed, 88 deletions(-) delete mode 100644 examples/python/simple_operators_check.py diff --git a/examples/python/simple_operators_check.py b/examples/python/simple_operators_check.py deleted file mode 100644 index caa8ec5f9..000000000 --- a/examples/python/simple_operators_check.py +++ /dev/null @@ -1,88 +0,0 @@ -from pathlib import Path -import sys - -sys.path.insert(0, str(Path(__file__).resolve().parent.parent.parent / "src" / "python")) - -import numpy as np -from numpy.typing import NDArray -from scipy import sparse - -from pymole import ( - Divergence1D, - Divergence2D, - Gradient1D, - Gradient2D, - Grid1D, - Laplacian1D, - Laplacian2D, -) - - -def report_grid(name: str, grid: Grid1D) -> None: - print(f"{name}: {grid}; cell size (resolution): {grid.cell_size:.4f}") - - -def report_nullity_check( - name: str, - op: sparse.csr_matrix, - expected_shape: tuple[int, int], - constant_field: NDArray[np.float64], -) -> None: - print(f"{name} shape: {op.shape}; expected {expected_shape}") - print( - f"||{name} @ constant||_inf = {np.linalg.norm(op @ constant_field, ord=np.inf)}" - ) - - -def check_1d_operators(accuracy_order: int = 2) -> None: - num_cells = 2 * accuracy_order + 1 - - grid = Grid1D(0.0, 1.0, num_cells) - - G = Gradient1D(grid, accuracy_order).matrix - D = Divergence1D(grid, accuracy_order).matrix - L = Laplacian1D(grid, accuracy_order).matrix - - f_constant = np.ones(num_cells + 2) - v_constant = np.ones(num_cells + 1) - - print("=== 1D operators ===") - report_grid("x grid", grid) - report_nullity_check("G", G, (num_cells + 1, num_cells + 2), f_constant) - report_nullity_check("D", D, (num_cells + 2, num_cells + 1), v_constant) - report_nullity_check("L", L, (num_cells + 2, num_cells + 2), f_constant) - print(f"D @ G == L? -> {abs(D @ G - L).max() < 1e-12}") - print() - - -def check_2d_operators(accuracy_order: int = 2) -> None: - x_grid = Grid1D(-1.0, 1.0, 2 * accuracy_order + 1) - y_grid = Grid1D(-0.5, 2.0, 2 * accuracy_order + 2) - - m = x_grid.num_cells - n = y_grid.num_cells - - G = Gradient2D(x_grid, y_grid, accuracy_order).matrix - D = Divergence2D(x_grid, y_grid, accuracy_order).matrix - L = Laplacian2D(x_grid, y_grid, accuracy_order).matrix - - f_constant = np.ones((m + 2) * (n + 2)) - v_constant = np.ones(2 * m * n + m + n) - - print("=== 2D operators ===") - report_grid("x grid", x_grid) - report_grid("y grid", y_grid) - report_nullity_check("G", G, (2 * m * n + m + n, (m + 2) * (n + 2)), f_constant) - report_nullity_check("D", D, ((m + 2) * (n + 2), 2 * m * n + m + n), v_constant) - report_nullity_check("L", L, ((m + 2) * (n + 2), (m + 2) * (n + 2)), f_constant) - print(f"D @ G == L? -> {abs(D @ G - L).max() < 1e-12}") - print() - - -def main() -> None: - check_1d_operators() - check_2d_operators() - - -if __name__ == "__main__": - main() From c58a93482a11c6cf65d08ed620b1e312292f203a Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 20 Sep 2026 11:31:18 -0700 Subject: [PATCH 51/66] feat: add core pymole python package modules and mimetic operators --- python/examples/elliptic1D.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/python/examples/elliptic1D.py b/python/examples/elliptic1D.py index 8a70a9f99..d7b864d49 100644 --- a/python/examples/elliptic1D.py +++ b/python/examples/elliptic1D.py @@ -8,9 +8,6 @@ import shutil import matplotlib.pyplot as plt -from pathlib import Path - -sys.path.insert(0, str(Path(__file__).resolve().parent.parent.parent)+'/src/python') from pymole import Grid, Gradient, Divergence, Laplacian, RobinBoundaryCondition west = 0 # Domain's limits From 82167235b8f201d1ce5eba3348ccff9810c32ec8 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 20 Sep 2026 16:29:11 -0700 Subject: [PATCH 52/66] feat: add pymole Python package with grid, boundary condition, and mimetic operator modules --- python/pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/pyproject.toml b/python/pyproject.toml index 8502a1157..68187b0fb 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -14,7 +14,7 @@ build-backend = "setuptools.build_meta" [project] name = "pymole" -version = "0.1.0" +version = "1.2.0" requires-python = ">=3.10" dependencies = [ "numpy", From ec50ffcb5dcf1efdd6f9ebf1c02e90a40b6ee096 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 20 Sep 2026 16:36:28 -0700 Subject: [PATCH 53/66] feat: add core mimetic operators, grid, and boundary condition modules for pymole --- python/build/lib/pymole/BoundaryCondition.py | 71 ++++++++++++++++++++ 1 file changed, 71 insertions(+) create mode 100644 python/build/lib/pymole/BoundaryCondition.py diff --git a/python/build/lib/pymole/BoundaryCondition.py b/python/build/lib/pymole/BoundaryCondition.py new file mode 100644 index 000000000..073d94ad4 --- /dev/null +++ b/python/build/lib/pymole/BoundaryCondition.py @@ -0,0 +1,71 @@ +from abc import ABC, abstractmethod + +import numpy as np +from scipy import sparse + +from .gradient import Gradient + +class BoundaryCondition(ABC): + @abstractmethod + def apply(self): + """Subclasses must implement this.""" + raise NotImplementedError + + def isImplemented(self) -> bool: + """Returns True if the boundary condition is implemented.""" + return self.apply != BoundaryCondition.apply + + @classmethod + def hasImplemented(cls, bc) -> bool: + """Returns True if the boundary condition is implemented.""" + return bc.apply != BoundaryCondition.apply + + def __init__(self, grid, accuracy_order: int = 2): + self.grid = grid + self.accuracy_order = accuracy_order + + if self.grid.ndim not in (1, 2, 3): + raise ValueError( + f"BoundaryCondition currently supports 1D, 2D, and 3D grids. " + f"Got ndim={self.grid.ndim}" + ) + + self._configure_dimensions() + self._matrix = None + + def _configure_dimensions(self) -> None: + if self.grid.ndim == 1: + self.m = self.grid.num_cells + self.n = None + self.o = None + elif self.grid.ndim == 2: + self.m, self.n = self.grid.num_cells + self.o = None + else: + self.m, self.n, self.o = self.grid.num_cells + + @property + def matrix(self): + if self._matrix is None: + self._matrix = self.apply() + return self._matrix + + def _spacing_for_dimension(self, axis: int) -> float: + spacing = self.grid.spacing + if self.grid.ndim == 1: + return spacing + if self.grid.ndim == 2: + return spacing[axis] + return spacing[axis] + + def _build_1d_robin_matrix(self, size: int, dx: float, a: float, b: float): + A = sparse.lil_matrix((size + 2, size + 2), dtype=float) + A[0, 0] = a + A[-1, -1] = a + + B = sparse.lil_matrix((size + 2, size + 1), dtype=float) + B[0, 0] = -b + B[-1, -1] = b + + G = Gradient(self.grid, self.accuracy_order).matrix + return (A + B @ G).tocsr() \ No newline at end of file From 59d971958ec02421907f97cbe4b3ec0ba40fc16e Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 20 Sep 2026 16:56:33 -0700 Subject: [PATCH 54/66] feat: add pymole Python package with mimetic operators and grid classes --- python/examples/elliptic1D.py | 45 ++++++++++++--------- python/src/pymole/RobinBoundaryCondition.py | 9 ++++- 2 files changed, 34 insertions(+), 20 deletions(-) diff --git a/python/examples/elliptic1D.py b/python/examples/elliptic1D.py index d7b864d49..91f509eca 100644 --- a/python/examples/elliptic1D.py +++ b/python/examples/elliptic1D.py @@ -1,8 +1,9 @@ -''' - SPDX-License-Identifier: GPL-3.0-or-later - © 2008-2024 San Diego State University Research Foundation (SDSURF). - See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. -''' +""" +SPDX-License-Identifier: GPL-3.0-or-later +© 2008-2024 San Diego State University Research Foundation (SDSURF). +See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +""" + import numpy as np import sys import shutil @@ -10,13 +11,13 @@ from pymole import Grid, Gradient, Divergence, Laplacian, RobinBoundaryCondition -west = 0 # Domain's limits +west = 0 # Domain's limits east = 1 -k = 6 # Operator's order of accuracy -m = 2*k+1 # Minimum number of cells to attain the desired accuracy +k = 6 # Operator's order of accuracy +m = 2 * k + 1 # Minimum number of cells to attain the desired accuracy -grid = Grid.generate(west, east, shape=m+1) +grid = Grid.generate(west, east, shape=m + 1) L = Laplacian(grid, accuracy_order=k).matrix @@ -26,16 +27,24 @@ BC = RobinBoundaryCondition(grid, k, dirichlet_coefficient=1.0, neumann_coefficient=1.0) L_BC = L + BC.matrix -x = np.r_[west, np.arange(west + grid.spacing/2, east - grid.spacing/2 + grid.spacing/2, grid.spacing), east] +x = np.r_[ + west, + np.arange( + west + grid.spacing / 2, + east - grid.spacing / 2 + grid.spacing / 2, + grid.spacing, + ), + east, +] U = np.exp(x)[:, None] -U[0] = 0 # West BC -U[-1] = 2*np.exp(1) # East BC +U[0] = 0 # West BC +U[-1] = 2 * np.exp(1) # East BC U = np.linalg.inv(L_BC.toarray()) @ U -plt.plot(x, U, 'o', label='Approximated') -plt.plot(x, np.exp(x), label='Analytical') -plt.legend(loc='upper left') +plt.plot(x, U, "o", label="Approximated") +plt.plot(x, np.exp(x), label="Analytical") +plt.legend(loc="upper left") plt.title("Poisson's equation with Robin BC") -plt.xlabel('x') -plt.ylabel('u(x)') -plt.show() \ No newline at end of file +plt.xlabel("x") +plt.ylabel("u(x)") +plt.show() diff --git a/python/src/pymole/RobinBoundaryCondition.py b/python/src/pymole/RobinBoundaryCondition.py index 852b3e1d8..a78bfe6b6 100644 --- a/python/src/pymole/RobinBoundaryCondition.py +++ b/python/src/pymole/RobinBoundaryCondition.py @@ -8,9 +8,14 @@ from .BoundaryCondition import BoundaryCondition - class RobinBoundaryCondition(BoundaryCondition): - def __init__(self, grid, accuracy_order: int = 2, dirichlet_coefficient: float = 1.0, neumann_coefficient: float = 1.0): + def __init__( + self, + grid, + accuracy_order: int = 2, + dirichlet_coefficient: float = 1.0, + neumann_coefficient: float = 1.0, + ): self.dirichlet_coefficient = dirichlet_coefficient self.neumann_coefficient = neumann_coefficient super().__init__(grid, accuracy_order) From 6178fb4c474832159ff5bf58f4e047b41f780cc0 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 20 Sep 2026 19:22:14 -0700 Subject: [PATCH 55/66] feat: add pymole Python package with mimetic operators and grid classes --- .github/scripts/filter_compiled_sources.py | 42 ++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 .github/scripts/filter_compiled_sources.py diff --git a/.github/scripts/filter_compiled_sources.py b/.github/scripts/filter_compiled_sources.py new file mode 100644 index 000000000..23b329b06 --- /dev/null +++ b/.github/scripts/filter_compiled_sources.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 +""" +Filter a list of candidate C/C++ source files down to only those that +actually appear in compile_commands.json (i.e. files CMake really +compiles). This automatically excludes old_version/ and any other +directory not wired into the CMake build, with no hardcoded paths. + +Usage: + python3 filter_compiled_sources.py +""" +import json +import os +import sys + + +def main(): + if len(sys.argv) != 4: + print("Usage: filter_compiled_sources.py ") + sys.exit(2) + + compile_commands_path, raw_files_path, filtered_files_path = sys.argv[1:4] + + with open(compile_commands_path) as f: + compiled = {os.path.realpath(entry["file"]) for entry in json.load(f)} + + with open(raw_files_path) as f: + candidates = [line.strip() for line in f if line.strip()] + + kept = [c for c in candidates if os.path.realpath(c) in compiled] + skipped = [c for c in candidates if os.path.realpath(c) not in compiled] + + if skipped: + print(f"Skipping {len(skipped)} file(s) not in compile_commands.json:") + for s in skipped: + print(f" - {s}") + + with open(filtered_files_path, "w") as f: + f.write("\n".join(kept)) + + +if __name__ == "__main__": + main() From b051f2cf84abd87c0f198b28cc2e4dbe276694ad Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 22 Sep 2026 19:19:18 -0700 Subject: [PATCH 56/66] feat: add core pymole package modules including Grid, Curl, and RobinBoundaryCondition --- python/build/lib/pymole/__init__.py | 19 +++++++++++++++++++ ...daryCondition.py => boundary_condition.py} | 0 ...ndition.py => robin_boundary_condition.py} | 0 3 files changed, 19 insertions(+) create mode 100644 python/build/lib/pymole/__init__.py rename python/src/pymole/{BoundaryCondition.py => boundary_condition.py} (100%) rename python/src/pymole/{RobinBoundaryCondition.py => robin_boundary_condition.py} (100%) diff --git a/python/build/lib/pymole/__init__.py b/python/build/lib/pymole/__init__.py new file mode 100644 index 000000000..e41afd6d0 --- /dev/null +++ b/python/build/lib/pymole/__init__.py @@ -0,0 +1,19 @@ +from .grid import Grid +from .curl import Curl +from .divergence import Divergence +from .gradient import Gradient +from .laplacian import Laplacian +from .BoundaryCondition import BoundaryCondition +from .RobinBoundaryCondition import RobinBoundaryCondition + +__version__ = "0.1.1" +__all__ = [ + "BoundaryCondition", + "Curl", + "Divergence", + "Gradient", + "Grid", + "Laplacian", + "RobinBoundaryCondition", +] + diff --git a/python/src/pymole/BoundaryCondition.py b/python/src/pymole/boundary_condition.py similarity index 100% rename from python/src/pymole/BoundaryCondition.py rename to python/src/pymole/boundary_condition.py diff --git a/python/src/pymole/RobinBoundaryCondition.py b/python/src/pymole/robin_boundary_condition.py similarity index 100% rename from python/src/pymole/RobinBoundaryCondition.py rename to python/src/pymole/robin_boundary_condition.py From 12c5cd0f4888a66caa654f44a35b59ba6520bd70 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Tue, 22 Sep 2026 22:22:33 -0700 Subject: [PATCH 57/66] feat: add core pymole python package modules for grids, boundary conditions, and mimetic operators --- python/examples/model_problem.py | 64 ++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 python/examples/model_problem.py diff --git a/python/examples/model_problem.py b/python/examples/model_problem.py new file mode 100644 index 000000000..32764b2c2 --- /dev/null +++ b/python/examples/model_problem.py @@ -0,0 +1,64 @@ +""" +SPDX-License-Identifier: GPL-3.0-or-later +© 2008-2024 San Diego State University Research Foundation (SDSURF). +See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +""" + +import numpy as np +import sys +import shutil +import matplotlib.pyplot as plt + +from pymole import Grid, Gradient, Divergence, Laplacian, RobinBoundaryCondition + +west = 0 # Domain's limits +east = 1 + +k = 6 # Operator's order of accuracy +m = 2 * k + 1 # Minimum number of cells to attain the desired accuracy +lambda_val = -1.0 +alpha = -np.exp(lambda_val) +beta = (np.exp(lambda_val)-1)/lambda_val + +grid = Grid.generate(west, east, shape=m + 1) +print(grid.spacing) +L = Laplacian(grid, accuracy_order=k).matrix + +terminal_width, _ = shutil.get_terminal_size() +np.set_printoptions(linewidth=terminal_width, suppress=True) + +B = RobinBoundaryCondition(grid, k, dirichlet_coefficient=alpha, neumann_coefficient=beta) +L_BC = L + B.matrix + +print(L_BC.toarray()) + + +x = np.r_[ + west, + np.arange( + west + grid.spacing / 2, + east - grid.spacing / 2 + grid.spacing / 2, + grid.spacing, + ), + east, +] + + +F = (-lambda_val * np.exp(lambda_val*x))/(np.exp(lambda_val)-1) +U = np.linalg.inv(L_BC.toarray()) @ -F + +""" + + +U[0] = 0 # West BC +U[-1] = 2 * np.exp(1) # East BC +U = np.linalg.inv(L_BC.toarray()) @ U +""" + +plt.plot(x, U, "o", label="Approximated") +#plt.plot(x, np.exp(x), label="Analytical") +plt.legend(loc="upper left") +plt.title("The Model Problem") +plt.xlabel("x") +plt.ylabel("u(x)") +plt.show() From 20e40701601e07169e6ec8d643b66080f8d3b8df Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 23 Sep 2026 14:37:00 -0700 Subject: [PATCH 58/66] feat: add RobinBoundaryCondition class for multi-dimensional grid support --- python/.gitkeep | 0 python/docs/.gitkeep | 0 python/examples/.gitkeep | 0 python/src/.gitkeep | 0 4 files changed, 0 insertions(+), 0 deletions(-) delete mode 100644 python/.gitkeep delete mode 100644 python/docs/.gitkeep delete mode 100644 python/examples/.gitkeep delete mode 100644 python/src/.gitkeep diff --git a/python/.gitkeep b/python/.gitkeep deleted file mode 100644 index e69de29bb..000000000 diff --git a/python/docs/.gitkeep b/python/docs/.gitkeep deleted file mode 100644 index e69de29bb..000000000 diff --git a/python/examples/.gitkeep b/python/examples/.gitkeep deleted file mode 100644 index e69de29bb..000000000 diff --git a/python/src/.gitkeep b/python/src/.gitkeep deleted file mode 100644 index e69de29bb..000000000 From 49b638ee5eec3f3a93d7e99e3e56cd17d6a48cb7 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 23 Sep 2026 14:39:13 -0700 Subject: [PATCH 59/66] chore: remove file --- python/docs/_build/html/.gitkeep | 1 - python/docs/_static/.gitkeep | 1 - 2 files changed, 2 deletions(-) delete mode 100644 python/docs/_build/html/.gitkeep delete mode 100644 python/docs/_static/.gitkeep diff --git a/python/docs/_build/html/.gitkeep b/python/docs/_build/html/.gitkeep deleted file mode 100644 index 8b1378917..000000000 --- a/python/docs/_build/html/.gitkeep +++ /dev/null @@ -1 +0,0 @@ - diff --git a/python/docs/_static/.gitkeep b/python/docs/_static/.gitkeep deleted file mode 100644 index 8b1378917..000000000 --- a/python/docs/_static/.gitkeep +++ /dev/null @@ -1 +0,0 @@ - From 3dcfc739da46149ba0c2d058fe351c7329016c5f Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 23 Sep 2026 14:52:56 -0700 Subject: [PATCH 60/66] feat: add grid, mimetic operators, and boundary condition modules to pymole package --- python/build/lib/pymole/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/build/lib/pymole/__init__.py b/python/build/lib/pymole/__init__.py index e41afd6d0..e43c001b9 100644 --- a/python/build/lib/pymole/__init__.py +++ b/python/build/lib/pymole/__init__.py @@ -3,8 +3,8 @@ from .divergence import Divergence from .gradient import Gradient from .laplacian import Laplacian -from .BoundaryCondition import BoundaryCondition -from .RobinBoundaryCondition import RobinBoundaryCondition +from .boundary_condition import BoundaryCondition +from .robin_boundary_condition import RobinBoundaryCondition __version__ = "0.1.1" __all__ = [ From 6a3a9f559640484cc245f12480a3c5a0f63e622d Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 23 Sep 2026 14:53:26 -0700 Subject: [PATCH 61/66] feat: add pymole python package with mimetic operators, boundary conditions, and grid implementations --- python/src/pymole/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/src/pymole/__init__.py b/python/src/pymole/__init__.py index e41afd6d0..e43c001b9 100644 --- a/python/src/pymole/__init__.py +++ b/python/src/pymole/__init__.py @@ -3,8 +3,8 @@ from .divergence import Divergence from .gradient import Gradient from .laplacian import Laplacian -from .BoundaryCondition import BoundaryCondition -from .RobinBoundaryCondition import RobinBoundaryCondition +from .boundary_condition import BoundaryCondition +from .robin_boundary_condition import RobinBoundaryCondition __version__ = "0.1.1" __all__ = [ From 115b3984529335a80870d93b5ad884513556e0b7 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 23 Sep 2026 18:14:18 -0700 Subject: [PATCH 62/66] feat: add pymole python package for mimetic operators, grids, and boundary conditions with an example model problem --- python/examples/model_problem.py | 93 ++++++++++++++++++++++++++------ 1 file changed, 77 insertions(+), 16 deletions(-) diff --git a/python/examples/model_problem.py b/python/examples/model_problem.py index 32764b2c2..d373ce6a3 100644 --- a/python/examples/model_problem.py +++ b/python/examples/model_problem.py @@ -16,9 +16,52 @@ k = 6 # Operator's order of accuracy m = 2 * k + 1 # Minimum number of cells to attain the desired accuracy -lambda_val = -1.0 -alpha = -np.exp(lambda_val) -beta = (np.exp(lambda_val)-1)/lambda_val +lam = -1.0 +alpha = -np.exp(lam) +beta = (np.exp(lam)-1)/lam + +# x grid +# +n = m - 2 +if n <= 0: + raise ValueError("n must be positive") + +h = 1.0 / n + +# Cell centers: h/2, 3h/2, ..., 1-h/2 +x_internal = (np.arange(n) + 0.5) * h + +# Include the two boundary points +x = np.concatenate(([0.0], x_internal, [1.0])) + +# Exact solution +# +exp_lam = np.exp(lam) +expm1_lam = np.expm1(lam) # exp(lam) - 1 +C1 = (lam - 1.0) / (2.0 * expm1_lam - lam * exp_lam) +C2 = ( + lam + - 1.0 + - exp_lam / expm1_lam + - expm1_lam * C1 + ) / (lam * exp_lam) +f = ( + np.exp(lam * x) / (lam * expm1_lam) + + C1 * x + + C2 + ) + +b = np.concatenate(( + [-1.0], + -lam * np.exp(lam * x) / np.expm1(lam), + [0.0] +)) + +#F = (-lambda_val * np.exp(lambda_val*x))/(np.exp(lambda_val)-1) + +print(b.shape) +#print(L_BC.toarray().shape) + grid = Grid.generate(west, east, shape=m + 1) print(grid.spacing) @@ -32,28 +75,45 @@ print(L_BC.toarray()) +#RHS = [-1,0] -x = np.r_[ - west, - np.arange( - west + grid.spacing / 2, - east - grid.spacing / 2 + grid.spacing / 2, - grid.spacing, - ), - east, -] +#x = np.r_[ +# RHS[0], +# np.arange( +# RHS[0] + grid.spacing / 2, +# RHS[1] + grid.spacing / 2 + grid.spacing / 2, +# grid.spacing, +# ), +# RHS[1], +#] +#x_internal = (np.arange(m) + 0.5) * grid.spacing +#x = np.r_[ +# RHS[0], +# x_internal, +# RHS[1], +#] +#print(RHS[0] + grid.spacing / 2) +#print(RHS[1] + grid.spacing / 2 + grid.spacing / 2) +#print(grid.spacing) -F = (-lambda_val * np.exp(lambda_val*x))/(np.exp(lambda_val)-1) -U = np.linalg.inv(L_BC.toarray()) @ -F +# +#print(x) +#print(x.shape) -""" +U = np.linalg.inv(L_BC.toarray()) @ -b + + +print(U.shape) +print(x.shape) +""" + U[0] = 0 # West BC U[-1] = 2 * np.exp(1) # East BC U = np.linalg.inv(L_BC.toarray()) @ U -""" + plt.plot(x, U, "o", label="Approximated") #plt.plot(x, np.exp(x), label="Analytical") @@ -62,3 +122,4 @@ plt.xlabel("x") plt.ylabel("u(x)") plt.show() +#""" \ No newline at end of file From 131e3291b4d0ebf7c731c5ae53d3198f6c8997f8 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Wed, 23 Sep 2026 21:50:35 -0700 Subject: [PATCH 63/66] feat: add Python modules for pymole mimetic operators, grids, and boundary conditions --- python/examples/model_problem.py | 98 ++++++++++---------------------- 1 file changed, 29 insertions(+), 69 deletions(-) diff --git a/python/examples/model_problem.py b/python/examples/model_problem.py index d373ce6a3..f2aee208c 100644 --- a/python/examples/model_problem.py +++ b/python/examples/model_problem.py @@ -20,19 +20,26 @@ alpha = -np.exp(lam) beta = (np.exp(lam)-1)/lam -# x grid -# -n = m - 2 -if n <= 0: - raise ValueError("n must be positive") +# Generate the discrete operator first +grid = Grid.generate(west, east, shape=m + 1) +#print(grid.spacing) +L = Laplacian(grid, accuracy_order=k).matrix -h = 1.0 / n +B = RobinBoundaryCondition(grid, k, dirichlet_coefficient=alpha, neumann_coefficient=beta) +L_BC = -L + B.matrix +A = L_BC.toarray() + +# Number of interior cell centers required by the matrix +n = A.shape[0] - 2 +h = (east - west) / n +# Staggered grid +# # Cell centers: h/2, 3h/2, ..., 1-h/2 -x_internal = (np.arange(n) + 0.5) * h +x_internal = west + (np.arange(n) + 0.5) * h # Include the two boundary points -x = np.concatenate(([0.0], x_internal, [1.0])) +x = np.concatenate(([west], x_internal, [east])) # Exact solution # @@ -45,81 +52,34 @@ - exp_lam / expm1_lam - expm1_lam * C1 ) / (lam * exp_lam) -f = ( +f_exact = ( np.exp(lam * x) / (lam * expm1_lam) + C1 * x + C2 ) +F_internal = ( + -lam + * np.exp(lam * x_internal) + / np.expm1(lam) +) + b = np.concatenate(( [-1.0], - -lam * np.exp(lam * x) / np.expm1(lam), + F_internal, [0.0] )) -#F = (-lambda_val * np.exp(lambda_val*x))/(np.exp(lambda_val)-1) - -print(b.shape) -#print(L_BC.toarray().shape) - - -grid = Grid.generate(west, east, shape=m + 1) -print(grid.spacing) -L = Laplacian(grid, accuracy_order=k).matrix +# Solve without explicitly forming the matrix inverse +f_approx = np.linalg.solve(A, b) terminal_width, _ = shutil.get_terminal_size() np.set_printoptions(linewidth=terminal_width, suppress=True) -B = RobinBoundaryCondition(grid, k, dirichlet_coefficient=alpha, neumann_coefficient=beta) -L_BC = L + B.matrix - -print(L_BC.toarray()) - -#RHS = [-1,0] - -#x = np.r_[ -# RHS[0], -# np.arange( -# RHS[0] + grid.spacing / 2, -# RHS[1] + grid.spacing / 2 + grid.spacing / 2, -# grid.spacing, -# ), -# RHS[1], -#] -#x_internal = (np.arange(m) + 0.5) * grid.spacing -#x = np.r_[ -# RHS[0], -# x_internal, -# RHS[1], -#] - -#print(RHS[0] + grid.spacing / 2) -#print(RHS[1] + grid.spacing / 2 + grid.spacing / 2) -#print(grid.spacing) - -# -#print(x) -#print(x.shape) - - - -U = np.linalg.inv(L_BC.toarray()) @ -b - - -print(U.shape) -print(x.shape) -""" - -U[0] = 0 # West BC -U[-1] = 2 * np.exp(1) # East BC -U = np.linalg.inv(L_BC.toarray()) @ U - - -plt.plot(x, U, "o", label="Approximated") -#plt.plot(x, np.exp(x), label="Analytical") +plt.plot(x, f_approx, "o", label="Approximated") +plt.plot(x, f_exact, label="Analytical") plt.legend(loc="upper left") -plt.title("The Model Problem") +plt.title(fr"The Model Problem ($\lambda$ = {lam})") plt.xlabel("x") -plt.ylabel("u(x)") +plt.ylabel("f(x)") plt.show() -#""" \ No newline at end of file From d91a3298e6feb6ad99d32fd8520ef45277e8be23 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Thu, 24 Sep 2026 14:44:28 -0700 Subject: [PATCH 64/66] feat: add pymole python package with mimetic operators, boundary conditions, and examples --- python/tests/.gitkeep | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 python/tests/.gitkeep diff --git a/python/tests/.gitkeep b/python/tests/.gitkeep deleted file mode 100644 index e69de29bb..000000000 From 3027ed9f4945c4f966059d980482147f3015a145 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Thu, 24 Sep 2026 14:52:03 -0700 Subject: [PATCH 65/66] feat: add pymole Python library with mimetic operators, documentation, and examples --- python/docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/docs/conf.py b/python/docs/conf.py index 79138f95f..cbc6ad5bd 100644 --- a/python/docs/conf.py +++ b/python/docs/conf.py @@ -30,4 +30,4 @@ # -- Options for HTML output ------------------------------------------------ html_theme = "alabaster" -html_static_path = ["_static"] +#html_static_path = ["_static"] From 7e7f2f31ff241f600c1ac9d8aa4d12b98263a0b9 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Thu, 24 Sep 2026 15:06:54 -0700 Subject: [PATCH 66/66] feat: add pymole library operators, boundary conditions, and model problem examples --- python/src/pymole/robin_boundary_condition.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/src/pymole/robin_boundary_condition.py b/python/src/pymole/robin_boundary_condition.py index a78bfe6b6..6788a7f0d 100644 --- a/python/src/pymole/robin_boundary_condition.py +++ b/python/src/pymole/robin_boundary_condition.py @@ -6,7 +6,7 @@ from scipy import sparse -from .BoundaryCondition import BoundaryCondition +from .boundary_condition import BoundaryCondition class RobinBoundaryCondition(BoundaryCondition): def __init__(