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() diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 140c760b7..a711f2b33 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -271,3 +271,42 @@ 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 numpy scipy + python -m pip install -e ./python + pwd + + - name: Run tests + run: | + pytest python/tests diff --git a/.github/workflows/python-docs.yml b/.github/workflows/python-docs.yml new file mode 100644 index 000000000..105849712 --- /dev/null +++ b/.github/workflows/python-docs.yml @@ -0,0 +1,56 @@ +name: Python Documentation + +on: + pull_request: + branches: + - '**' + push: + branches: + - '**' + 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: | + 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 + cd python + pwd + if [ -f docs/requirements.txt ]; then pip install -r docs/requirements.txt; fi + + - name: Build documentation + run: | + pwd + cd python + pwd + sphinx-build -W --keep-going -b html docs docs/_build/html + + - name: Deploy documentation + if: github.event_name == 'push' + uses: peaceiris/actions-gh-pages@v4 + with: + github_token: ${{ secrets.GITHUB_TOKEN }} + publish_dir: python/docs/_build/html diff --git a/python/README.md b/python/README.md new file mode 100644 index 000000000..bb83ebb76 --- /dev/null +++ b/python/README.md @@ -0,0 +1,65 @@ +# 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. + +## 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. 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 the repository root, `mole`, use the following command to install the Python 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 same root directory `mole`). 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. Invoke as + + ```sh + python python/examples/elliptic1D.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 diff --git a/python/build/lib/pymole/__init__.py b/python/build/lib/pymole/__init__.py new file mode 100644 index 000000000..e43c001b9 --- /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 .boundary_condition import BoundaryCondition +from .robin_boundary_condition import RobinBoundaryCondition + +__version__ = "0.1.1" +__all__ = [ + "BoundaryCondition", + "Curl", + "Divergence", + "Gradient", + "Grid", + "Laplacian", + "RobinBoundaryCondition", +] + diff --git a/python/docs/conf.py b/python/docs/conf.py new file mode 100644 index 000000000..cbc6ad5bd --- /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"] diff --git a/python/docs/index.rst b/python/docs/index.rst new file mode 100644 index 000000000..e49e0cf57 --- /dev/null +++ b/python/docs/index.rst @@ -0,0 +1,8 @@ +Documentation for pymole, the Python MOLE module +================================================ + +Welcome to the pymole documentation. + +.. toctree:: + :maxdepth: 2 + :caption: Contents: 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 diff --git a/python/examples/elliptic1D.py b/python/examples/elliptic1D.py new file mode 100644 index 000000000..91f509eca --- /dev/null +++ b/python/examples/elliptic1D.py @@ -0,0 +1,50 @@ +""" +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 + +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) + +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, +] +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() diff --git a/python/examples/model_problem.py b/python/examples/model_problem.py new file mode 100644 index 000000000..f2aee208c --- /dev/null +++ b/python/examples/model_problem.py @@ -0,0 +1,85 @@ +""" +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 +lam = -1.0 +alpha = -np.exp(lam) +beta = (np.exp(lam)-1)/lam + +# Generate the discrete operator first +grid = Grid.generate(west, east, shape=m + 1) +#print(grid.spacing) +L = Laplacian(grid, accuracy_order=k).matrix + +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 = west + (np.arange(n) + 0.5) * h + +# Include the two boundary points +x = np.concatenate(([west], x_internal, [east])) + +# 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_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], + F_internal, + [0.0] +)) + +# 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) + +plt.plot(x, f_approx, "o", label="Approximated") +plt.plot(x, f_exact, label="Analytical") +plt.legend(loc="upper left") +plt.title(fr"The Model Problem ($\lambda$ = {lam})") +plt.xlabel("x") +plt.ylabel("f(x)") +plt.show() diff --git a/python/pyproject.toml b/python/pyproject.toml new file mode 100644 index 000000000..68187b0fb --- /dev/null +++ b/python/pyproject.toml @@ -0,0 +1,31 @@ +[tool.ruff] +src = ["src"] +fix = true + +[tool.ruff.lint] +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 = "1.2.0" +requires-python = ">=3.10" +dependencies = [ + "numpy", + "scipy", + "matplotlib", +] + +[project.optional-dependencies] +test = [ + "pytest", +] + +[tool.setuptools.packages.find] +where = ["src"] diff --git a/python/src/pymole/__init__.py b/python/src/pymole/__init__.py new file mode 100644 index 000000000..e43c001b9 --- /dev/null +++ b/python/src/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 .boundary_condition import BoundaryCondition +from .robin_boundary_condition import RobinBoundaryCondition + +__version__ = "0.1.1" +__all__ = [ + "BoundaryCondition", + "Curl", + "Divergence", + "Gradient", + "Grid", + "Laplacian", + "RobinBoundaryCondition", +] + diff --git a/python/src/pymole/boundary_condition.py b/python/src/pymole/boundary_condition.py new file mode 100644 index 000000000..073d94ad4 --- /dev/null +++ b/python/src/pymole/boundary_condition.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 diff --git a/python/src/pymole/curl.py b/python/src/pymole/curl.py new file mode 100644 index 000000000..fff4d8c44 --- /dev/null +++ b/python/src/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/python/src/pymole/divergence.py b/python/src/pymole/divergence.py new file mode 100644 index 000000000..939c8b9be --- /dev/null +++ b/python/src/pymole/divergence.py @@ -0,0 +1,200 @@ +''' + 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 not in (2, 4, 6): + raise ValueError( + f"Unsupported order of accuracy: {self.accuracy_order}" + ) + + match self.grid.ndim: + + case 1: + 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 _: + raise ValueError( + f"Unsupported grid dimension: {self.grid.ndim}" + ) + + def _construct_1d(self): + + m = self.grid.num_cells + + # 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"points 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, + ): + + D = sparse.lil_matrix((m + 2, m + 1), dtype=float) + + 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) + + 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] + + case _: + raise ValueError( + f"Unsupported order of accuracy: {self.accuracy_order}" + ) + + return (D / dx).tocsr() + + 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/python/src/pymole/gradient.py b/python/src/pymole/gradient.py new file mode 100644 index 000000000..5d602a3b4 --- /dev/null +++ b/python/src/pymole/gradient.py @@ -0,0 +1,315 @@ +''' + 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 + + @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 not in (2, 4, 6): + raise ValueError( + f"Unsupported order of accuracy: {self.accuracy_order}" + ) + + match self.grid.ndim: + + case 1: + 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) + + case _: + raise ValueError( + f"Unsupported grid dimension: {self.grid.ndim}" + ) + + def _construct_1d(self): + """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( + f"Grid requires at least " + f"{self.minimum_cells_for_order(self.accuracy_order)} " + f"points for order {self.accuracy_order}" + ) + + m = self.grid.num_cells + return self._build_1d_gradient_mimetic(m, 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 from nodal scalar values.""" + + 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." + ) + + Gx = self._build_2d_gradient_x(m, n, dx) + Gy = self._build_2d_gradient_y(m, n, dy) + + 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).""" + + 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 square 1D central-difference gradient used by 3D helpers.""" + + 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.""" + + 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 G.tocsr() + + def _build_2d_gradient_y(self, m: int, n: int, dy: float): + """Builds 2D gradient operator in y-direction.""" + + 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 G.tocsr() + + 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/python/src/pymole/grid.py b/python/src/pymole/grid.py new file mode 100644 index 000000000..062f60598 --- /dev/null +++ b/python/src/pymole/grid.py @@ -0,0 +1,193 @@ +''' + 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, Optional + +class Grid: + """ + 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 + @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, 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): + """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, ...]], + p2: Union[float, int, np.ndarray, List[float], Tuple[float, ...]], + 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' + ) -> 'Grid': + """ + 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). + 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 Grid object encapsulating the coordinate arrays. + """ + 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) + + 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 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: + 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)] + + if dims == 1: + return Grid(coords[0]) + + 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 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] + 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] (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/python/src/pymole/laplacian.py b/python/src/pymole/laplacian.py new file mode 100644 index 000000000..c87ae6f13 --- /dev/null +++ b/python/src/pymole/laplacian.py @@ -0,0 +1,41 @@ +''' + 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 +from pymole.gradient import Gradient + +from pymole import Gradient +from pymole import Divergence + +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/python/src/pymole/mimetic_operator.py b/python/src/pymole/mimetic_operator.py new file mode 100644 index 000000000..465b33e0e --- /dev/null +++ b/python/src/pymole/mimetic_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/python/src/pymole/requirements.txt b/python/src/pymole/requirements.txt new file mode 100644 index 000000000..77909df50 --- /dev/null +++ b/python/src/pymole/requirements.txt @@ -0,0 +1,3 @@ +numpy +scipy +pytest diff --git a/python/src/pymole/robin_boundary_condition.py b/python/src/pymole/robin_boundary_condition.py new file mode 100644 index 000000000..6788a7f0d --- /dev/null +++ b/python/src/pymole/robin_boundary_condition.py @@ -0,0 +1,65 @@ +''' + 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 .boundary_condition 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}") diff --git a/python/tests/__init__.py b/python/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/python/tests/test_curl.py b/python/tests/test_curl.py new file mode 100644 index 000000000..5de405f1a --- /dev/null +++ b/python/tests/test_curl.py @@ -0,0 +1,31 @@ +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/python/tests/test_divergence.py b/python/tests/test_divergence.py new file mode 100644 index 000000000..8aa5db2bf --- /dev/null +++ b/python/tests/test_divergence.py @@ -0,0 +1,168 @@ +''' + 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, Divergence + +class TestDivergence: + + def test_unsupported_accuracy_order(self): + + grid = Grid.generate( + 0, + 10, + shape=10, + ) + + with pytest.raises(ValueError): + Divergence( + grid, + accuracy_order=6, + ) + + def test_grid_too_small_1d(self): + + grid = Grid.generate( + 0, + 10, + shape=5, + ) + + with pytest.raises(ValueError): + Divergence( + grid, + 4, + ) + + 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, + ) diff --git a/python/tests/test_gradient.py b/python/tests/test_gradient.py new file mode 100644 index 000000000..b81544295 --- /dev/null +++ b/python/tests/test_gradient.py @@ -0,0 +1,283 @@ +''' + 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=8, + ) + + 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+1, + ) + + 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.num_cells + total_nodes = (m + 1) * (n + 1) + + assert Gx.shape == (n * (m + 1), total_nodes) + assert Gy.shape == (m * (n + 1), 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, + ).matrix + + n = len(grid.x) + + u = np.ones(n+1) + + 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_gradient_accuracy(self): + + grid = Grid.generate( + 0.0, + 10.0, + shape=101, + ) + + x = np.linspace(0, 10, 101) + xSq = x ** 2 + dx = x[1] - x[0] + expected_grad_xSq = np.gradient(xSq, dx) + + Gx = Gradient( + grid, + 2, + ).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], + atol=1e-1, + ) + + 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) diff --git a/python/tests/test_grid_utility.py b/python/tests/test_grid_utility.py new file mode 100644 index 000000000..12bbf9ee8 --- /dev/null +++ b/python/tests/test_grid_utility.py @@ -0,0 +1,163 @@ +''' + 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 pathlib import Path +import pytest +import numpy as np +import sys + +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) + 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.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.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.x, expected) + assert grid.x.shape == (4,) + +def test_shape_grid_2d(): + # Test 2D grid generation using shape + 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) + + # 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 + + # 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 + 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) + + 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 + 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 + 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) + +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]) diff --git a/python/tests/test_laplacian.py b/python/tests/test_laplacian.py new file mode 100644 index 000000000..d86977a52 --- /dev/null +++ b/python/tests/test_laplacian.py @@ -0,0 +1,94 @@ +''' + 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, Divergence, 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 + 1) + + 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) + 1) + 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 + x = np.append(x, x[-1] + grid.spacing) + 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() + + 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)