diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 5eb1b20a7..60885bc19 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -23,9 +23,11 @@ permissions: jobs: build_docs: runs-on: ubuntu-latest + env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN}} OMP_NUM_THREADS: 2 + steps: - name: Checkout uses: actions/checkout@v6 @@ -38,17 +40,7 @@ jobs: python-version: '3.10' - name: Install non-Python dependencies on Ubuntu - uses: awalsh128/cache-apt-pkgs-action@latest - with: - packages: gfortran openmpi-bin libopenmpi-dev libhdf5-openmpi-dev - version: 1.0 - execute_install_scripts: true - - - name: Reconfigure non-Python dependencies on Ubuntu - run: | - sudo apt-get update - sudo apt-get install --reinstall openmpi-bin libhdf5-openmpi-dev liblapack-dev libblas-dev - sudo apt install graphviz pandoc + uses: ./.github/actions/ubuntu_installations - name: Print information on MPI and HDF5 libraries run: | @@ -70,6 +62,10 @@ jobs: pip install .[test] pip freeze + - name: Install non-Python dependencies for Documentation + run: | + sudo apt install graphviz pandoc + - name: Install Python dependencies for Documentation run: | pip install -r docs/requirements.txt diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 1b4c2ff7b..ffd83e900 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -40,7 +40,8 @@ jobs: OMP_NUM_THREADS: 2 steps: - - uses: actions/checkout@v6 + - name: Checkout + uses: actions/checkout@v6 with: submodules: true diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 000000000..3390dfeed --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,104 @@ +# Project: PSYDAC + +## Overview +Python 3 library for isogeometric analysis (IGA). Solve general systems of partial +differential equations (PDEs) in weak form, defined using the domain-specific +language provided by SymPDE. Supports finite element exterior calculus (FEEC) +with tensor-product spline spaces. Handles multi-patch geometries in various ways; +usually broken-FEEC a.k.a. CONGA (conforming/non-conforming Galerkin). +Python code is automatically generated for the assembly of user-defined functionals, +linear forms, and bilinear forms. This Python code is then accelerated to C/Fortran +speed using Pyccel. The library enables large parallel computations on distributed- +memory supercomputers using MPI and OpenMP. + +## Tech stack +- Python 3.10+, type hints required in new code +- meson-python (Meson build backend) for building +- pip with submodules (for igakit) +- pytest + pytest-cov + pytest-mpi + pytest-xdist for testing +- sympde for symbolic definition of weak formulations +- pyccel for translation of Python kernels to Fortran or C +- mpi4py for MPI parallelization +- h5py for parallel I/O +- petsc4py for direct linear solvers (optional install) +- sphinx + numpydoc + nbsphinx + myst-parser for documentation + (see `docs/requirements.txt`; build with `make -C docs html`) + +## Code standards +- Follow PEP 8 style guide whenever possible +- Docstrings of public functions and classes follow Numpydoc conventions +- pylint + black + isort for linting and formatting +- Computational kernels to be accelerated with Pyccel named as `_kernels.py` +- Minimize code duplication. Never add code that already exists +- Strive for concise, clean, and human-readable code. Avoid useless verbosity +- Self-explanatory names for variables, functions, and classes +- Do not reassign value to an existing variable, especially if the type changes +- Use short comments (one-liners or inline) to explain "why" rather than "what" or "how" +- In docstrings, use single backticks for names and short code, + and double backticks only for code expressions +- Write new code in clean style (e.g. no superfluous parentheses, f-strings), + even if the surrounding legacy code does not +- Do not add author credits in code or comments; + credit goes in commit messages and PR descriptions + +## Testing conventions +- Each library subpackage contains a `tests/` folder with an `__init__.py` file +- Unit tests are part of the library and shipped with it +- Unit tests can be run with `psydac test` CLI command (see `README.md`) +- Parallel tests are marked with `@pytest.mark.mpi` (from pytest-mpi) and run with + `psydac test --mpi`; serial tests run with `psydac test` +- Test file names follow `test_.py` +- Test function names follow `test_` or `test__` +- Parametrize unit tests with `@pytest.mark.parametrize` to minimize code duplication +- Keep run time of tests at a minimum +- Aim for 100% coverage on newly committed code +- Every bugfix comes with a regression test that fails before the fix +- Use distinct sizes in each direction (e.g. 3 x 5 cells, not 4 x 4), + so that swapped indices are detected +- Write files only to the `tmp_path` fixture (`mpi_tmp_path` in MPI tests), + never to the working directory +- After changing a module, also run the tests of the subpackages which use it + (e.g. `psydac/api/tests` after changing `psydac/cad`) + +## File structure +- psydac/api - high-level Python interface +- psydac/cad - computer-aided design (CAD) functionality +- psydac/cmd - command-line interface (CLI) commands +- psydac/core - splines functionality +- psydac/ddm - MPI decomposition of domain, vectors of spline coefficients, and matrices +- psydac/feec - finite element exterior calculus (FEEC) +- psydac/fem - "finite element method" middle-level Python interface (spaces and fields) +- psydac/linalg - "linear algebra" low-level Python interface (vectors, linear operators, iterative solvers) +- psydac/mapping - spline and NURBS mappings, and a gallery of discrete mappings +- psydac/polar - implementation of polar splines for H1 spaces +- psydac/pyccel - distillation of an old version of Pyccel for Python code generation (to be removed) +- psydac/utilities - various generic functions used across the library +- docs - Sphinx documentation +- examples - Various examples of library usage (mostly Jupyter notebooks) +- scripts - Various scripts for library maintenance +- subprojects - Git submodules + +## Workflow +- PRs are squash-merged and the PR description becomes the commit message: + keep it complete, and write "(fixes #N)" next to each fix +- Every PR adds entries to `CHANGELOG.md` under Added/Changed/Fixed/Removed +- Format of the CHANGELOG entries: + . `- # : `, followed by `, see #` if there is one + . `- [DEVELOPER] ` for changes which only concern developers +- Issues describe bugs as they appear on `devel`, with a reproducer that runs on `devel` +- Do not commit changes to submodule pointers (`subprojects/`) unless intended + +## Always do +- Before claiming done, run these checks on the Python files you changed: + . Run `python -m black --check $(git diff --name-only --diff-filter=d devel -- "*.py")` + . Run `python -m isort --check $(git diff --name-only --diff-filter=d devel -- "*.py")` + . Run `pylint --disable=all --enable=unused-import $(git diff --name-only --diff-filter=d devel -- "*.py")` +- Until the code base is reformatted, black and isort fail on most existing files: + do not reformat code outside your change. The pylint check must pass. +- Unless specified differently, assume process with `rank=0` in MPI communicator is root +- Use root MPI process to perform serial operations (e.g. open/save image) + +## Never do +- Commit secrets or `.env` files +- Add new dependencies without team discussion +- Use `eval` function together with `input` diff --git a/CHANGELOG.md b/CHANGELOG.md index 48039f532..cb92b760e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,10 @@ All notable changes to this project will be documented in this file. ## Unreleased ### Added + +- #527 : Add constructor `Geometry.from_file` with option `domain` to reuse an existing SymPDE domain, after checking that it matches the file +- #527 : Add function `get_available_mappings` and option `mpi_dims_mask` to `discrete_mapping` in `psydac.mapping.discrete_gallery` +- #527 : Add property `mpi_dims_mask` to `DomainDecomposition` - #596 : Add `topetsc` method to `psydac.linalg.basic.ZeroOperator` - #596 : Add `out` argument to `array_to_psydac` and `vec_topetsc` - #576 : Add module `psydac.feec.polar.conga_projections` with broken FEEC polar projections in 2D @@ -12,21 +16,42 @@ All notable changes to this project will be documented in this file. - #576 : Add module `psydac.utilities.parallel_utils` for parallel execution and gathering variable-length arrays - #576 : Add module `psydac.utilities.operators` with class `Laplacian` used in some old examples - #577 : Add an installation configuration option to choose the backend language +- #567 : Improve `psydac test` command (with several new features) +- #565 : Expand editable install info in `README.md` - [DEVELOPER] Create action `install_petsc4py` to install PETSc & `petsc4py` w/ complex support +- [DEVELOPER] Configure Pylint in `pyproject.toml` +- [DEVELOPER] Configure isort in `pyproject.toml`, with the profile of black +- [DEVELOPER] Add `AGENTS.md` with developer rules ### Fixed +- #527 : Require `sympde==0.20.0` which changes how multipatch interfaces are defined +- #527 : Fix `psydac.cad.elevate` and `refine` for NURBS mappings, see #608 +- #527 : Keep the NURBS weights in `psydac.cad.translate`, see #609 +- #527 : Fix type checks of input mappings in module `psydac.cad.cad`, see #610 +- #527 : Fix `psydac.cad.elevate` and `refine` for distributed mappings, see #611 +- #527 : Check the weights field in `NurbsMapping.__init__`, see #616 +- #527 : Fix reading multipatch geometries with `pdim=1`, see #617 +- #527 : Fix `psydac.cad.elevate` and `refine` for mappings with `pdim=1`, see #618 +- #527 : Fix `psydac.cad.refine` when inserting a value which is already a knot, see #619 +- #527 : Keep the `mpi_dims_mask` of a `DomainDecomposition` in `TensorFemSpace.add_refined_space` and `psydac.cad.cad.refine`, see #622 +- #527 : Fix `local_ncells` of the decomposition returned by `DomainDecomposition.refine`, see #623 +- #527 : Raise a `ValueError` in `Geometry.export` if a patch has no spline mapping, instead of an `AttributeError` after creating the file, see #631 +- #527 : Use the physical dimension of the mappings in `Geometry.from_topological_domain` and `from_discrete_mapping` +- #527 : Pair the patches of a geometry file with the domain interiors by name, not by position +- #527 : Avoid setting the `shape` of NumPy arrays (deprecated since NumPy 2.5) in `KroneckerLinearSolver` +- #527 : Fix error message in `psydac.api.discretization.get_max_degree_of_one_space` - #596 : Allow export to VTK with a prescribed grid containing exact domain bounds - #596 : Fix global geometric projectors with complex `dtype` - #576 : Fix bug in `TensorFemSpace.eval_field` caused by round-off at MPI subdomain boundaries - #576 : Require `sympde==0.19.3` which fixes a bug in the linearity checks - #579 : Require `h5py>=3.16` which installs correctly with `setuptools>=81.0` +- #576 : Fix bug in `TensorFemSpace.eval_field` caused by round-off at MPI subdomain boundaries - #579 : Don't run postprocessing unit tests with `pytest-xdist` because `h5py` is not thread-safe - #579 : Return error code on failure of the `psydac test` and `psydac compile` commands - #577 : Fix installation following release of Pyccel 2.2 - #571 : Fix correct application of the sum factorization algorithm -- #570 : Optimize PSYDAC logo -- #565 : Expand editable install info in `README.md` +- #567 : Fix parallel creation of folder `__psydac__` in `psydac.api.fem_bilinear_form` - #566 : Fix command `psydac test --mpi` on Ubuntu machines - [DEVELOPER] Add missing 'description' properties (required!) to our GitHub actions - [DEVELOPER] Update CI installation of `petsc4py` after release of `setuptools` 81.0 @@ -35,6 +60,14 @@ All notable changes to this project will be documented in this file. ### Changed +- #527 : Pass a domain decomposition `ddm` and physical dimension `pdim` to `Geometry.__init__`; remove arguments `ncells`, `periodic`, `filename`, `comm` and `mpi_dims_mask` +- #527 : Use the domain decomposition of the mapping's space in `Geometry.from_discrete_mapping`; remove arguments `comm` and `mpi_dims_mask` +- #527 : Key `Geometry.mappings` by the names of the domain interiors +- #527 : Return Cartesian (not homogeneous) control points from the NURBS functions in `psydac.cad.gallery`, see #608 +- #527 : Keep the given SymPDE domain in `discretize(domain, filename=...)`, after checking that it matches the geometry file +- #527 : Do not import the submodules of `psydac.cad` in its `__init__.py` +- #527 : Use latest version of Igakit (commit dalcinl/igakit@2a0c9ee of 2026/10/14) which supports NumPy >= 2.5 +- #527 : Raise a `ValueError` with an explanation, instead of an `AssertionError`, in `psydac.fem.partitioning.partition_coefficients` when an MPI process owns too few coefficients - #576 : Use latest version of Igakit (commit dalcinl/igakit@92ee097 of 2026/07/24) which supports NumPy >= 2.4 - #595 : Use PETSc 3.25.5 whose Python bindings `petsc4py` are built correctly with `cython>=3` - #580 : Use PETSc 3.25.0 whose Python bindings `petsc4py` install correctly with `setuptools>=81.0` @@ -42,6 +75,7 @@ All notable changes to this project will be documented in this file. - #579 : Require `numpy>=2.1` to support Python >= 3.10 - #579 : Require `pytest>=9.0` and use `pytest.toml` instead of `pytest.ini` for Pytest configuration - #579 : Move coverage configuration from `pyproject.toml` to `psydac/pytest.toml` +- #570 : Optimize PSYDAC logo - [DEVELOPER] Update GitHub Actions for repository checkout and Python setup - [DEVELOPER] Rename actions: `macos/ubuntu_install` -> `macos/ubuntu_installations` - [DEVELOPER] Do not check file changes to trigger testing workflow on PRs @@ -52,6 +86,9 @@ All notable changes to this project will be documented in this file. ### Removed +- #527 : Remove method `Geometry.read`, use constructor `Geometry.from_file` instead +- #527 : Remove property `Geometry.is_parallel` + ## [1.0.0] - 2026-01-19 The first official release on PyPI. A complete overhaul since version 0.1. diff --git a/examples/old_examples/poisson_2d_mapping.py b/examples/old_examples/poisson_2d_mapping.py index 565fcdebb..c64eaec03 100644 --- a/examples/old_examples/poisson_2d_mapping.py +++ b/examples/old_examples/poisson_2d_mapping.py @@ -632,7 +632,7 @@ def main(*, test_case, ncells, degree, nquads, map_discrete = SplineMapping.from_mapping(V, map_analytic) # Write discrete geometry to HDF5 file t0 = time() - geometry = Geometry.from_discrete_mapping(map_discrete, comm=mpi_comm) + geometry = Geometry.from_discrete_mapping(map_discrete) geometry.export('geo.h5') t1 = time() timing['export'] += t1-t0 diff --git a/examples/performance/compare_3d_matrix_assembly_speed.py b/examples/performance/compare_3d_matrix_assembly_speed.py index ebb04bf7f..ac2c68227 100644 --- a/examples/performance/compare_3d_matrix_assembly_speed.py +++ b/examples/performance/compare_3d_matrix_assembly_speed.py @@ -80,7 +80,7 @@ def make_square_torus_geometry_3d(ncells, degree, comm=None): mapping = SquareTorus('S') map_discrete = SplineMapping.from_mapping(V_h, mapping.get_callable_mapping()) - geometry = Geometry.from_discrete_mapping(map_discrete, comm=comm) + geometry = Geometry.from_discrete_mapping(map_discrete) if mpi_rank == 0: if not os.path.isdir('geometry'): diff --git a/psydac/api/ast/tests/test_poisson.py b/psydac/api/ast/tests/test_poisson.py index d6470a32f..c6025cd32 100644 --- a/psydac/api/ast/tests/test_poisson.py +++ b/psydac/api/ast/tests/test_poisson.py @@ -12,8 +12,7 @@ from sympde.calculus import grad, dot from sympde.topology import ScalarFunctionSpace from sympde.topology import elements_of, LogicalExpr -from sympde.topology import Square -from sympde.topology import Mapping, IdentityMapping, PolarMapping +from sympde.topology import Domain from sympde.expr import integral from sympde.expr import LinearForm from sympde.expr import BilinearForm @@ -38,9 +37,8 @@ #============================================================================== def test_codegen(): - domain = Square() - M = Mapping('M',2) - domain = M(domain) + domain = Domain.from_file(filename) + M = domain.mapping V = ScalarFunctionSpace('V', domain) u,v = elements_of(V, names='u,v') diff --git a/psydac/api/ast/tests/test_system_1.py b/psydac/api/ast/tests/test_system_1.py index 95fc81744..3694b5503 100644 --- a/psydac/api/ast/tests/test_system_1.py +++ b/psydac/api/ast/tests/test_system_1.py @@ -12,7 +12,7 @@ from sympde.topology import VectorFunctionSpace from sympde.topology import element_of from sympde.topology import Square -from sympde.topology import Mapping#, IdentityMapping, PolarMapping +from sympde.topology import Domain, Mapping from sympde.expr import integral from sympde.expr import LinearForm from sympde.expr import BilinearForm @@ -61,7 +61,9 @@ def test_codegen(): h1norm_F = SemiNorm(error, domain, kind='h1') # Create computational domain from topological domain - domain_h = discretize(domain, filename=filename) + # The spaces are defined on the logical domain, but the geometry needs + # the mapped domain defined in the file + domain_h = discretize(Domain.from_file(filename), filename=filename) # Discrete spaces Vh = discretize(V, domain_h) diff --git a/psydac/api/ast/tests/test_system_2.py b/psydac/api/ast/tests/test_system_2.py index 72ea2dbcd..3aecbddeb 100644 --- a/psydac/api/ast/tests/test_system_2.py +++ b/psydac/api/ast/tests/test_system_2.py @@ -12,7 +12,7 @@ from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace from sympde.topology import element_of from sympde.topology import Square -from sympde.topology import Mapping +from sympde.topology import Domain, Mapping from sympde.expr import integral from sympde.expr import LinearForm from sympde.expr import BilinearForm @@ -59,7 +59,9 @@ def test_codegen(): l = LinearForm((q,v), int_0(f1*q[0]+f2*q[1]+v)) # Create computational domain from topological domain - domain_h = discretize(domain, filename=filename) + # The spaces are defined on the logical domain, but the geometry needs + # the mapped domain defined in the file + domain_h = discretize(Domain.from_file(filename), filename=filename) # Discrete spaces Vh = discretize(V1*V2, domain_h) diff --git a/psydac/api/ast/tests/test_system_3.py b/psydac/api/ast/tests/test_system_3.py index d5bc276a4..576df45b4 100644 --- a/psydac/api/ast/tests/test_system_3.py +++ b/psydac/api/ast/tests/test_system_3.py @@ -12,7 +12,7 @@ from sympde.topology import VectorFunctionSpace from sympde.topology import element_of from sympde.topology import Square -from sympde.topology import Mapping#, IdentityMapping, PolarMapping +from sympde.topology import Domain, Mapping from sympde.expr import integral from sympde.expr import BilinearForm from sympde.expr.evaluation import TerminalExpr @@ -54,7 +54,9 @@ def test_codegen(): b = BilinearForm((u,v), int_0(div(v)*div(u))) # Create computational domain from topological domain - domain_h = discretize(domain, filename=filename) + # The spaces are defined on the logical domain, but the geometry needs + # the mapped domain defined in the file + domain_h = discretize(Domain.from_file(filename), filename=filename) # Discrete spaces Vh = discretize(V, domain_h) diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index a7a1a5e31..294f002de 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -110,7 +110,7 @@ def get_max_degree_of_one_space(Vh): Returns ------- list[int] - The maximum polynomial degre of Vh with respect to each coordinate. + The maximum polynomial degree of Vh with respect to each coordinate. """ @@ -125,7 +125,7 @@ def get_max_degree_of_one_space(Vh): return [max(p) for p in zip(*degree)] else: - raise TypeError(f'Type({V}) not understood') + raise TypeError(f'Type {type(Vh).__name__} not understood') def get_max_degree(*spaces): @@ -454,11 +454,8 @@ def discretize_space(V, domain_h, *, degree=None, multiplicity=None, knots=None, connectivity = construct_connectivity(domain) if isinstance(domain_h, Geometry) and all(domain_h.mappings.values()): - # from a discrete geoemtry - if interiors[0].name in domain_h.mappings: - mappings = [domain_h.mappings[inter.name] for inter in interiors] - else: - mappings = [domain_h.mappings[inter.logical_domain.name] for inter in interiors] + # from a discrete geometry + mappings = [domain_h.mappings[inter.name] for inter in interiors] # Get all the FEM spaces from the mapping and convert their coeff_space at the dtype needed spaces = [change_dtype(m.space, dtype) for m in mappings] @@ -578,7 +575,7 @@ def discretize_domain(domain, *, filename=None, ncells=None, periodic=None, comm raise ValueError("Cannot provide both 'filename' and 'ncells'") elif filename: - return Geometry(filename=filename, comm=comm, mpi_dims_mask=mpi_dims_mask) + return Geometry.from_file(filename, domain=domain, comm=comm, mpi_dims_mask=mpi_dims_mask) elif ncells: return Geometry.from_topological_domain(domain, ncells, periodic=periodic, comm=comm, mpi_dims_mask=mpi_dims_mask) diff --git a/psydac/api/postprocessing.py b/psydac/api/postprocessing.py index 55d088b0e..8a97b6ad1 100644 --- a/psydac/api/postprocessing.py +++ b/psydac/api/postprocessing.py @@ -1042,15 +1042,16 @@ def _process_domain(self): spl_maps = domain_h.mappings if domain_h.mappings is not None else {} if isinstance(domain.interior, InteriorDomain): - self._mappings[domain.name] = spl_maps.get(domain.logical_domain.name, domain.mapping) + self._mappings[domain.name] = spl_maps.get(domain.interior.name, + domain.mapping) else: if isinstance(domain.mapping, MultiPatchMapping): for interior in domain.interior.as_tuple(): - self._mappings[interior.name] = spl_maps.get(interior.logical_domain.name, \ + self._mappings[interior.name] = spl_maps.get(interior.name, \ domain.mapping.mappings[interior.logical_domain]) else: for interior in domain.interior.as_tuple(): - self._mappings[interior.name] = spl_maps.get(interior.logical_domain.name, \ + self._mappings[interior.name] = spl_maps.get(interior.name, \ interior.mapping) else: domain = self._domain diff --git a/psydac/api/tests/test_2d_complex.py b/psydac/api/tests/test_2d_complex.py index 85c2479dc..1a637776f 100644 --- a/psydac/api/tests/test_2d_complex.py +++ b/psydac/api/tests/test_2d_complex.py @@ -9,8 +9,7 @@ from mpi4py import MPI import pytest -import numpy as np -from sympy import pi, cos, sin, symbols, conjugate, exp +from sympy import pi, cos, sin, symbols, exp from sympy import Tuple, Matrix from sympy import lambdify @@ -22,7 +21,7 @@ from sympde.topology import NormalVector from sympde.topology import Union from sympde.topology import Domain, Square -from sympde.topology import IdentityMapping, AffineMapping, PolarMapping +from sympde.topology import IdentityMapping from sympde.expr import BilinearForm, LinearForm, integral from sympde.expr import Norm, SemiNorm from sympde.expr import find, EssentialBC @@ -363,7 +362,7 @@ def test_complex_poisson_2d_multipatch(): A = Square('A',bounds1=(0, 0.5), bounds2=(0, 1)) B = Square('B',bounds1=(0.5, 1.), bounds2=(0, 1)) - domain = Domain.join([A, B], [((0, 0, 1), (1, 0, -1))], 'domain') + domain = Domain.join([A, B], [((0, 0, 1), (1, 0, -1), 1)], 'domain') x, y = domain.coordinates @@ -518,8 +517,8 @@ def teardown_function(): else: from psydac.fem.plotting_utilities import get_plotting_grid, get_grid_vals - from psydac.fem.plotting_utilities import get_patch_knots_gridlines, my_small_plot - from psydac.api.tests.build_domain import build_11_patch_pretzel + from psydac.fem.plotting_utilities import my_small_plot + from psydac.api.tests.build_domain import build_11_patch_pretzel from psydac.feec.pull_push import pull_2d_hcurl domain = build_11_patch_pretzel() diff --git a/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py b/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py index 336066c71..2ae7ea7b5 100644 --- a/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py +++ b/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py @@ -7,25 +7,22 @@ from pathlib import Path import pytest -import numpy as np from mpi4py import MPI from sympy import pi, sin, cos, Tuple, Matrix -from sympde.calculus import grad, dot, curl, cross +from sympde.calculus import dot, curl, cross from sympde.calculus import minus, plus from sympde.topology import VectorFunctionSpace from sympde.topology import elements_of from sympde.topology import NormalVector -from sympde.topology import Square, Domain -from sympde.topology import IdentityMapping, PolarMapping +from sympde.topology import Domain from sympde.expr.expr import LinearForm, BilinearForm from sympde.expr.expr import integral from sympde.expr.expr import Norm -from sympde.expr.equation import find, EssentialBC +from sympde.expr.equation import find from psydac.api.discretization import discretize from psydac.api.tests.build_domain import build_11_patch_pretzel, build_2_patch_annulus -from psydac.fem.basic import FemField from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL from psydac.feec.pull_push import pull_2d_hcurl @@ -204,10 +201,10 @@ def teardown_function(): if __name__ == '__main__': - from collections import OrderedDict - from sympy import lambdify + from collections import OrderedDict + from sympy import lambdify from psydac.fem.plotting_utilities import get_plotting_grid, get_grid_vals - from psydac.fem.plotting_utilities import get_patch_knots_gridlines, my_small_plot + from psydac.fem.plotting_utilities import my_small_plot domain = build_11_patch_pretzel() x,y = domain.coordinates diff --git a/psydac/api/tests/test_2d_multipatch_mapping_poisson.py b/psydac/api/tests/test_2d_multipatch_mapping_poisson.py index f730f6bc6..3f7a0de19 100644 --- a/psydac/api/tests/test_2d_multipatch_mapping_poisson.py +++ b/psydac/api/tests/test_2d_multipatch_mapping_poisson.py @@ -15,7 +15,7 @@ from sympde.calculus import minus, plus from sympde.topology import ScalarFunctionSpace from sympde.topology import elements_of -from sympde.topology import NormalVector, Union +from sympde.topology import NormalVector from sympde.topology import Square, Domain from sympde.topology import IdentityMapping, PolarMapping, AffineMapping from sympde.expr.expr import LinearForm, BilinearForm @@ -131,7 +131,7 @@ def test_poisson_2d_2_patches_dirichlet_1(): assert ( abs(h1_error - expected_h1_error) < 1e-7 ) #------------------------------------------------------------------------------ -def test_poisson_2d_3_patches_dirichlet_2(): +def test_poisson_2d_3_patches_dirichlet_2(tmp_path): mapping_1 = IdentityMapping('M1', 2) mapping_2 = PolarMapping ('M2', 2, c1 = 0., c2 = 0.5, rmin = 0., rmax=1.) @@ -141,12 +141,12 @@ def test_poisson_2d_3_patches_dirichlet_2(): B = Square('B',bounds1=(0.5, 1.), bounds2=(0, np.pi)) C = Square('C',bounds1=(0.5, 1.), bounds2=(np.pi-0.5, np.pi + 1)) - D1 = mapping_1(A) - D2 = mapping_2(B) - D3 = mapping_3(C) + D1 = mapping_1(A) + D2 = mapping_2(B) + D3 = mapping_3(C) - connectivity = [((0,1,1),(1,1,-1)), ((1,1,1),(2,1,-1))] patches = [D1, D2, D3] + connectivity = [((0, 1, 1), (1, 1,-1), 1), ((1, 1, 1), (2, 1,-1), 1)] domain = Domain.join(patches, connectivity, 'domain') x,y = domain.coordinates @@ -155,7 +155,7 @@ def test_poisson_2d_3_patches_dirichlet_2(): l2_error, h1_error, uh = run_poisson_2d(solution, f, domain, ncells=[2**2,2**2], degree=[2,2]) - plot_fn=f'uh_multipatch_poisson.pdf' + plot_fn = str(tmp_path / 'uh_multipatch_poisson.pdf') plot_field(fem_field=uh, Vh=uh.space, domain=domain, title='uh', filename=plot_fn, hide_plot=True) expected_l2_error = 0.0019402242901236006 @@ -252,8 +252,11 @@ def test_poisson_2d_4_patch_dirichlet_0(): D3 = mapping_3(C) D4 = mapping_4(D) - connectivity = [((0,1,1),(1,1,-1)), ((2,1,1),(3,1,-1)), ((0,0,1),(2,0,-1)),((1,0,1),(3,0,-1))] patches = [D1, D2, D3, D4] + connectivity = [((0, 1, 1), (1, 1,-1), 1), + ((2, 1, 1), (3, 1,-1), 1), + ((0, 0, 1), (2, 0,-1), 1), + ((1, 0, 1), (3, 0,-1), 1)] domain = Domain.join(patches, connectivity, 'domain') x,y = domain.coordinates diff --git a/psydac/api/tests/test_2d_multipatch_poisson.py b/psydac/api/tests/test_2d_multipatch_poisson.py index 297975eba..b4f230ebb 100644 --- a/psydac/api/tests/test_2d_multipatch_poisson.py +++ b/psydac/api/tests/test_2d_multipatch_poisson.py @@ -78,8 +78,8 @@ def test_poisson_2d_2_patch_dirichlet_0(): A = Square('A', bounds1=(0, 0.5), bounds2=(0, 1)) B = Square('B', bounds1=(0.5, 1), bounds2=(0, 1)) - connectivity = [((0,0,1), (1,0,-1), 1)] patches = [A, B] + connectivity = [((0, 0, 1), (1, 0,-1), 1)] domain = Domain.join(patches, connectivity, 'domain') x,y = domain.coordinates @@ -100,8 +100,8 @@ def test_poisson_2d_2_patch_dirichlet_1(): A = Square('A', bounds1=(0, 0.5), bounds2=(0, 1)) B = Square('B', bounds1=(0.5, 1), bounds2=(0, 1)) - connectivity = [((0,0,1), (1,0,-1), 1)] patches = [A, B] + connectivity = [((0, 0, 1), (1, 0,-1), 1)] domain = Domain.join(patches, connectivity, 'domain') x,y = domain.coordinates @@ -121,8 +121,8 @@ def test_poisson_2d_2_patch_dirichlet_2(): A = Square('A', bounds1=(0, 0.5), bounds2=(0, 1)) B = Square('B', bounds1=(0.5, 1), bounds2=(0, 1)) - connectivity = [((0,0,1), (1,0,-1), 1)] patches = [A, B] + connectivity = [((0, 0, 1), (1, 0,-1), 1)] domain = Domain.join(patches, connectivity, 'domain') x,y = domain.coordinates @@ -150,8 +150,8 @@ def test_poisson_2d_2_patch_dirichlet_3(): D1 = M1(A) D2 = M2(B) - connectivity = [((0,0,1), (1,0,1), 1)] patches = [D1, D2] + connectivity = [((0, 0, 1), (1, 0, 1), 1)] domain = Domain.join(patches, connectivity, 'domain') x,y = domain.coordinates @@ -182,8 +182,8 @@ def test_poisson_2d_2_patch_dirichlet_4(): D1 = M1(A) D2 = M2(B) - connectivity = [((0,0,-1), (1,0,-1), 1)] - patches = [D1,D2] + patches = [D1, D2] + connectivity = [((0, 0, -1), (1, 0, -1), 1)] domain = Domain.join(patches, connectivity, 'domain') x,y = domain.coordinates diff --git a/psydac/api/tests/test_2d_navier_stokes.py b/psydac/api/tests/test_2d_navier_stokes.py index 8118f277a..bbd98ff77 100644 --- a/psydac/api/tests/test_2d_navier_stokes.py +++ b/psydac/api/tests/test_2d_navier_stokes.py @@ -8,17 +8,11 @@ import pytest import numpy as np -from sympy import pi, cos, sin, sqrt, exp, ImmutableDenseMatrix as Matrix, Tuple, lambdify +from sympy import pi, cos, sin, exp, ImmutableDenseMatrix as Matrix, Tuple from scipy.sparse.linalg import spsolve -from scipy.sparse.linalg import gmres as sp_gmres -from scipy.sparse.linalg import minres as sp_minres -from scipy.sparse.linalg import cg as sp_cg -from scipy.sparse.linalg import bicg as sp_bicg -from scipy.sparse.linalg import bicgstab as sp_bicgstab - -from sympde.calculus import grad, dot, inner, div, curl, cross -from sympde.calculus import Transpose, laplace -from sympde.topology import NormalVector + +from sympde.calculus import grad, dot, inner, div +from sympde.calculus import Transpose from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace from sympde.topology import ProductSpace from sympde.topology import element_of, elements_of @@ -27,7 +21,6 @@ from sympde.expr import Norm from sympde.expr import find, EssentialBC from sympde.core import Constant -from sympde.expr import TerminalExpr from sympde.expr import linearize from psydac.api.essential_bc import apply_essential_bc @@ -39,7 +32,7 @@ from psydac.linalg.stencil import * from psydac.linalg.block import * from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL -from psydac.utilities.utils import refine_array_1d, animate_field, split_space, split_field +from psydac.utilities.utils import animate_field, split_space, split_field from psydac.linalg.solvers import inverse from mpi4py import MPI @@ -467,9 +460,7 @@ def teardown_function(): else: - import matplotlib.pyplot as plt from matplotlib import animation - from time import time Tf = 3. dt_h = 0.05 @@ -479,7 +470,7 @@ def teardown_function(): solutions, p_h, domain, domain_h = run_time_dependent_navier_stokes_2d(filename, dt_h=dt_h, nt=nt, scipy=False) domain = domain.logical_domain - mapping = domain_h.mappings['patch_0'] + mapping, = domain_h.mappings.values() anim = animate_field(solutions, domain, mapping, res=(150,150), progress=True) anim.save('animated_fields_{}_{}.mp4'.format(str(Tf).replace('.','_'), str(dt_h).replace('.','_')), writer=animation.FFMpegWriter(fps=60)) diff --git a/psydac/api/tests/test_discretization.py b/psydac/api/tests/test_discretization.py new file mode 100644 index 000000000..27920c275 --- /dev/null +++ b/psydac/api/tests/test_discretization.py @@ -0,0 +1,60 @@ +#---------------------------------------------------------------------------# +# This file is part of PSYDAC which is released under MIT License. See the # +# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE # +# for full license details. # +#---------------------------------------------------------------------------# +import numpy as np +import pytest + +from psydac.api.discretization import get_max_degree, get_max_degree_of_one_space +from psydac.ddm.cart import DomainDecomposition +from psydac.fem.splines import SplineSpace +from psydac.fem.tensor import TensorFemSpace +from psydac.fem.vector import VectorFemSpace + +#============================================================================== +def make_tensor_space_2d(degree: tuple[int, int], + domain_decomposition: DomainDecomposition) -> TensorFemSpace: + """ + Create a 2D tensor-product spline space on a uniform grid. + + Parameters + ---------- + degree : tuple[int, int] + The spline degree along each direction. + + domain_decomposition : DomainDecomposition + The decomposition of the 2D grid, which also gives its number of cells. + + Returns + ------- + TensorFemSpace + The spline space. + """ + spaces = [SplineSpace(degree=p, grid=np.linspace(0.0, 1.0, n + 1)) + for p, n in zip(degree, domain_decomposition.ncells)] + return TensorFemSpace(domain_decomposition, *spaces) + +#============================================================================== +def test_get_max_degree_of_one_space() -> None: + + domain_decomposition = DomainDecomposition([3, 5], [False, False]) + V1 = make_tensor_space_2d((2, 3), domain_decomposition) + V2 = make_tensor_space_2d((3, 1), domain_decomposition) + + assert list(get_max_degree_of_one_space(V1)) == [2, 3] + assert list(get_max_degree_of_one_space(VectorFemSpace(V1, V2))) == [3, 3] + + with pytest.raises(TypeError, match='SplineSpace'): + get_max_degree_of_one_space(V1.spaces[0]) + +#============================================================================== +def test_get_max_degree() -> None: + + domain_decomposition = DomainDecomposition([3, 5], [False, False]) + V1 = make_tensor_space_2d((2, 3), domain_decomposition) + V2 = make_tensor_space_2d((3, 1), domain_decomposition) + V3 = make_tensor_space_2d((1, 4), domain_decomposition) + + assert list(get_max_degree(V1, V2)) == [3, 3] + assert list(get_max_degree(VectorFemSpace(V1, V2), V3)) == [3, 4] diff --git a/psydac/api/tests/test_postprocessing.py b/psydac/api/tests/test_postprocessing.py index 73cc941f9..815b7140d 100644 --- a/psydac/api/tests/test_postprocessing.py +++ b/psydac/api/tests/test_postprocessing.py @@ -17,6 +17,7 @@ from psydac.api.discretization import discretize from psydac.fem.basic import FemField from psydac.fem.tensor import TensorFemSpace +from psydac.mapping.discrete import SplineMapping from psydac.utilities.utils import refine_array_1d from psydac.feec.pull_push import (push_2d_hcurl, push_2d_h1, @@ -57,7 +58,7 @@ def build_2_mapped_squares(): D2 = mapping_2(B) patches = [D1, D2] - connectivity = [((0,1,1),(1,1,-1))] + connectivity = [((0, 1, 1), (1, 1,-1), 1)] return Domain.join(patches, connectivity, 'domain') @@ -66,7 +67,7 @@ def build_2_squares(): B = Square('B',bounds1=(0.5, 1.), bounds2=(np.pi/2, np.pi)) patches = [A, B] - connectivity = [((0,1,1),(1,1,-1))] + connectivity = [((0, 1, 1), (1, 1,-1), 1)] return Domain.join(patches, connectivity, 'domain') @@ -75,10 +76,9 @@ def build_2_cubes(): B = Cube('B',bounds1=(0.5, 1.), bounds2=(np.pi/2, np.pi), bounds3=(0, 1)) patches = [A, B] - connectivity = [((0,1,1),(1,1,-1))] + connectivity = [((0, 1, 1), (1, 1,-1), (1, 1, 1))] return Domain.join(patches, connectivity, 'domain') - ############################################################################### # Output Manager tests # ############################################################################### @@ -495,6 +495,9 @@ def test_reconstruct_DerhamSequence_discrete_domain(geometry, seq, dtype): fields_file='test_reconstruct_DerhamSequence_discrete_domain.h5' ) + # The post-processor must use the spline mappings of the geometry file + assert all(isinstance(F, SplineMapping) for F in Pm._mappings.values()) + Om2 = OutputManager( 'test_reconstruct_DerhamSequence_discrete_domain_2.yml', 'test_reconstruct_DerhamSequence_discrete_domain.h5' @@ -538,8 +541,8 @@ def test_reconstruct_multipatch(dtype): A = Square('A',bounds1=bounds1, bounds2=bounds2_A) B = Square('B',bounds1=bounds1, bounds2=bounds2_B) - connectivity = [((0,1,1),(1,1,-1))] - patches = [A,B] + patches = [A, B] + connectivity = [((0, 1, 1), (1, 1,-1), 1)] domain = Domain.join(patches, connectivity, 'domain') Va = ScalarFunctionSpace('Va', A) diff --git a/psydac/cad/__init__.py b/psydac/cad/__init__.py index 6f39d9a03..419109b64 100644 --- a/psydac/cad/__init__.py +++ b/psydac/cad/__init__.py @@ -3,10 +3,3 @@ # LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE # # for full license details. # #---------------------------------------------------------------------------# - -__all__ = ['geometry'] - -from psydac.cad import geometry -from psydac.cad import cad -from psydac.cad import gallery -from psydac.cad import utils diff --git a/psydac/cad/cad.py b/psydac/cad/cad.py index abc39baab..1c3556682 100644 --- a/psydac/cad/cad.py +++ b/psydac/cad/cad.py @@ -3,35 +3,67 @@ # LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE # # for full license details. # #---------------------------------------------------------------------------# +from typing import Callable, Iterable + import numpy as np +from mpi4py import MPI from psydac.fem.splines import SplineSpace from psydac.fem.tensor import TensorFemSpace from psydac.fem.basic import FemField +from psydac.linalg.stencil import StencilVector from psydac.mapping.discrete import SplineMapping, NurbsMapping from psydac.ddm.cart import DomainDecomposition #============================================================================== -def translate(mapping, displ): +def translate(mapping : SplineMapping, displ : Iterable[float]): + """ + Translate a CAD geometry by a given vector displacement. + + Translate a single-patch CAD geometry (given as a spline or NURBS mapping) + by adding a given displacement vector to the control points of the mapping. + The weights of a NURBS mapping are not changed. + + Parameters + ---------- + mapping : SplineMapping + The discrete mapping to be translated, which represents a CAD geometry. + + displ : Iterable[float] + The vector displacement by which to translate the geometry. + + Returns + ------- + SplineMapping + A new discrete mapping representing the translated geometry. It is a + NurbsMapping if the input mapping is a NurbsMapping. + """ + assert isinstance(mapping, SplineMapping) + assert isinstance(displ, Iterable) + displ = np.array(displ) - assert( mapping.pdim == len(displ) ) + assert mapping.pdim == len(displ) pdim = mapping.pdim space = mapping.space control_points = mapping.control_points - fields = [FemField( space ) for d in range( pdim )] + fields = [FemField(space) for d in range(pdim)] # Get spline coefficients for each coordinate X_i starts = space.coeff_space.starts ends = space.coeff_space.ends - idx_to = tuple( slice( s, e+1 ) for s,e in zip( starts, ends ) ) - for i,field in enumerate( fields ): + idx_to = tuple(slice(s, e+1) for s, e in zip(starts, ends)) + for i, field in enumerate(fields): idx_from = tuple(list(idx_to)+[i]) field.coeffs[idx_to] = control_points[idx_from] + displ[i] field.coeffs.update_ghost_regions() - return SplineMapping( *fields ) + # The control points are Cartesian, hence the weights are not affected + if isinstance(mapping, NurbsMapping): + return NurbsMapping(*fields, mapping.weights_field.copy()) + + return SplineMapping(*fields) #============================================================================== def elevate(mapping, axis, times): @@ -41,68 +73,15 @@ def elevate(mapping, axis, times): Note: we are using igakit for the moment, until we implement the elevation degree algorithm in psydac """ - try: - from igakit.nurbs import NURBS - except: - raise ImportError('Could not find igakit.') - - assert( isinstance(mapping, (SplineSpace, NurbsMapping)) ) - assert( isinstance(times, int) ) - assert( isinstance(axis, int) ) - - space = mapping.space - domain_decomposition = space.domain_decomposition - pdim = mapping.pdim - - knots = [V.knots for V in space.spaces] - degree = [V.degree for V in space.spaces] - shape = [V.nbasis for V in space.spaces] - points = np.zeros(shape+[mapping.pdim]) - for i,f in enumerate( mapping._fields ): - points[...,i] = f._coeffs.toarray().reshape(shape) - - weights = None - if isinstance(mapping, NurbsMapping): - weights = mapping._weights_field._coeffs.toarray().reshape(shape) - - for i in range(pdim): - points[...,i] /= weights[...] - - # degree elevation using igakit - nrb = NURBS(knots, points, weights=weights) - nrb = nrb.clone().elevate(axis, times) - - spaces = [SplineSpace(degree=p, knots=u) for p,u in zip( nrb.degree, nrb.knots )] - space = TensorFemSpace( domain_decomposition, *spaces ) - fields = [FemField( space ) for d in range( pdim )] - - # Get spline coefficients for each coordinate X_i - starts = space.coeff_space.starts - ends = space.coeff_space.ends - idx_to = tuple( slice( s, e+1 ) for s,e in zip( starts, ends ) ) - for i,field in enumerate( fields ): - idx_from = tuple(list(idx_to)+[i]) - idw_from = tuple(idx_to) - if isinstance(mapping, NurbsMapping): - field.coeffs[idx_to] = nrb.points[idx_from] * nrb.weights[idw_from] - - else: - field.coeffs[idx_to] = nrb.points[idx_from] - - field.coeffs.update_ghost_regions() - - if isinstance(mapping, NurbsMapping): - weights_field = FemField( space ) + assert isinstance(mapping, SplineMapping) + assert isinstance(times, int) + assert isinstance(axis, int) - idx_from = idx_to - weights_field.coeffs[idx_to] = nrb.weights[idx_from] - weights_field.coeffs.update_ghost_regions() + # The decomposition does not change, hence each process only needs the + # coefficients owned by its neighbours that are stored in its ghost regions + nrb = _to_igakit(mapping, _local_array).elevate(axis, times) - fields.append( weights_field ) - - return NurbsMapping( *fields ) - - return SplineMapping( *fields ) + return _from_igakit(nrb, mapping, mapping.space.domain_decomposition) #============================================================================== @@ -114,81 +93,119 @@ def refine(mapping, axis, values): Note: we are using igakit for the moment, until we implement the knot insertion algorithm in psydac """ - try: - from igakit.nurbs import NURBS - except: - raise ImportError('Could not find igakit.') - - assert( isinstance(mapping, (SplineSpace, NurbsMapping)) ) - assert( isinstance(values, (list, tuple)) ) - assert( isinstance(axis, int) ) - - space = mapping.space - domain_decomposition = space.domain_decomposition - pdim = mapping.pdim - - knots = [V.knots for V in space.spaces] - degree = [V.degree for V in space.spaces] - shape = [V.nbasis for V in space.spaces] - points = np.zeros(shape+[mapping.pdim]) - for i,f in enumerate( mapping._fields ): - points[...,i] = f._coeffs.toarray().reshape(shape) - - weights = None - if isinstance(mapping, NurbsMapping): - weights = mapping._weights_field._coeffs.toarray().reshape(shape) - - for i in range(pdim): - points[...,i] /= weights[...] - - # degree elevation using igakit - nrb = NURBS(knots, points, weights=weights) - nrb = nrb.clone().refine(axis, values) + assert isinstance(mapping, SplineMapping) + assert isinstance(values, (list, tuple)) + assert isinstance(axis, int) - spaces = [SplineSpace(degree=p, knots=u) for p,u in zip( nrb.degree, nrb.knots )] + # The new balanced decomposition may move the subdomain boundaries far + # away, hence each process needs all the coefficients + nrb = _to_igakit(mapping, _global_array).refine(axis, values) - ncells = list(domain_decomposition.ncells) - ncells[axis] += len(values) - domain_decomposition = DomainDecomposition(ncells, domain_decomposition.periods, comm=domain_decomposition.comm) + return _from_igakit(nrb, mapping) - space = TensorFemSpace( domain_decomposition, *spaces ) - fields = [FemField( space ) for d in range( pdim )] +#============================================================================== +def _to_igakit(mapping: SplineMapping, to_array: Callable[[StencilVector], np.ndarray]): + """ + Create an igakit NURBS object from a spline or NURBS mapping. + + Parameters + ---------- + mapping : SplineMapping + The input mapping (possibly a NurbsMapping). + + to_array : Callable[[StencilVector], np.ndarray] + The function which returns the coefficients of a field as an array + with the global shape, which must be correct where it is used. + + Returns + ------- + igakit.nurbs.NURBS + The igakit object with the Cartesian control points and the weights + of the mapping. + """ + try: + from igakit.nurbs import NURBS + except ImportError as err: + raise ImportError('Could not find igakit.') from err - # Get spline coefficients for each coordinate X_i - starts = space.coeff_space.starts - ends = space.coeff_space.ends - idx_to = tuple( slice( s, e+1 ) for s,e in zip( starts, ends ) ) - for i,field in enumerate( fields ): - idx_from = tuple(list(idx_to)+[i]) - idw_from = tuple(idx_to) - if isinstance(mapping, NurbsMapping): - field.coeffs[idx_to] = nrb.points[idx_from] * nrb.weights[idw_from] + knots = [V.knots for V in mapping.space.spaces] + shape = [V.nbasis for V in mapping.space.spaces] - else: - field.coeffs[idx_to] = nrb.points[idx_from] + # igakit requires at least 2 coordinates, the extra one is zero + points = np.zeros(shape + [max(mapping.pdim, 2)]) + for i, f in enumerate(mapping.fields): + points[..., i] = to_array(f.coeffs).reshape(shape) + weights = None if isinstance(mapping, NurbsMapping): - weights_field = FemField( space ) - - idx_from = idx_to - weights_field.coeffs[idx_to] = nrb.weights[idx_from] - weights_field.coeffs.update_ghost_regions() - - fields.append( weights_field ) + weights = to_array(mapping.weights_field.coeffs).reshape(shape) - return NurbsMapping( *fields ) + # igakit expects Cartesian control points + return NURBS(knots, points, weights=weights) - return SplineMapping( *fields ) - - - -###################################### -if __name__ == '__main__': - from psydac.cad.geometry import Geometry - - geo = Geometry('square_0.h5') - mapping = geo.patches[0] - new = translate(mapping, [1., 0., 0.]) +def _from_igakit(nrb, + mapping: SplineMapping, + domain_decomposition: DomainDecomposition | None = None) -> SplineMapping: + """ + Create a mapping of the same type and dimension as a given mapping, from + an igakit NURBS object. + + Parameters + ---------- + nrb : igakit.nurbs.NURBS + The igakit object, whose control points and weights are known on + every process. + + mapping : SplineMapping + The mapping whose type (SplineMapping or NurbsMapping), physical + dimension, periodicity, and MPI communicator are used. + + domain_decomposition : DomainDecomposition, optional + The decomposition of the new mapping's domain. If not given, a new + balanced decomposition is created for the new number of cells. + + Returns + ------- + SplineMapping + The new mapping. It is a NurbsMapping if `mapping` is a NurbsMapping. + """ + spaces = [SplineSpace(degree=p, knots=u) for p, u in zip(nrb.degree, nrb.knots)] + + if domain_decomposition is not None: + ddm = domain_decomposition + else: + # The number of cells is given by the spaces, because inserting a value + # which is already a knot only increases its multiplicity (see #619) + old_ddm = mapping.space.domain_decomposition + ddm = DomainDecomposition([W.ncells for W in spaces], old_ddm.periods, comm=old_ddm.comm, + mpi_dims_mask=old_ddm.mpi_dims_mask) + + space = TensorFemSpace(ddm, *spaces) + + arrays = [nrb.points[..., i] for i in range(mapping.pdim)] + is_nurbs = isinstance(mapping, NurbsMapping) + if is_nurbs: + arrays.append(nrb.weights) + + # Each process copies the coefficients that it owns + idx = tuple(slice(s, e + 1) for s, e in zip(space.coeff_space.starts, space.coeff_space.ends)) + fields = [FemField(space) for _ in arrays] + for field, array in zip(fields, arrays): + field.coeffs[idx] = array[idx] + field.coeffs.update_ghost_regions() - geo = Geometry(patches=[mapping, new]) - geo.export('square_mp.h5') + return NurbsMapping(*fields) if is_nurbs else SplineMapping(*fields) + +def _local_array(coeffs: StencilVector) -> np.ndarray: + """Return the array of coefficients owned by the process and its neighbours.""" + coeffs.update_ghost_regions() + # In parallel, the array is zero outside of the block and the ghost regions + return coeffs.toarray(with_pads=True) + +def _global_array(coeffs: StencilVector) -> np.ndarray: + """Return the global array of coefficients on every process.""" + array = coeffs.toarray() + # In parallel, toarray() returns zeros outside of the block owned by the process + if coeffs.space.parallel: + coeffs.space.cart.comm.Allreduce(MPI.IN_PLACE, array, op=MPI.SUM) + return array diff --git a/psydac/cad/gallery.py b/psydac/cad/gallery.py index f9f2d7732..6cd71a344 100644 --- a/psydac/cad/gallery.py +++ b/psydac/cad/gallery.py @@ -3,13 +3,11 @@ # LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE # # for full license details. # #---------------------------------------------------------------------------# -import string -import random import numpy as np from psydac.fem.splines import SplineSpace from psydac.fem.tensor import TensorFemSpace -from psydac.mapping.discrete import SplineMapping, NurbsMapping +from psydac.mapping.discrete import NurbsMapping from psydac.cad.utils import plot_mapping #============================================================================== @@ -18,6 +16,8 @@ def quart_circle( rmin=0.5, rmax=1.0, center=None ): degrees = (2, 1) pdim = 2 + s = 1./np.sqrt(2) + knots = [[0.0 , 0.0 , 0.0 , 1.0 , 1.0 , 1.0], [0.0 , 0.0 , 1.0 , 1.0] ] @@ -38,16 +38,13 @@ def quart_circle( rmin=0.5, rmax=1.0, center=None ): weights = np.zeros((3,2)) j = 0 weights[0,j] = 1.0 - weights[1,j] = 0.707106781187 + weights[1,j] = s weights[2,j] = 1.0 j = 1 weights[0,j] = 1.0 - weights[1,j] = 0.707106781187 + weights[1,j] = s weights[2,j] = 1.0 - for i in range(pdim): - points[...,i] *= weights[...] - return degrees, knots, points, weights @@ -57,6 +54,8 @@ def annulus(rmin=0.5, rmax=1.0, center=None): degrees = (1, 2) pdim = 2 + s = 1./np.sqrt(2) + knots = [[0.0 , 0.0 , 1.0 , 1.0], [0.0 , 0.0 , 0.0 , 0.25 , 0.25 , 0.5 , 0.5 , 0.75 , 0.75 , 1.0 , 1.0 , 1.0] ] @@ -89,28 +88,25 @@ def annulus(rmin=0.5, rmax=1.0, center=None): weights = np.zeros((2,9)) j = 0 weights[j,0] = 1.0 - weights[j,1] = 0.707106781187 + weights[j,1] = s weights[j,2] = 1.0 - weights[j,3] = 0.707106781187 + weights[j,3] = s weights[j,4] = 1.0 - weights[j,5] = 0.707106781187 + weights[j,5] = s weights[j,6] = 1.0 - weights[j,7] = 0.707106781187 + weights[j,7] = s weights[j,8] = 1.0 j = 1 weights[j,0] = 1.0 - weights[j,1] = 0.707106781187 + weights[j,1] = s weights[j,2] = 1.0 - weights[j,3] = 0.707106781187 + weights[j,3] = s weights[j,4] = 1.0 - weights[j,5] = 0.707106781187 + weights[j,5] = s weights[j,6] = 1.0 - weights[j,7] = 0.707106781187 + weights[j,7] = s weights[j,8] = 1.0 - for i in range(pdim): - points[...,i] *= weights[...] - return degrees, knots, points, weights @@ -150,9 +146,6 @@ def circle(radius=1.0, center=None): weights[1,2] = s weights[2,2] = 1. - for i in range(pdim): - points[...,i] *= weights[...] - return degrees, knots, points, weights #============================================================================== diff --git a/psydac/cad/geometry.py b/psydac/cad/geometry.py index cf89c85d1..ba905cac7 100644 --- a/psydac/cad/geometry.py +++ b/psydac/cad/geometry.py @@ -8,20 +8,19 @@ # the topology i.e. connectivity, boundaries # For the moment, it is used as a container, that can be loaded from a file # (hdf5) -from itertools import product -from collections import abc -import string -import random -import yaml import os -import string -import random -import warnings +from typing import Iterable +from itertools import chain import numpy as np import h5py +import yaml from mpi4py import MPI +from sympde.topology import Domain, Line, Square, Cube, NCubeInterior, Mapping, NCube +from sympde.topology.basic import Union +from sympde.topology.callable_mapping import BasicCallableMapping + from psydac.fem.splines import SplineSpace from psydac.fem.tensor import TensorFemSpace from psydac.fem.partitioning import create_cart, construct_connectivity, construct_interface_spaces @@ -29,174 +28,359 @@ from psydac.linalg.block import BlockVectorSpace, BlockVector from psydac.ddm.cart import DomainDecomposition, MultiPatchDomainDecomposition +__all__ = ( + 'Geometry', + 'export_nurbs_to_hdf5', + 'import_geopdes_to_nurbs', + 'refine_knots', + 'refine_nurbs', +) -from sympde.topology import Domain, Interface, Line, Square, Cube, NCubeInterior, Mapping, NCube -from sympde.topology.basic import Union +NoneType = type(None) #============================================================================== class Geometry: """ Distributed discrete geometry that works for single and multiple patches. - The Geometry object can be created in two ways: - - case 1 : through a geometry file whos name can be given to the constructor - - case 2 : provide the ncells, the periodicity and the mapping objects of each patch. - - Parameters - ---------- - domain : Sympde.topology.Domain - The symbolic domain to be discretized. - ncells : list | tuple | dict - The number of cells of the discretized topological domain in each direction. + A Geometry object is usually created by calling + `~psydac.api.discretization.discretize` on a SymPDE domain: + `discretize(domain, filename=...)` calls `Geometry.from_file`, and + `discretize(domain, ncells=...)` calls `Geometry.from_topological_domain`. - periodic : list | tuple | dict - The periodicity of the topological domain in each direction. + It can also be created directly with one of the class methods: - mappings : dict - The Mapping of each patch. + - `Geometry.from_file`: from a geometry file (single or multi-patch). + - `Geometry.from_discrete_mapping`: from a + `~psydac.mapping.discrete.SplineMapping` (single patch). + - `Geometry.from_topological_domain`: from a topological domain, the number + of cells, and the periodicity (single or multi-patch). - filename: str - The path to the geometry file. - - comm: MPI.Comm - MPI intra-communicator. - - mpi_dims_mask: list of bool - True if the dimension is to be used in the domain decomposition (=default for each dimension). - If mpi_dims_mask[i]=False, the i-th dimension will not be decomposed. - + The constructor `__init__` is meant for internal use by these class methods. """ _ldim = None _pdim = None _patches = [] _topology = None + def __init__(self, + domain : Domain, + *, + ddm : DomainDecomposition | MultiPatchDomainDecomposition, + pdim : int, + mappings : dict[str, BasicCallableMapping | None] = None): + """ + Create a Geometry instance from an existing domain decomposition. + + This constructor is meant for internal use: the class methods should be + preferred. The number of cells and the periodicity of each patch are + obtained from the domain decomposition. + + Parameters + ---------- + domain : Sympde.topology.Domain + The symbolic topological domain to be discretized. + + ddm : DomainDecomposition | MultiPatchDomainDecomposition + The decomposition of the domain across the MPI processes: a + `DomainDecomposition` for a single patch, otherwise a + `MultiPatchDomainDecomposition` with one decomposition per patch. + + pdim : int + Number of physical dimensions of the Geometry object (pdim >= ldim). + + mappings : dict[str, BasicCallableMapping | None], optional + The callable mapping of each patch, keyed by the names of the domain + interiors, or None if the patch has no callable mapping (default: + no patch has one). Each mapping must have physical dimension `pdim`. + If a mapping is a `SplineMapping` (or a `NurbsMapping`), its space + must be defined on the domain decomposition of its patch, i.e. on + the same object as `ddm` (single patch) or as the corresponding + element of `ddm.domains` (multiple patches). + """ + # Type checks + assert isinstance(domain, Domain) + assert isinstance(ddm, (DomainDecomposition, MultiPatchDomainDecomposition)) + assert isinstance(pdim, int) + assert isinstance(mappings, (NoneType, dict)) + + # Extract info from domain + ldim : int = domain.dim + interior_names : list = domain.interior_names + + # Check sanity of pdim + assert pdim >= ldim + + # Check that the domain decomposition matches the domain, and get the + # decomposition, number of cells, and periodicity of each patch + if len(domain) == 1: + assert isinstance(ddm, DomainDecomposition) + patch_ddms = [ddm] + patch_ncells = [ddm.ncells] + patch_periods = [ddm.periods] + else: + assert isinstance(ddm, MultiPatchDomainDecomposition) + assert len(ddm.domains) == len(domain) + patch_ddms = ddm.domains + patch_ncells = ddm.ncells + patch_periods = ddm.periods + + assert all(len(n) == ldim for n in patch_ncells) + + # Store the number of cells and the periodicity as read-only tuples + ncells = {name: tuple(int(ni) for ni in n) + for name, n in zip(interior_names, patch_ncells)} + periodic = {name: tuple(bool(pi) for pi in p) + for name, p in zip(interior_names, patch_periods)} + + # Check sanity of mappings. A spline mapping must be defined on the + # domain decomposition of its patch + if mappings is None: + mappings = {name : None for name in interior_names} + else: + assert set(mappings.keys()) == set(interior_names) + assert all(isinstance(m, (BasicCallableMapping, NoneType)) for m in mappings.values()) + assert all(m.pdim == pdim for m in mappings.values() if m is not None) + assert all(mappings[name].space.domain_decomposition is patch_ddm + for name, patch_ddm in zip(interior_names, patch_ddms) + if isinstance(mappings[name], SplineMapping)) + + # Many functions access the mappings by position: use the order of + # the domain interiors + mappings = {name: mappings[name] for name in interior_names} + + # Add attributes to the new object + self._domain = domain + self._ldim = ldim + self._pdim = pdim + self._ncells = ncells + self._mappings = mappings + self._periodic = periodic + self._comm = ddm.comm + self._ddm = ddm + self._cart = None + #-------------------------------------------------------------------------- - # Option [1]: from a (domain, mappings) or a file + # Option [1]: from a file #-------------------------------------------------------------------------- - def __init__(self, domain=None, ncells=None, periodic=None, mappings=None, - filename=None, comm=None, mpi_dims_mask=None): - - # ... read the geometry if the filename is given - if filename is not None: - self.read(filename, comm=comm, mpi_dims_mask=mpi_dims_mask) - - elif domain is not None: - assert isinstance(domain, Domain) - assert isinstance(ncells, dict) - assert isinstance(mappings, dict) - if periodic is not None: - assert isinstance(periodic, dict) - - # ... check sanity - interior_names = domain.interior_names - mappings_keys = sorted(list(mappings.keys())) - - assert sorted(interior_names) == mappings_keys - # ... - - if periodic is None: - periodic = {patch: [False]*len(ncells_i) for patch, ncells_i in ncells.items()} - - self._domain = domain - self._ldim = domain.dim - self._pdim = domain.dim # TODO must be given => only dim is defined for a Domain - self._ncells = ncells - self._periodic = periodic - self._mappings = mappings - self._cart = None - self._is_parallel = comm is not None - - if len(domain) == 1: - #name = domain.name - name = interior_names[0] - self._ddm = DomainDecomposition(ncells[name], periodic[name], comm=comm, mpi_dims_mask=mpi_dims_mask) - else: - ncells = [ncells[itr] for itr in interior_names] - periodic = [periodic[itr] for itr in interior_names] - self._ddm = MultiPatchDomainDecomposition(ncells, periodic, comm=comm) + @classmethod + def from_file(cls, + filename : str, + *, + domain : Domain = None, + comm : MPI.Intracomm = None, + mpi_dims_mask : Iterable[bool] = None): - else: - raise ValueError('Wrong input') - # ... + """ + Create a Geometry instance from an HDF5 input file in Psydac format. - self._comm = comm + Parameters + ---------- + filename: str + The path to the geometry file. + + domain: sympde.topology.Domain, optional + The topological domain of the geometry. It must match the domain + defined in the file: same dimension, patch names, patches and + interfaces. Each patch must have a non-analytical mapping, to which + the corresponding spline mapping is attached; hence the same domain + should not be used with different geometry files. If not given, the + domain is created from the file. + + comm: MPI.Intracomm, optional + The MPI intra-communicator. + + mpi_dims_mask: Iterable[bool], optional + True if the dimension is to be used in the domain decomposition + (=default for each dimension). If mpi_dims_mask[i]=False, the i-th + dimension will not be decomposed. + + Returns + ------- + Geometry + The new instance. + + Raises + ------ + ValueError + If the given domain does not match the domain defined in the file, + or if one of its patches has no mapping or an analytical mapping. + """ + # Check the extension of the file + _, ext = os.path.splitext(filename) + if ext != '.h5': + raise ValueError('> Only h5 files are supported') + + # Read the topological domain, or check the given one against it + file_domain = Domain.from_file(filename) + if domain is None: + geometry_domain = file_domain + else: + _check_domain_matches_file(domain, file_domain, filename) + geometry_domain = domain + + parallel = comm is not None and comm.size > 1 + kwargs = dict(driver='mpio', comm=comm) if parallel else {} + with h5py.File(filename, mode='r', **kwargs) as h5: + yml = yaml.load(h5['geometry.yml'][()], Loader=yaml.SafeLoader) + if len(yml['patches']) == 0: + raise ValueError("Input file contains no patches.") + + # Pair the patches with the domain interiors by name, because the + # orders may differ (e.g. 'patch_10' is sorted before 'patch_2') + yml['patches'] = _sort_patches(yml['patches'], geometry_domain) + spaces = _read_patch_spaces(h5, yml) + names = geometry_domain.interior_names + ncells = {name: [W.ncells for W in s] for name, s in zip(names, spaces)} + periodic = {name: [W.periodic for W in s] for name, s in zip(names, spaces)} + + ddm = _make_ddm(geometry_domain, ncells, periodic, + comm=comm, mpi_dims_mask=mpi_dims_mask) + mappings = _read_spline_mappings(h5, yml, geometry_domain, ddm, spaces) + + # Add the spline mappings to the symbolic mappings of the domain + for itr in _get_interiors(geometry_domain): + itr.mapping.set_callable_mapping(mappings[itr.name]) + + return cls(geometry_domain, ddm=ddm, pdim=yml['pdim'], mappings=mappings) #-------------------------------------------------------------------------- # Option [2]: from a discrete mapping #-------------------------------------------------------------------------- @classmethod - def from_discrete_mapping(cls, mapping, *, comm=None, mpi_dims_mask=None, name=None): - """Create a geometry from one discrete mapping. + def from_discrete_mapping(cls, mapping, *, name=None): + """ + Create a single-patch Geometry instance from one discrete mapping. + + The geometry uses the domain decomposition of the mapping's spline + space, hence its distribution across the MPI processes is that of the + mapping. Parameters ---------- mapping : SplineMapping - The Mapping from the unit square to the physical domain. + The mapping from the unit square to the physical domain. - comm : MPI.Comm - MPI intra-communicator. - - mpi_dims_mask: list of bool - True if the dimension is to be used in the domain decomposition (=default for each dimension). - If mpi_dims_mask[i]=False, the i-th dimension will not be decomposed. - - name : string - Optional name for the Mapping that will be created. - Needed to avoid conflicts in case several mappings are created + name : str + Optional name for the symbolic Mapping that will be created. + Needed to avoid conflicts in case several mappings are created. + + Returns + ------- + Geometry + The new instance. """ mapping_name = name if name else 'mapping' - dim = mapping.ldim - M = Mapping(mapping_name, dim = dim) + dim = mapping.ldim + M = Mapping(mapping_name, ldim=dim, pdim=mapping.pdim) # this is a symbolic mapping domain = M(NCube(name = 'Omega', dim = dim, min_coords = [0.] * dim, max_coords = [1.] * dim)) M.set_callable_mapping(mapping) - mappings = {domain.name: mapping} - ncells = {domain.name: mapping.space.domain_decomposition.ncells} - periodic = {domain.name: mapping.space.domain_decomposition.periods} - - return Geometry(domain=domain, ncells=ncells, periodic=periodic, mappings=mappings, comm=comm, mpi_dims_mask=mpi_dims_mask) + return Geometry(domain = domain, + ddm = mapping.space.domain_decomposition, + pdim = mapping.pdim, + mappings = {domain.name: mapping}) #-------------------------------------------------------------------------- # Option [3]: discrete topological line/square/cube #-------------------------------------------------------------------------- @classmethod def from_topological_domain(cls, domain, ncells, *, periodic=None, comm=None, mpi_dims_mask=None): + """ + Create a Geometry instance from a topological domain, without discrete mappings. + + The domain may have symbolic mappings (e.g. analytical mappings), and + callable mappings attached to them are kept (e.g. the discrete mappings + of a domain returned by `Geometry.from_file`), but no discrete mappings + are created. All the patches of the domain must be `NCubeInterior` + objects of the same dimension, i.e. lines, squares, or cubes. The + physical dimension is given by the mappings of the patches if they + have any, and otherwise it equals the logical dimension. + + Parameters + ---------- + domain : sympde.topology.Domain + The single- or multi-patch topological domain, with or without + symbolic mappings. + + ncells : Iterable[int] | dict[str, Iterable[int]] + The number of cells along each direction. A list or tuple applies + to all patches; a dictionary gives the number of cells of each + patch, keyed by the names of the domain interiors. + + periodic : Iterable[bool] | dict[str, Iterable[bool]], optional + The periodicity along each direction, given in the same way as + `ncells` (default: no periodicity). A warning is issued if a + multipatch domain is periodic along any direction. + + comm : MPI.Intracomm, optional + The MPI intra-communicator. + + mpi_dims_mask : Iterable[bool], optional + True if the dimension is to be used in the domain decomposition + (=default for each dimension). If mpi_dims_mask[i]=False, the i-th + dimension will not be decomposed. Only used for a single patch. + + Returns + ------- + Geometry + The new instance. + + Raises + ------ + TypeError + If a patch of the domain is not an `NCubeInterior`. + + ValueError + If the patches of the domain have different physical dimensions. + """ + assert isinstance(domain, Domain) + interior = domain.interior if not isinstance(interior, Union): interior = [interior] for itr in interior: if not isinstance(itr, NCubeInterior): - msg = "Topological domain must be an NCube;"\ + msg = "The topological domain of each patch must be an NCube;"\ " got {} instead.".format(type(itr)) raise TypeError(msg) - mappings = {itr.name:None for itr in interior} + mappings = {itr.name : None for itr in interior} + + # The physical dimension is given by the mappings of the patches, if any + pdims = {itr.mapping.pdim if itr.mapping is not None else itr.dim for itr in interior} + if len(pdims) > 1: + raise ValueError(f'The patches of the domain have different physical dimensions: {pdims}') + pdim = pdims.pop() if isinstance(ncells, (list, tuple)): - ncells = {itr.name:ncells for itr in interior} + ncells = {itr.name : ncells for itr in interior} if periodic is None: - periodic = [False]*domain.dim + periodic = [False] * domain.dim else: if len(interior) > 1 and True in periodic: + import warnings msg = "Discretizing a multipatch domain with a periodic flag is not advised -- continue at your own risk." # [MCP 18.12.2025] the following line may be causing a strange error in the CI (MPI tests for macos-14/Python 3.10) # warnings.warn(msg, Warning) warnings.warn(msg, UserWarning) - if isinstance(periodic, (list, tuple)): - periodic = {itr.name:periodic for itr in interior} + periodic = {itr.name : periodic for itr in interior} - geo = Geometry(domain=domain, mappings=mappings, ncells=ncells, periodic=periodic, comm=comm, mpi_dims_mask=mpi_dims_mask) + ddm = _make_ddm(domain, ncells, periodic, + comm=comm, mpi_dims_mask=mpi_dims_mask) - return geo + return Geometry(domain = domain, + ddm = ddm, + pdim = pdim, + mappings = mappings) #-------------------------------------------------------------------------- @property @@ -227,10 +411,6 @@ def domain(self): def ddm(self): return self._ddm - @property - def is_parallel(self): - return self._is_parallel - @property def mappings(self): return self._mappings @@ -238,167 +418,6 @@ def mappings(self): def __len__(self): return len(self.domain) - def read(self, filename, comm=None, mpi_dims_mask=None): - # ... check extension of the file - basename, ext = os.path.splitext(filename) - if not(ext == '.h5'): - raise ValueError('> Only h5 files are supported') - # ... - - # read the topological domain - domain = Domain.from_file(filename) - connectivity = construct_connectivity(domain) - - if len(domain)==1: - interiors = [domain.interior] - else: - interiors = list(domain.interior.args) - - if not(comm is None): - kwargs = dict( driver='mpio', comm=comm ) if comm.size > 1 else {} - - else: - kwargs = {} - - h5 = h5py.File( filename, mode='r', **kwargs ) - yml = yaml.load( h5['geometry.yml'][()], Loader=yaml.SafeLoader ) - - ldim = yml['ldim'] - pdim = yml['pdim'] - - n_patches = len( yml['patches'] ) - - # ... - if n_patches == 0: - - h5.close() - raise ValueError( "Input file contains no patches." ) - # ... - - # ... read patchs - mappings = {} - ncells = {} - periodic = {} - spaces = [None]*n_patches - for i_patch in range( n_patches ): - - item = yml['patches'][i_patch] - patch_name = item['name'] - mapping_id = item['mapping_id'] - dtype = item['type'] - patch = h5[mapping_id] - if dtype in ['SplineMapping', 'NurbsMapping']: - - degree = [int (p) for p in patch.attrs['degree' ]] - periodic_i = [bool(b) for b in patch.attrs['periodic']] - knots = [patch['knots_{}'.format(d)][:] for d in range( ldim )] - space_i = [SplineSpace( degree=p, knots=k, periodic=P ) - for p,k,P in zip( degree, knots, periodic_i )] - - spaces[i_patch] = space_i - - ncells [interiors[i_patch].name] = [sp.ncells for sp in space_i] - periodic[interiors[i_patch].name] = periodic_i - - self._cart = None - if n_patches == 1: - self._ddm = DomainDecomposition(ncells[domain.name], periodic[domain.name], comm=comm, mpi_dims_mask=mpi_dims_mask) - ddms = [self._ddm] - else: - ncells_ = [ncells[itr.name] for itr in interiors] - periodic = [periodic[itr.name] for itr in interiors] - self._ddm = MultiPatchDomainDecomposition(ncells_, periodic, comm=comm) - ddms = self._ddm.domains - - carts = create_cart(ddms, spaces) - g_spaces = {inter:TensorFemSpace( ddms[i], *spaces[i], cart=carts[i]) for i,inter in enumerate(interiors)} - - for i,j in connectivity: - ((axis_i, ext_i), (axis_j , ext_j)) = connectivity[i, j] - minus = interiors[i] - plus = interiors[j] - max_ncells = [max(ni,nj) for ni,nj in zip(ncells[minus.name],ncells[plus.name])] - g_spaces[minus].add_refined_space(ncells=max_ncells) - g_spaces[plus].add_refined_space(ncells=max_ncells) - - # ... construct interface spaces - construct_interface_spaces(self._ddm, g_spaces, carts, interiors, connectivity) - - for i_patch in range( n_patches ): - - item = yml['patches'][i_patch] - patch_name = item['name'] - mapping_id = item['mapping_id'] - dtype = item['type'] - patch = h5[mapping_id] - space_i = spaces[i_patch] - if dtype in ['SplineMapping', 'NurbsMapping']: - tensor_space = g_spaces[interiors[i_patch]] - - if dtype == 'SplineMapping': - mapping = SplineMapping.from_control_points( tensor_space, - patch['points'][..., :pdim] ) - - elif dtype == 'NurbsMapping': - mapping = NurbsMapping.from_control_points_weights( tensor_space, - patch['points'][..., :pdim], - patch['weights'] ) - - mapping.set_name( item['name'] ) - mappings[patch_name] = mapping - - if n_patches>1: - coeffs = [[e._coeffs for e in mapping._fields] for mapping in mappings.values()] - spaces = [[coeffs_ij.space for coeffs_ij in coeffs_i] for coeffs_i in coeffs] - spaces = [BlockVectorSpace(*space) for space in spaces] - w_spaces = [sp.spaces[0] for sp in spaces] - space = BlockVectorSpace(*spaces, connectivity=connectivity) - w_space = BlockVectorSpace(*w_spaces, connectivity=connectivity) - v = BlockVector(space) - w = BlockVector(w_space) - mapping_list = list(mappings.values()) - for i in range(n_patches): - for j in range(len(coeffs[i])): - v[i][j] = coeffs[i][j] - - mapping = mapping_list[i] - if isinstance(mapping, NurbsMapping): - w[i] = mapping.weights_field.coeffs - else: - w[i] = v[i][0].space.zeros() - - v.update_ghost_regions() - w.update_ghost_regions() - - else: - mapping = list(mappings.values())[0] - for f in mapping._fields: - f.coeffs.update_ghost_regions() - - if isinstance(mapping, NurbsMapping): - mapping.weights_field.coeffs.update_ghost_regions() - - - # ... close the h5 file - h5.close() - # ... - - # Add spline callable mappings to domain undefined mappings - # NOTE: We assume that interiors and mappings.values() use the same ordering - for patch, F in zip(interiors, mappings.values()): - patch.mapping.set_callable_mapping(F) - - # ... - self._ldim = ldim - self._pdim = pdim - self._mappings = mappings - self._domain = domain - self._comm = comm - self._ncells = ncells - self._periodic = periodic - self._is_parallel = comm is not None - # ... - def export( self, filename ): """ Parameters @@ -417,9 +436,13 @@ def export( self, filename ): yml['ldim'] = self.ldim yml['pdim'] = self.pdim + # Only spline mappings can be written to file + missing = [name for name, mapping in self.mappings.items() + if not isinstance(mapping, SplineMapping)] + if missing: + raise ValueError(f'Cannot export patches without a spline mapping: {missing}') + # ... information about the patches - if not( self.mappings ): - raise ValueError('No mappings were found') patches_info = [] i_mapping = 0 @@ -437,7 +460,6 @@ def export( self, filename ): yml['patches'] = patches_info # ... - # ... topology topo_yml = self.domain.todict() # ... @@ -883,3 +905,274 @@ def _read_patch(lines, i_patch, n_lines_per_patch, list_begin_line): nrb = NURBS(knots, control=points, weights=W) return nrb + +#============================================================================== +def _get_interiors(domain : Domain) -> list: + """Return the list of interior patches of a (multi-patch) domain.""" + return [domain.interior] if len(domain) == 1 else list(domain.interior.args) + +def _sort_patches(patches : list[dict], domain : Domain) -> list[dict]: + """ + Sort the patches of a geometry file in the order of the domain interiors. + + Each patch is identified by its name, which is either the name of the + interior (e.g. 'mapping_0(patch_0)', as written by `Geometry.export`), or + the name of its logical domain (e.g. 'patch_0', as written by + `export_multipatch_nurbs_to_hdf5`). + + Parameters + ---------- + patches : list[dict] + The patches listed in the `geometry.yml` section of the file. + + domain : sympde.topology.Domain + The topological domain read from the file. + + Returns + ------- + list[dict] + The same patches, in the order of the domain interiors. + + Raises + ------ + ValueError + If the patch of an interior cannot be identified uniquely. + """ + by_name = {patch['name']: patch for patch in patches} + if len(by_name) != len(patches): + raise ValueError('The patches of the geometry file must have distinct names') + + sorted_patches = [] + for itr in _get_interiors(domain): + names = {itr.name} if itr.logical_domain is None else {itr.name, itr.logical_domain.name} + matches = [by_name[name] for name in names if name in by_name] + if len(matches) != 1: + raise ValueError(f'Cannot identify the patch of interior {itr.name} in the geometry file') + sorted_patches.append(matches[0]) + + return sorted_patches + +def _get_interface_keys(domain : Domain) -> set: + """Return the interfaces of a domain as hashable tuples, with orientation.""" + interfaces = domain.interfaces + if interfaces is None: + return set() + if not isinstance(interfaces, Union): + interfaces = [interfaces] + return {(i.minus.domain.name, i.minus.axis, i.minus.ext, + i.plus .domain.name, i.plus .axis, i.plus .ext, i.ornt) for i in interfaces} + +def _check_domain_matches_file(domain : Domain, file_domain : Domain, filename : str): + """Raise a ValueError if the domain cannot be used with the geometry file.""" + + # Properties are compared one by one, because domain == file_domain raises + # an error for multi-patch domains (see pyccel/sympde#197) + + # Analytical mappings would be used in assembly instead of the spline ones + for patch in _get_interiors(domain): + if patch.mapping is None or patch.mapping.is_analytical: + raise ValueError(f"Patch {patch.name} must have a non-analytical " + f"mapping to be used with geometry file {filename}") + + if domain.dim != file_domain.dim: + raise ValueError(f"Domain dimension {domain.dim} does not match " + f"dimension {file_domain.dim} in file {filename}") + + if domain.interior_names != file_domain.interior_names: + raise ValueError(f"Patch names {domain.interior_names} do not match " + f"{file_domain.interior_names} in file {filename}") + + # Patch equality ignores the parametric bounds (see pyccel/sympde#198) + for patch, file_patch in zip(_get_interiors(domain), _get_interiors(file_domain)): + bounds = (patch .logical_domain.min_coords, patch .logical_domain.max_coords) + file_bounds = (file_patch.logical_domain.min_coords, file_patch.logical_domain.max_coords) + if bounds != file_bounds: + raise ValueError(f"Parametric bounds {bounds} of patch {patch.name} " + f"do not match {file_bounds} in file {filename}") + + if _get_interface_keys(domain) != _get_interface_keys(file_domain): + raise ValueError(f"Interfaces of the domain do not match those in file {filename}") + +#============================================================================== +def _make_ddm(domain : Domain, + ncells : dict[str, Iterable[int]], + periodic : dict[str, Iterable[bool]] | None = None, + comm : MPI.Intracomm | None = None, + mpi_dims_mask : Iterable[bool] | None = None, + ) -> DomainDecomposition | MultiPatchDomainDecomposition: + """ + Create the domain decomposition of a single- or multi-patch domain. + + Parameters + ---------- + domain : Sympde.topology.Domain + The symbolic topological domain to be decomposed. + + ncells : dict[str, Iterable[int]] + The number of cells of each patch along each direction. + + periodic : dict[str, Iterable[bool]], optional + The periodicity of each patch along each direction (default: False). + + comm : MPI.Intracomm, optional + The MPI intra-communicator. + + mpi_dims_mask : Iterable[bool], optional + True if the dimension is to be used in the domain decomposition + (=default for each dimension). Only used for a single patch. + + Returns + ------- + DomainDecomposition | MultiPatchDomainDecomposition + The decomposition of a single patch, or of multiple patches. + """ + assert isinstance(domain, Domain) + assert isinstance(ncells, dict) + assert isinstance(periodic, (NoneType, dict)) + assert isinstance(comm, (NoneType, MPI.Intracomm)) + assert isinstance(mpi_dims_mask, (NoneType, Iterable)) + + ldim = domain.dim + interior_names = domain.interior_names + + # Check sanity of ncells + assert set(ncells.keys()) == set(interior_names) + assert all(len(n) == ldim for n in ncells.values()) + assert all(isinstance(ni, (int, np.integer)) for ni in chain(*ncells.values())) + assert all(ni > 0 for ni in chain(*ncells.values())) + + # Check sanity of periodic (no periodicity by default) + if periodic is None: + periods = {name: (False,) * ldim for name in interior_names} + else: + periods = periodic + assert set(periods.keys()) == set(interior_names) + assert all(len(p) == ldim for p in periods.values()) + assert all(isinstance(pi, bool) for pi in chain(*periods.values())) + + # Check sanity of mpi_dims_mask + if mpi_dims_mask is not None: + assert len(mpi_dims_mask) == ldim + assert all(isinstance(mask, bool) for mask in mpi_dims_mask) + + if len(domain) == 1: + name = interior_names[0] + return DomainDecomposition(ncells[name], periods[name], comm=comm, + mpi_dims_mask=mpi_dims_mask) + + return MultiPatchDomainDecomposition([ncells [name] for name in interior_names], + [periods[name] for name in interior_names], + comm=comm) + +#============================================================================== +def _read_patch_spaces(h5 : h5py.File, yml : dict) -> list[list[SplineSpace]]: + """ + Read the 1D spline spaces of each patch from an open geometry file. + + Parameters + ---------- + h5 : h5py.File + The open geometry file. + + yml : dict + The content of its `geometry.yml` section. + + Returns + ------- + list[list[SplineSpace]] + The 1D spline spaces along each direction, for each patch. + """ + spaces = [] + for item in yml['patches']: + if item['type'] not in ('SplineMapping', 'NurbsMapping'): + raise ValueError(f"Mapping type {item['type']} of patch {item['name']} " + "is not supported") + patch = h5[item['mapping_id']] + degree = [int (p) for p in patch.attrs['degree' ]] + periodic = [bool(b) for b in patch.attrs['periodic']] + knots = [patch[f'knots_{d}'][:] for d in range(yml['ldim'])] + spaces.append([SplineSpace(degree=p, knots=k, periodic=P) + for p, k, P in zip(degree, knots, periodic)]) + return spaces + +def _read_spline_mappings(h5 : h5py.File, + yml : dict, + domain : Domain, + ddm : DomainDecomposition | MultiPatchDomainDecomposition, + spaces : list[list[SplineSpace]], + ) -> dict[str, SplineMapping]: + """ + Create the spline mappings of all patches from an open geometry file, on + the given domain decomposition. + + Parameters + ---------- + h5 : h5py.File + The open geometry file. + + yml : dict + The content of its `geometry.yml` section. + + domain : Sympde.topology.Domain + The topological domain, whose interiors follow the order of the patches. + + ddm : DomainDecomposition | MultiPatchDomainDecomposition + The decomposition of the domain. + + spaces : list[list[SplineSpace]] + The 1D spline spaces along each direction, for each patch. + + Returns + ------- + dict[str, SplineMapping] + The spline mappings, keyed by the names of the domain interiors. + """ + interiors = _get_interiors(domain) + connectivity = construct_connectivity(domain) + ddms = [ddm] if len(domain) == 1 else ddm.domains + pdim = yml['pdim'] + + carts = create_cart(ddms, spaces) + g_spaces = {itr: TensorFemSpace(ddms[i], *spaces[i], cart=carts[i]) + for i, itr in enumerate(interiors)} + + # Add refined spaces on both sides of each interface + for i, j in connectivity: + max_ncells = [max(Wi.ncells, Wj.ncells) for Wi, Wj in zip(spaces[i], spaces[j])] + g_spaces[interiors[i]].add_refined_space(ncells=max_ncells) + g_spaces[interiors[j]].add_refined_space(ncells=max_ncells) + + construct_interface_spaces(ddm, g_spaces, carts, interiors, connectivity) + + mappings = {} + for itr, item in zip(interiors, yml['patches']): + patch = h5[item['mapping_id']] + points = patch['points'][..., :pdim] + if item['type'] == 'SplineMapping': + mapping = SplineMapping.from_control_points(g_spaces[itr], points) + else: + mapping = NurbsMapping.from_control_points_weights(g_spaces[itr], points, + patch['weights']) + mapping.set_name(item['name']) + mappings[itr.name] = mapping + + # Update the ghost regions within each patch and across interfaces + if len(domain) > 1: + # One block vector per component, with one block per patch which + # references the coefficients of the mapping + mapping_list = list(mappings.values()) + components = [[m.fields[k].coeffs for m in mapping_list] for k in range(pdim)] + if any(isinstance(m, NurbsMapping) for m in mapping_list): + components.append([m.weights_field.coeffs if isinstance(m, NurbsMapping) + else m.fields[0].coeffs.space.zeros() for m in mapping_list]) + for blocks in components: + V = BlockVectorSpace(*[b.space for b in blocks], connectivity=connectivity) + BlockVector(V, blocks=blocks).update_ghost_regions() + else: + mapping, = mappings.values() + for f in mapping.fields: + f.coeffs.update_ghost_regions() + if isinstance(mapping, NurbsMapping): + mapping.weights_field.coeffs.update_ghost_regions() + + return mappings diff --git a/psydac/cad/mesh/multipatch/magnet.h5 b/psydac/cad/mesh/multipatch/magnet.h5 index c18c688ec..9e497c2a0 100644 Binary files a/psydac/cad/mesh/multipatch/magnet.h5 and b/psydac/cad/mesh/multipatch/magnet.h5 differ diff --git a/psydac/cad/mesh/multipatch/plate_with_hole_mp.h5 b/psydac/cad/mesh/multipatch/plate_with_hole_mp.h5 index 168f193ce..4c7b26c90 100644 Binary files a/psydac/cad/mesh/multipatch/plate_with_hole_mp.h5 and b/psydac/cad/mesh/multipatch/plate_with_hole_mp.h5 differ diff --git a/psydac/cad/mesh/multipatch/plate_with_hole_mp_6.h5 b/psydac/cad/mesh/multipatch/plate_with_hole_mp_6.h5 index f6bee5b7f..6e6cf7456 100644 Binary files a/psydac/cad/mesh/multipatch/plate_with_hole_mp_6.h5 and b/psydac/cad/mesh/multipatch/plate_with_hole_mp_6.h5 differ diff --git a/psydac/cad/mesh/multipatch/plate_with_hole_mp_7.h5 b/psydac/cad/mesh/multipatch/plate_with_hole_mp_7.h5 index cf56bfcad..fca6c4265 100644 Binary files a/psydac/cad/mesh/multipatch/plate_with_hole_mp_7.h5 and b/psydac/cad/mesh/multipatch/plate_with_hole_mp_7.h5 differ diff --git a/psydac/cad/mesh/multipatch/square.h5 b/psydac/cad/mesh/multipatch/square.h5 index ae3f87a27..781c05feb 100644 Binary files a/psydac/cad/mesh/multipatch/square.h5 and b/psydac/cad/mesh/multipatch/square.h5 differ diff --git a/psydac/cad/mesh/multipatch/square_repeated_knots.h5 b/psydac/cad/mesh/multipatch/square_repeated_knots.h5 index 471b157f3..a17343290 100644 Binary files a/psydac/cad/mesh/multipatch/square_repeated_knots.h5 and b/psydac/cad/mesh/multipatch/square_repeated_knots.h5 differ diff --git a/psydac/cad/tests/test_cad.py b/psydac/cad/tests/test_cad.py new file mode 100644 index 000000000..9340c2995 --- /dev/null +++ b/psydac/cad/tests/test_cad.py @@ -0,0 +1,193 @@ +#---------------------------------------------------------------------------# +# This file is part of PSYDAC which is released under MIT License. See the # +# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE # +# for full license details. # +#---------------------------------------------------------------------------# +from itertools import product + +import numpy as np +import pytest +from mpi4py import MPI + +from psydac.cad.cad import elevate, refine, translate +from psydac.cad.gallery import annulus, circle, quart_circle +from psydac.cad.tests.test_gallery import make_nurbs_mapping +from psydac.ddm.cart import DomainDecomposition +from psydac.fem.splines import SplineSpace +from psydac.fem.tensor import TensorFemSpace +from psydac.mapping.discrete import NurbsMapping, SplineMapping +from psydac.mapping.discrete_gallery import discrete_mapping + +def make_curve(*, nurbs: bool = False, pdim: int = 1) -> SplineMapping: + """Create a nonlinear spline or NURBS curve in `pdim` dimensions.""" + W = SplineSpace(degree=2, knots=[0.0, 0.0, 0.0, 0.3, 0.7, 1.0, 1.0, 1.0]) + V = TensorFemSpace(DomainDecomposition([W.ncells], [False]), W) + points = np.array([[0.0, 0.0, 1.0], + [0.2, 0.5, 0.8], + [0.9, 0.7, 0.3], + [1.5, 0.4, -0.2], + [2.0, -0.1, 0.0]])[:, :pdim] + if nurbs: + weights = np.array([1.0, 0.5, 2.0, 0.7, 1.0]) + return NurbsMapping.from_control_points_weights(V, points, weights) + return SplineMapping.from_control_points(V, points) + +def make_surface(*, pdim: int = 2) -> NurbsMapping: + """Create a NURBS surface in `pdim` dimensions from the quarter annulus, + lifted to z = x * y if pdim = 3.""" + degrees, knots, points, weights = quart_circle() + if pdim == 3: + z = points[..., 0] * points[..., 1] + points = np.concatenate([points, z[..., None]], axis=-1) + return make_nurbs_mapping(degrees, knots, points, weights) + +def make_volume() -> SplineMapping: + """Create a nonlinear spline volume in 3D.""" + return discrete_mapping('collela', ncells=[2, 2, 2], degree=[2, 2, 2]) + +mappings_1d = pytest.mark.parametrize('make_mapping', [ + make_curve, + lambda: make_curve(nurbs=True), + lambda: make_curve(pdim=2), + lambda: make_curve(nurbs=True, pdim=3), +], ids=['spline', 'nurbs', 'spline_curve_2d', 'nurbs_curve_3d']) + +# A spline mapping, NURBS mappings whose weights vary along axis 0 +# (quart_circle), axis 1 (annulus), or both (circle), and a NURBS surface +mappings_2d = pytest.mark.parametrize('make_mapping', [ + lambda: discrete_mapping('collela', ncells=[4, 4], degree=[2, 2]), + lambda: make_nurbs_mapping(*quart_circle()), + lambda: make_nurbs_mapping(*annulus()), + lambda: make_nurbs_mapping(*circle()), + lambda: make_surface(pdim=3), +], ids=['spline', 'quart_circle', 'annulus', 'circle', 'surface_3d']) + +#============================================================================== +def assert_same_geometry(F, G) -> None: + """Check that two mappings agree on a uniform grid of the logical domain of G.""" + eta = list(product(np.linspace(0.0, 1.0, 9), repeat=G.ldim)) + assert np.allclose([F(*e) for e in eta], [G(*e) for e in eta], rtol=0, atol=1e-14) + +def check_translate(F) -> None: + """Check that translation moves the geometry by a fixed displacement.""" + displ = np.array([1.0, -2.0, 0.5])[:F.pdim] + G = translate(F, displ) + + assert type(G) is type(F) + assert_same_geometry(lambda *eta: np.asarray(F(*eta)) + displ, G) + +def check_elevate(F, axis: int) -> None: + """Check that degree elevation along an axis preserves the geometry.""" + G = elevate(F, axis=axis, times=1) + + expected_degree = list(F.space.degree) + expected_degree[axis] += 1 + assert type(G) is type(F) + assert list(G.space.degree) == expected_degree + assert_same_geometry(F, G) + +def check_refine(F, axis: int, values: tuple[float, ...] = (0.1, 0.35, 0.6, 0.9)) -> None: + """Check that knot insertion along an axis preserves the geometry.""" + # The default values avoid the double knots of annulus, which already have multiplicity p + G = refine(F, axis=axis, values=values) + + expected_breaks = np.union1d(F.space.spaces[axis].breaks, values) + assert type(G) is type(F) + assert np.allclose(G.space.spaces[axis].breaks, expected_breaks, rtol=0, atol=1e-15) + assert_same_geometry(F, G) + +#============================================================================== +@mappings_1d +def test_translate_1d(make_mapping) -> None: + check_translate(make_mapping()) + +@mappings_2d +def test_translate_2d(make_mapping) -> None: + check_translate(make_mapping()) + +#============================================================================== +@mappings_1d +def test_elevate_1d(make_mapping) -> None: + check_elevate(make_mapping(), axis=0) + +@mappings_2d +@pytest.mark.parametrize('axis', [0, 1]) +def test_elevate_2d(axis: int, make_mapping) -> None: + check_elevate(make_mapping(), axis) + +@pytest.mark.parametrize('axis', [0, 1, 2]) +def test_elevate_3d(axis: int) -> None: + check_elevate(make_volume(), axis) + +#============================================================================== +@mappings_1d +def test_refine_1d(make_mapping) -> None: + check_refine(make_mapping(), axis=0) + +# Values which only increase the multiplicity of a knot create no new cells +@pytest.mark.parametrize('values', [(0.3,), (0.5, 0.5)], ids=['existing_knot', 'repeated_value']) +def test_refine_multiplicity(values: tuple[float, ...]) -> None: + check_refine(make_curve(nurbs=True), axis=0, values=values) + +@mappings_2d +@pytest.mark.parametrize('axis', [0, 1]) +def test_refine_2d(axis: int, make_mapping) -> None: + check_refine(make_mapping(), axis) + +@pytest.mark.parametrize('axis', [0, 1, 2]) +def test_refine_3d(axis: int) -> None: + check_refine(make_volume(), axis) + +#============================================================================== +def make_mapping_parallel(comm: MPI.Comm, mpi_dims_mask: list[bool] | None = None): + return discrete_mapping('collela', ncells=[8, 8], degree=[2, 2], comm=comm, + mpi_dims_mask=mpi_dims_mask) + +def assert_same_local_coeffs(F, G) -> None: + """Check that the local data (owned and ghost coefficients) of a distributed + mapping G agrees with the corresponding part of a serial mapping F.""" + for f, g in zip(F.fields, G.fields): + V = g.coeffs.space + everything = (slice(None),) * V.ndim + local_part = tuple(slice(s - m*p, e + m*p + 1) + for s, e, p, m in zip(V.starts, V.ends, V.pads, V.shifts)) + assert np.allclose(g.coeffs[everything], f.coeffs[local_part], + rtol=0, atol=1e-14) + +#============================================================================== +@pytest.mark.mpi +@pytest.mark.parametrize('axis', [0, 1]) +def test_elevate_parallel(axis: int) -> None: + + F = make_mapping_parallel(MPI.COMM_WORLD) + F_serial = make_mapping_parallel(MPI.COMM_SELF) + + G = elevate(F, axis=axis, times=1) + G_serial = elevate(F_serial, axis=axis, times=1) + + assert_same_local_coeffs(G_serial, G) + +#============================================================================== +# The shifted values move the subdomain boundaries by more than p cells +@pytest.mark.mpi +@pytest.mark.parametrize('axis', [0, 1]) +@pytest.mark.parametrize('values', [ + [0.1, 0.35, 0.6, 0.9], + [0.01, 0.04, 0.07, 0.1, 0.13, 0.16, 0.19, 0.22], +], ids=['spread', 'shifted']) +@pytest.mark.parametrize('mpi_dims_mask', [None, [False, True], [True, False]]) +def test_refine_parallel(mpi_dims_mask: list[bool] | None, values: list, axis: int) -> None: + + F = make_mapping_parallel(MPI.COMM_WORLD, mpi_dims_mask) + F_serial = make_mapping_parallel(MPI.COMM_SELF) + + G = refine(F, axis=axis, values=values) + G_serial = refine(F_serial, axis=axis, values=values) + + assert_same_local_coeffs(G_serial, G) + + # The directions excluded by the mask are not decomposed (see #622) + ddm = G.space.domain_decomposition + assert ddm.mpi_dims_mask == F.space.domain_decomposition.mpi_dims_mask + if mpi_dims_mask is not None: + assert all(n == 1 for n, use_dim in zip(ddm.nprocs, mpi_dims_mask) if not use_dim) diff --git a/psydac/cad/tests/test_gallery.py b/psydac/cad/tests/test_gallery.py new file mode 100644 index 000000000..65ae359ea --- /dev/null +++ b/psydac/cad/tests/test_gallery.py @@ -0,0 +1,41 @@ +#---------------------------------------------------------------------------# +# This file is part of PSYDAC which is released under MIT License. See the # +# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE # +# for full license details. # +#---------------------------------------------------------------------------# +import numpy as np +import pytest + +from psydac.cad.gallery import annulus, circle, quart_circle +from psydac.ddm.cart import DomainDecomposition +from psydac.fem.splines import SplineSpace +from psydac.fem.tensor import TensorFemSpace +from psydac.mapping.discrete import NurbsMapping + +#============================================================================== +def make_nurbs_mapping(degrees, knots, points, weights) -> NurbsMapping: + """Create a NURBS mapping from the output of a gallery function.""" + spaces = [SplineSpace(knots=k, degree=p) for k, p in zip(knots, degrees)] + ncells = [len(V.breaks) - 1 for V in spaces] + domain_decomposition = DomainDecomposition(ncells, [False] * len(spaces)) + space = TensorFemSpace(domain_decomposition, *spaces) + return NurbsMapping.from_control_points_weights(space, points, weights) + +#============================================================================== +# Each circular edge of the logical square is given as (axis, value, radius) +@pytest.mark.parametrize('center', [None, (1.0, -2.0)]) +@pytest.mark.parametrize(('gallery_function', 'kwargs', 'circular_edges'), [ + (quart_circle, dict(rmin=0.5, rmax=1.0), [(1, 0.0, 0.5), (1, 1.0, 1.0)]), + (annulus, dict(rmin=0.5, rmax=1.0), [(0, 0.0, 0.5), (0, 1.0, 1.0)]), + (circle, dict(radius=2.0), [(0, 0.0, 2.0), (0, 1.0, 2.0), (1, 0.0, 2.0), (1, 1.0, 2.0)]), +], ids=['quart_circle', 'annulus', 'circle']) +def test_gallery_circular_edges(gallery_function, kwargs: dict, circular_edges: list, center) -> None: + + F = make_nurbs_mapping(*gallery_function(**kwargs, center=center)) + x0 = np.zeros(2) if center is None else np.array(center) + + t = np.linspace(0.0, 1.0, 17) + for axis, value, radius in circular_edges: + eta = [(value, ti) if axis == 0 else (ti, value) for ti in t] + r = [np.linalg.norm(F(*e) - x0) for e in eta] + assert np.allclose(r, radius, rtol=0, atol=1e-14) diff --git a/psydac/cad/tests/test_geometry.py b/psydac/cad/tests/test_geometry.py index 6671d1870..2a73935bd 100644 --- a/psydac/cad/tests/test_geometry.py +++ b/psydac/cad/tests/test_geometry.py @@ -4,67 +4,128 @@ # for full license details. # #---------------------------------------------------------------------------# import os +from itertools import product import pytest import numpy as np from mpi4py import MPI -from sympde.topology import Domain, Line, Square, Cube, Mapping +from sympde.topology import Domain, Line, Square, Cube, Mapping, IdentityMapping from psydac.cad.geometry import Geometry, export_nurbs_to_hdf5, refine_nurbs from psydac.cad.geometry import import_geopdes_to_nurbs from psydac.cad.cad import elevate, refine -from psydac.cad.gallery import quart_circle, circle +from psydac.cad.gallery import quart_circle +from psydac.cad.multipatch import export_multipatch_nurbs_to_hdf5 +from psydac.cad.tests.test_cad import assert_same_geometry, make_curve, make_surface from psydac.mapping.discrete import SplineMapping, NurbsMapping from psydac.mapping.discrete_gallery import discrete_mapping from psydac.fem.splines import SplineSpace from psydac.fem.tensor import TensorFemSpace from psydac.utilities.utils import refine_array_1d -from psydac.ddm.cart import DomainDecomposition +from psydac.ddm.cart import (DomainDecomposition, + MultiPatchDomainDecomposition) -base_dir = os.path.dirname(os.path.realpath(__file__)) -#============================================================================== -def test_geometry_2d_1(): +import psydac.cad.mesh as mesh_mod - ncells = [1,1] - degree = [2,2] - # create an identity mapping - mapping = discrete_mapping('identity', ncells=ncells, degree=degree) - - # create a topological domain - F = Mapping('F', dim=2) - domain = F(Square(name='Omega')) - - # associate the mapping to the topological domain - mappings = {domain.name: mapping} - - # Define ncells as a dict - ncells = {domain.name:ncells} +base_dir = os.path.dirname(os.path.realpath(__file__)) - # create a geometry from a topological domain and the dict of mappings - geo = Geometry(domain=domain, ncells=ncells, mappings=mappings) +# Unit patch of each logical dimension +PATCHES = {1: Line, 2: Square} +#============================================================================== +def _assert_geometries_equal(geo: Geometry, ref: Geometry) -> None: + """Check that two geometries are equal, comparing their patches by position.""" + assert geo.ldim == ref.ldim + assert geo.pdim == ref.pdim + assert list(geo.ncells.values()) == list(ref.ncells.values()) + assert [list(p) for p in geo.periodic.values()] == [list(p) for p in ref.periodic.values()] + assert len(geo.mappings) == len(ref.mappings) + + for F, F_ref in zip(geo.mappings.values(), ref.mappings.values()): + assert type(F) is type(F_ref) + assert list(F.space.degree) == list(F_ref.space.degree) + assert all(np.array_equal(k, k_ref) for k, k_ref in zip(F.space.knots, F_ref.space.knots)) + + fields = [*F .fields, F .weights_field] if isinstance(F, NurbsMapping) else F .fields + fields_ref = [*F_ref.fields, F_ref.weights_field] if isinstance(F, NurbsMapping) else F_ref.fields + assert all(np.array_equal(f.coeffs.toarray(), f_ref.coeffs.toarray()) + for f, f_ref in zip(fields, fields_ref)) + +def _check_round_trips(geo: Geometry, mapping: SplineMapping, tmp_path) -> Geometry: + """Check that a single-patch geometry is preserved by export/read and + by from_discrete_mapping, and return the geometry read from file.""" + filename = str(tmp_path / 'geo.h5') + geo.export(filename) + geo_read = Geometry.from_file(filename) + _assert_geometries_equal(geo_read, geo) + + # A geometry read from file is exported again without changes + filename_again = str(tmp_path / 'geo_again.h5') + geo_read.export(filename_again) + _assert_geometries_equal(Geometry.from_file(filename_again), geo) + + _assert_geometries_equal(Geometry.from_discrete_mapping(mapping), geo) + + return geo_read + +def _make_identity_mapping(ddm: DomainDecomposition, degree: list[int], shift: float = 0.0) -> SplineMapping: + """Create a spline mapping on the decomposition, equal to the identity + shifted along x. Its control points are given by the Greville abscissae.""" + spaces = [SplineSpace(degree=p, grid=np.linspace(0.0, 1.0, n + 1)) + for p, n in zip(degree, ddm.ncells)] + grevilles = [W.greville for W in spaces] + grevilles[0] = grevilles[0] + shift + points = np.stack(np.meshgrid(*grevilles, indexing='ij'), axis=-1) + return SplineMapping.from_control_points(TensorFemSpace(ddm, *spaces), points) - # export the geometry - geo.export('geo.h5') +#============================================================================== +@pytest.mark.parametrize('ldim', [1, 2]) +def test_geometry_identity_spline(ldim: int, tmp_path) -> None: - # read it again - geo_0 = Geometry(filename='geo.h5') + ddm = DomainDecomposition([2, 3][:ldim], [False] * ldim) + mapping = _make_identity_mapping(ddm, degree=[2] * ldim) + domain = Mapping('F', dim=ldim)(PATCHES[ldim](name='Omega')) - # export it again - geo_0.export('geo_0.h5') + geo = Geometry(domain, ddm=ddm, pdim=ldim, mappings={domain.name: mapping}) + geo_read = _check_round_trips(geo, mapping, tmp_path) - # create a geometry from a discrete mapping - geo_1 = Geometry.from_discrete_mapping(mapping) + # The mapping read from file is the identity + F_read, = geo_read.mappings.values() + eta = list(product(np.linspace(0.0, 1.0, 5), repeat=ldim)) + assert np.allclose([F_read(*e) for e in eta], eta, rtol=0, atol=1e-15) - # export it - geo_1.export('geo_1.h5') +#============================================================================== +@pytest.mark.parametrize('make_mapping', [ + lambda: make_curve(nurbs=True), + lambda: make_curve(nurbs=True, pdim=2), + lambda: make_curve(nurbs=True, pdim=3), + lambda: make_surface(pdim=3), +], ids=['line', 'curve_2d', 'curve_3d', 'surface_3d']) +def test_geometry_nurbs_manifold(make_mapping, tmp_path) -> None: + + mapping = elevate(make_mapping(), axis=0, times=1) + mapping = refine(mapping, axis=0, values=[0.1, 0.5]) + ldim, pdim = mapping.ldim, mapping.pdim + domain = Mapping('F', ldim=ldim, pdim=pdim)(PATCHES[ldim](name='Omega')) + + geo = Geometry(domain, ddm=mapping.space.domain_decomposition, pdim=pdim, + mappings={domain.name: mapping}) + geo_read = _check_round_trips(geo, mapping, tmp_path) + + # The mapping read from file is the original one + F_read, = geo_read.mappings.values() + assert_same_geometry(F_read, mapping) + + # The symbolic mapping has the physical dimension of the discrete one + assert Geometry.from_discrete_mapping(mapping).domain.mapping.pdim == pdim #============================================================================== -def test_geometry_2d_2(): +def test_geometry_nurbs_quarter_annulus(tmp_path) -> None: # create a nurbs mapping - degrees, knots, points, weights = quart_circle( rmin=0.5, rmax=1.0, center=None ) + rmin, rmax = 0.5, 1.0 + degrees, knots, points, weights = quart_circle( rmin=rmin, rmax=rmax, center=None ) # Create tensor spline space, distributed spaces = [SplineSpace( knots=k, degree=p ) for k,p in zip(knots, degrees)] @@ -86,98 +147,21 @@ def test_geometry_2d_2(): # associate the mapping to the topological domain mappings = {domain.name: mapping} - # Define ncells as a dict - ncells = {domain.name:[len(space.breaks)-1 for space in mapping.space.spaces]} - - periodic = {domain.name:[space.periodic for space in mapping.space.spaces]} - # create a geometry from a topological domain and the dict of mappings - geo = Geometry(domain=domain, ncells=ncells, periodic=periodic, mappings=mappings) - - # export the geometry - geo.export('quart_circle.h5') + geo = Geometry(domain, ddm=mapping.space.domain_decomposition, pdim=2, + mappings=mappings) - # read it again - geo_0 = Geometry(filename='quart_circle.h5') + geo_read = _check_round_trips(geo, mapping, tmp_path) - # export it again - geo_0.export('quart_circle_0.h5') - - # create a geometry from a discrete mapping - geo_1 = Geometry.from_discrete_mapping(mapping) - - # export it - geo_1.export('quart_circle_1.h5') - -#============================================================================== -# TODO to be removed -def test_geometry_2d_3(): - - # create a nurbs mapping - degrees, knots, points, weights = quart_circle( rmin=0.5, rmax=1.0, center=None ) - - # Create tensor spline space, distributed - spaces = [SplineSpace( knots=k, degree=p ) for k,p in zip(knots, degrees)] - ncells = [len(space.breaks)-1 for space in spaces] - domain_decomposition = DomainDecomposition(ncells=ncells, periods=[False]*2, comm=None) - - space = TensorFemSpace( domain_decomposition, *spaces ) - - mapping = NurbsMapping.from_control_points_weights( space, points, weights ) - - mapping = elevate( mapping, axis=1, times=1 ) - - n = 8 - t = np.linspace(0, 1, n+1)[1:-1] - - # TODO allow for 1d numpy array - t = list(t) - - for axis in [0, 1]: - mapping = refine( mapping, axis=axis, values=t ) - - # create a geometry from a discrete mapping - geo = Geometry.from_discrete_mapping(mapping) - - # export it - geo.export('quart_circle.h5') - -#============================================================================== -# TODO to be removed -def test_geometry_2d_4(): - - # create a nurbs mapping - radius = np.sqrt(2)/2. - degrees, knots, points, weights = circle( radius=radius, center=None ) - - # Create tensor spline space, distributed - spaces = [SplineSpace( knots=k, degree=p ) for k,p in zip(knots, degrees)] - ncells = [len(space.breaks)-1 for space in spaces] - domain_decomposition = DomainDecomposition(ncells=ncells, periods=[False]*2, comm=None) - - space = TensorFemSpace( domain_decomposition, *spaces ) - - mapping = NurbsMapping.from_control_points_weights( space, points, weights ) - - n = 8 -# n = 32 - t = np.linspace(0, 1, n+1)[1:-1] - - # TODO allow for 1d numpy array - t = list(t) - - for axis in [0, 1]: - mapping = refine( mapping, axis=axis, values=t ) - - # create a geometry from a discrete mapping - geo = Geometry.from_discrete_mapping(mapping) - - # export it - geo.export('circle.h5') + # The mapping read from file is a quarter annulus, radial along axis 1 + F_read = list(geo_read.mappings.values())[0] + t = np.linspace(0.0, 1.0, 9) + assert np.allclose([np.hypot(*F_read(e1, 0.0)) for e1 in t], rmin, rtol=0, atol=1e-14) + assert np.allclose([np.hypot(*F_read(e1, 1.0)) for e1 in t], rmax, rtol=0, atol=1e-14) #============================================================================== @pytest.mark.mpi -def test_geometry_with_mpi_dims_mask(): +def test_geometry_with_mpi_dims_mask(mpi_tmp_path): comm = MPI.COMM_WORLD rank = comm.rank @@ -199,27 +183,20 @@ def test_geometry_with_mpi_dims_mask(): # associate the mapping to the topological domain mappings = {domain.name: mapping} - # Define d_ncells as a dict - d_ncells = {domain.name: ncells} - # Create a geometry from a topological domain and the dict of mappings # Here we allow for any distribution of the domain: mpi_dims_mask is not passed - geo = Geometry(domain=domain, ncells=d_ncells, mappings=mappings, comm=comm) - geo.export('geo_mpi_dims.h5') + geo = Geometry(domain, ddm=mapping.space.domain_decomposition, pdim=3, + mappings=mappings) + filename = str(mpi_tmp_path / 'geo_mpi_dims.h5') + geo.export(filename) # Read geometry file in parallel, but using mpi_dims_mask - geo_from_file = Geometry(filename='geo_mpi_dims.h5', comm=comm, mpi_dims_mask=mpi_dims_mask) + geo_from_file = Geometry.from_file(filename=filename, comm=comm, mpi_dims_mask=mpi_dims_mask) # Verify that the domain is distributed as expected assert geo_from_file.ddm.starts == expected_starts assert geo_from_file.ddm.ends == expected_ends - # Safely remove the file - comm.Barrier() - if rank == 0: - os.remove('geo_mpi_dims.h5') - - # ============================================================================== @pytest.mark.mpi def test_from_discrete_mapping(): @@ -227,18 +204,20 @@ def test_from_discrete_mapping(): comm = MPI.COMM_WORLD rank = comm.rank size = comm.size - mpi_dims_mask = [False, False, True] # We swill verify that this has an effect + mpi_dims_mask = [False, False, True] # We will verify that this has an effect ncells = [4, 8, 2 * size] # Each process should have two cells along x3 - degree = [3, 3, 3] + degree = [3, 3, 2] # Each process must own at least `degree` cells expected_starts = (0, 0, 2 * rank) expected_ends = (3, 7, 2 * rank + 1) - # Create a mapping - mapping = discrete_mapping('identity', ncells=ncells, degree=degree) + # Create a mapping using mpi_dims_mask + mapping = discrete_mapping('identity', ncells=ncells, degree=degree, + mpi_dims_mask=mpi_dims_mask) - # Create geometry from the mapping using mpi_dims_mask - geo_from_mapping = Geometry.from_discrete_mapping(mapping, comm=comm, mpi_dims_mask=mpi_dims_mask) + # The geometry uses the domain decomposition of the mapping + geo_from_mapping = Geometry.from_discrete_mapping(mapping) + assert geo_from_mapping.ddm is mapping.space.domain_decomposition # Verify that the domain is distributed as expected assert geo_from_mapping.ddm.starts == expected_starts @@ -268,10 +247,188 @@ def test_from_topological_domain(): assert geo_from_domain.ddm.starts == expected_starts assert geo_from_domain.ddm.ends == expected_ends +# ============================================================================== +def _join_squares(pdims: tuple[int, int]) -> Domain: + """Join two unit squares, with mappings of the given physical dimensions.""" + A = Mapping('FA', ldim=2, pdim=pdims[0])(Square('A')) + B = Mapping('FB', ldim=2, pdim=pdims[1])(Square('B')) + return Domain.join([A, B], [((0, 0, 1), (1, 0, -1), 1)], 'Omega') + +@pytest.mark.parametrize(('make_domain_2d', 'pdim'), [ + (lambda: Square('Omega'), 2), + (lambda: Mapping('F', dim=2)(Square('Omega')), 2), + (lambda: Mapping('F', ldim=2, pdim=3)(Square('Omega')), 3), + (lambda: _join_squares((3, 3)), 3), +], ids=['no-mapping', 'mapping', 'surface-3d', 'multipatch-surface-3d']) +def test_from_topological_domain_pdim(make_domain_2d, pdim: int) -> None: + + # The physical dimension is given by the mappings of the patches, if any + geo = Geometry.from_topological_domain(make_domain_2d(), [3, 5]) + + assert geo.ldim == 2 + assert geo.pdim == pdim + +def test_from_topological_domain_different_pdims() -> None: + + with pytest.raises(ValueError, match='different physical dimensions'): + Geometry.from_topological_domain(_join_squares((2, 3)), [3, 5]) + +@pytest.mark.parametrize('npatches', [1, 2]) +def test_geometry_export_without_spline_mappings(npatches: int, tmp_path) -> None: + + geo = Geometry.from_topological_domain(_make_domain(npatches), [3, 5]) + filename = tmp_path / 'geo.h5' + + # The error is raised before the file is created + with pytest.raises(ValueError, match='without a spline mapping'): + geo.export(str(filename)) + assert not filename.exists() + +# ============================================================================== +@pytest.mark.parametrize('ldim', [1, 2]) +@pytest.mark.parametrize('npatches', [1, 2]) +def test_geometry_init_without_mappings(npatches: int, ldim: int) -> None: + + domain = _make_domain(npatches, ldim=ldim) + + ncells = [3, 5][:ldim] + if npatches == 1: + ddm = DomainDecomposition(ncells, [False] * ldim) + else: + ddm = MultiPatchDomainDecomposition([ncells] * npatches, + [[False] * ldim] * npatches) + + geo = Geometry(domain, ddm=ddm, pdim=ldim) + + # The number of cells and the periodicity are obtained from the decomposition + assert geo.ncells == {name: tuple(ncells) for name in domain.interior_names} + assert geo.periodic == {name: (False,) * ldim for name in domain.interior_names} + + assert geo.mappings == {name: None for name in domain.interior_names} + +#============================================================================== +def _export_strip(filename: str, ncells_x: list[int]) -> None: + """ + Export a strip of bilinear patches glued along x, where patch k maps the + unit square onto [k, k+1] x [0, 1] with ncells_x[k] cells along x. + """ + from igakit.nurbs import NURBS + + nurbs = [] + for k, nx in enumerate(ncells_x): + points = np.zeros((nx + 1, 2, 2)) + points[..., 0] = np.linspace(k, k + 1, nx + 1)[:, None] + points[..., 1] = [0, 1] + knots_x = np.concatenate([[0.], np.linspace(0, 1, nx + 1), [1.]]) + nurbs.append(NURBS([knots_x, [0, 0, 1, 1]], points)) + connectivity = {(k, k + 1): ((0, 1), (0, -1)) for k in range(len(ncells_x) - 1)} + + export_multipatch_nurbs_to_hdf5(filename, nurbs, connectivity) + +# SymPDE sorts the interiors by name (patch_0, patch_1, patch_10, patch_2, ...), +# while the file keeps the order patch_0, ..., patch_10 +@pytest.mark.parametrize('ncells_x', [[1] * 11, [1] * 10 + [2]], ids=['same-sizes', 'patch_10-refined']) +def test_geometry_from_file_more_than_10_patches(ncells_x: list[int], tmp_path) -> None: + + filename = str(tmp_path / 'strip.h5') + _export_strip(filename, ncells_x) + geo = Geometry.from_file(filename) + + # Each patch must be paired with its interior by name, not by position + for itr in geo.domain.interior.args: + k = int(itr.logical_domain.name.removeprefix('patch_')) + x, _ = geo.mappings[itr.name](0.5, 0.5) + assert np.isclose(x, k + 0.5), f'interior {itr.name} got the mapping of patch_{round(x - 0.5)}' + assert geo.ncells[itr.name] == (ncells_x[k], 1) + +#============================================================================== +def test_geometry_from_file_multipatch() -> None: + + mesh_dir = os.path.dirname(mesh_mod.__file__) + filename = os.path.join(mesh_dir, 'multipatch', 'magnet.h5') + geo = Geometry.from_file(filename) + names = geo.domain.interior_names + + # Patch information is keyed by the interior names, in the same order + assert list(geo.ncells) == names + assert list(geo.periodic) == names + assert list(geo.mappings) == names + assert all(isinstance(n, tuple) and len(n) == 2 for n in geo.ncells.values()) + assert all(p == (False, False) for p in geo.periodic.values()) + + # Each spline mapping is defined on the decomposition of its patch + for F, patch_ddm in zip(geo.mappings.values(), geo.ddm.domains): + assert F.space.domain_decomposition is patch_ddm + +#============================================================================== +def _make_domain(npatches: int, *, ldim: int = 2, ornt: int = 1) -> Domain: + """Create a domain made of one or two unit lines or squares, with generic mappings.""" + Patch = PATCHES[ldim] + if npatches == 1: + return Mapping('G', dim=ldim)(Patch('P')) + + A = Mapping('GA', dim=ldim)(Patch('PA')) + B = Mapping('GB', dim=ldim)(Patch('PB')) + # 1D interfaces have no orientation + return Domain.join([A, B], [((0, 0, 1), (1, 0, -1), ornt if ldim > 1 else None)], 'Omega') + +def _export_geometry(domain: Domain, filename: str) -> None: + """Export a spline geometry on the domain, with patches side by side along x.""" + names = domain.interior_names + ldim = domain.dim + ncells = [2, 3][:ldim] + if len(domain) == 1: + ddm = DomainDecomposition(ncells, [False] * ldim) + patch_ddms = [ddm] + else: + ddm = MultiPatchDomainDecomposition([ncells] * len(names), + [[False] * ldim] * len(names)) + patch_ddms = ddm.domains + + # Spline mappings must be defined on the decomposition of their patch + mappings = {name: _make_identity_mapping(patch_ddm, degree=[2] * ldim, shift=i) + for i, (name, patch_ddm) in enumerate(zip(names, patch_ddms))} + + Geometry(domain, ddm=ddm, pdim=ldim, mappings=mappings).export(filename) + +#============================================================================== +@pytest.mark.parametrize('ldim', [1, 2]) +@pytest.mark.parametrize('npatches', [1, 2]) +def test_geometry_from_file_with_domain(npatches: int, ldim: int, tmp_path) -> None: + + domain = _make_domain(npatches, ldim=ldim) + filename = str(tmp_path / 'geo.h5') + _export_geometry(domain, filename) + + geo = Geometry.from_file(filename, domain=domain) + + # The given domain is kept, and the spline mappings are attached to it + assert geo.domain is domain + patches = [domain.interior] if npatches == 1 else domain.interior.args + for patch, F in zip(patches, geo.mappings.values()): + assert patch.mapping.get_callable_mapping() is F + +#============================================================================== +@pytest.mark.parametrize(('npatches', 'make_wrong_domain', 'message'), [ + (1, lambda: IdentityMapping('G', dim=2)(Square('P')), 'non-analytical'), + (1, lambda: Square('P'), 'non-analytical'), + (1, lambda: Mapping('G', dim=3)(Cube('P')), 'dimension'), + (1, lambda: Mapping('H', dim=2)(Square('P')), 'Patch names'), + (1, lambda: Mapping('G', dim=2)(Square('P', bounds1=(0, 2))), 'Parametric bounds'), + (2, lambda: _make_domain(2, ornt=-1), 'Interfaces'), +], ids=['analytical', 'no-mapping', 'dimension', 'names', 'bounds', 'orientation']) +def test_geometry_from_file_with_wrong_domain(npatches: int, make_wrong_domain, message: str, tmp_path) -> None: + + filename = str(tmp_path / 'geo.h5') + _export_geometry(_make_domain(npatches), filename) + + with pytest.raises(ValueError, match=message): + Geometry.from_file(filename, domain=make_wrong_domain()) + #============================================================================== @pytest.mark.parametrize( 'ncells', [[8,8], [12,12], [14,14]] ) @pytest.mark.parametrize( 'degree', [[2,2], [3,2], [2,3], [3,3], [4,4]] ) -def test_export_nurbs_to_hdf5(ncells, degree): +def test_export_nurbs_to_hdf5(ncells, degree, tmp_path): # create pipe geometry from igakit.cad import circle, ruled, bilinear, join @@ -284,11 +441,11 @@ def test_export_nurbs_to_hdf5(ncells, degree): # refine the nurbs object new_pipe = refine_nurbs(pipe, ncells=ncells, degree=degree) - filename = "pipe.h5" + filename = str(tmp_path / "pipe.h5") export_nurbs_to_hdf5(filename, new_pipe) # read the geometry - geo = Geometry(filename=filename) + geo = Geometry.from_file(filename) domain = geo.domain min_coords = domain.logical_domain.min_coords @@ -300,7 +457,7 @@ def test_export_nurbs_to_hdf5(ncells, degree): assert abs(max_coords[0] - pipe.breaks(0)[-1])<1e-15 assert abs(max_coords[1] - pipe.breaks(1)[-1])<1e-15 - mapping = geo.mappings[domain.logical_domain.name] + mapping = geo.mappings[domain.interior.name] assert isinstance(mapping, NurbsMapping) @@ -324,8 +481,7 @@ def test_export_nurbs_to_hdf5(ncells, degree): #============================================================================== @pytest.mark.parametrize( 'ncells', [[8,8], [12,12], [14,14]] ) @pytest.mark.parametrize( 'degree', [[2,2], [3,2], [2,3], [3,3], [4,4]] ) -def test_import_geopdes_to_nurbs(ncells, degree): - +def test_import_geopdes_to_nurbs(ncells, degree, tmp_path): filename = os.path.join(base_dir, "geo_Lshaped_C1.txt") L_shaped = import_geopdes_to_nurbs(filename) @@ -333,11 +489,11 @@ def test_import_geopdes_to_nurbs(ncells, degree): # refine the nurbs object L_shaped = refine_nurbs(L_shaped, ncells=ncells, degree=degree) - filename = "L_shaped.h5" + filename = str(tmp_path / "L_shaped.h5") export_nurbs_to_hdf5(filename, L_shaped) # read the geometry - geo = Geometry(filename=filename) + geo = Geometry.from_file(filename) domain = geo.domain min_coords = domain.logical_domain.min_coords @@ -349,7 +505,7 @@ def test_import_geopdes_to_nurbs(ncells, degree): assert abs(max_coords[0] - L_shaped.breaks(0)[-1])<1e-15 assert abs(max_coords[1] - L_shaped.breaks(1)[-1])<1e-15 - mapping = geo.mappings[domain.logical_domain.name] + mapping = geo.mappings[domain.interior.name] space = mapping.space knots = space.knots @@ -361,39 +517,14 @@ def test_import_geopdes_to_nurbs(ncells, degree): if isinstance(mapping, NurbsMapping): assert np.allclose(L_shaped.weights.flatten(), mapping._weights_field.coeffs.toarray(), 1e-15, 1e-15) -#============================================================================== -@pytest.mark.xfail -def test_geometry_1(): - - line = Geometry.as_line(ncells=[10]) - square = Geometry.as_square(ncells=[10, 10]) - cube = Geometry.as_cube(ncells=[10, 10, 10]) - #============================================================================== # CLEAN UP SYMPY NAMESPACE #============================================================================== def teardown_module(): - import os from sympy.core import cache cache.clear_cache() - # Remove HDF5 files generated by Geometry.export() - filenames = [ - 'geo.h5', - 'geo_0.h5', - 'geo_1.h5', - 'quart_circle.h5', - 'quart_circle_0.h5', - 'quart_circle_1.h5', - 'circle.h5', - 'pipe.h5', - 'L_shaped.h5', - ] - for fname in filenames: - if os.path.exists(fname): - os.remove(fname) - def teardown_function(): from sympy.core import cache cache.clear_cache() diff --git a/psydac/ddm/cart.py b/psydac/ddm/cart.py index a43bceb9b..79d5b5e95 100644 --- a/psydac/ddm/cart.py +++ b/psydac/ddm/cart.py @@ -214,6 +214,7 @@ def __init__(self, ncells, periods, comm=None, global_comm=None, num_threads=Non self._ncells = tuple ( ncells ) self._periods = tuple ( periods ) + self._mpi_dims_mask = None if mpi_dims_mask is None else tuple(mpi_dims_mask) self._comm = comm self._global_comm = comm if global_comm is None else global_comm self._comm_cart = comm @@ -300,6 +301,11 @@ def ncells( self ): def periods( self ): return self._periods + @property + def mpi_dims_mask(self) -> tuple[bool, ...] | None: + """The directions which may be decomposed, or None if not specified (all directions).""" + return self._mpi_dims_mask + @property def size( self ): return self._size @@ -398,9 +404,10 @@ def refine(self, ncells, global_element_starts, global_element_ends): assert len( ncells ) == len( self.ncells ) assert all(nc>=snc for nc, snc in zip(ncells, self.ncells)) + # The mask is needed to recover the same process grid (see #622) domain = DomainDecomposition(self.ncells, self.periods, comm=self.comm, global_comm=self.global_comm, num_threads=self.num_threads, - size=self.size) + size=self.size, mpi_dims_mask=self.mpi_dims_mask) domain._ncells = tuple ( ncells ) # Store arrays with all the starts and ends along each direction for every process @@ -412,7 +419,7 @@ def refine(self, ncells, global_element_starts, global_element_ends): domain._starts = tuple( domain._global_element_starts[axis][c] for axis,c in zip(range(self._ndims), self._coords) ) domain._ends = tuple( domain._global_element_ends [axis][c] for axis,c in zip(range(self._ndims), self._coords) ) - domain._local_ncells = tuple(e-s+1 for s,e in zip(self._starts, self._ends)) + domain._local_ncells = tuple(e-s+1 for s,e in zip(domain._starts, domain._ends)) return domain #================================================================================== diff --git a/psydac/ddm/tests/test_cart.py b/psydac/ddm/tests/test_cart.py new file mode 100644 index 000000000..6187225e7 --- /dev/null +++ b/psydac/ddm/tests/test_cart.py @@ -0,0 +1,43 @@ +#---------------------------------------------------------------------------# +# This file is part of PSYDAC which is released under MIT License. See the # +# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE # +# for full license details. # +#---------------------------------------------------------------------------# +import pytest +from mpi4py import MPI + +from psydac.ddm.cart import DomainDecomposition + +masks = pytest.mark.parametrize('mpi_dims_mask', [None, [False, True], [True, False]]) + +#============================================================================== +@masks +def test_DomainDecomposition_mpi_dims_mask(mpi_dims_mask) -> None: + + ddm = DomainDecomposition([3, 5], [False, False], mpi_dims_mask=mpi_dims_mask) + + expected = None if mpi_dims_mask is None else tuple(mpi_dims_mask) + assert ddm.mpi_dims_mask == expected + +#============================================================================== +@pytest.mark.mpi +@masks +def test_DomainDecomposition_refine(mpi_dims_mask) -> None: + + # Each direction can be decomposed among all processes + comm = MPI.COMM_WORLD + ncells = [2 * comm.size, 3 * comm.size] + ddm = DomainDecomposition(ncells, [False, False], comm=comm, mpi_dims_mask=mpi_dims_mask) + + # Split each cell into 2 x 2 cells, keeping the subdomain boundaries + global_starts = [2 * s for s in ddm.global_element_starts] + global_ends = [2 * e + 1 for e in ddm.global_element_ends] + new_ddm = ddm.refine([2 * n for n in ncells], global_starts, global_ends) + + # The refined decomposition has the same process grid (see #622) + assert new_ddm.mpi_dims_mask == ddm.mpi_dims_mask + assert list(new_ddm.nprocs) == list(ddm.nprocs) + assert list(new_ddm.comm_cart.Get_topo()[0]) == list(ddm.nprocs) + assert new_ddm.starts == tuple(2 * s for s in ddm.starts) + assert new_ddm.ends == tuple(2 * e + 1 for e in ddm.ends) + assert new_ddm.local_ncells == tuple(2 * n for n in ddm.local_ncells) # see #623 diff --git a/psydac/feec/polar/examples/poisson_2d.py b/psydac/feec/polar/examples/poisson_2d.py index 688f821d1..1f1558d62 100644 --- a/psydac/feec/polar/examples/poisson_2d.py +++ b/psydac/feec/polar/examples/poisson_2d.py @@ -19,18 +19,14 @@ `mpirun -n 6 python poisson_2d.py -S -d 3 3 -t disk -D 0.2 -m 'C0conga'` """ -from dataclasses import dataclass -from time import sleep, time +# Imports are local to keep the CLI fast (e.g. --help) +# pylint: disable=import-outside-toplevel -import numpy as np -from mpi4py import MPI -from pyccel import lambdify -from sympy import Rational, cos, pi, sin, symbols +from dataclasses import dataclass +# Base classes are needed at class definition from psydac.feec.polar.examples.polar_model_2d import PolarModel2D from psydac.linalg.basic import LinearOperator -from psydac.linalg.stencil import StencilMatrix, StencilVector -from psydac.utilities.operators import Laplacian # ============================================================================== @@ -56,6 +52,8 @@ class Poisson2D(PolarModel2D): """ def __init__(self, domain_log, mapping, phi_log, rho_log): + from pyccel import lambdify + super().__init__(domain_log, mapping) self._domain_log = domain_log @@ -93,8 +91,10 @@ def disk(R, shift_D): and source term. """ + import numpy as np from sympde.topology.analytical_mapping import TargetMapping from sympde.topology.domain import Square + from sympy import cos, pi, sin, symbols domain_log = Square("Omega", bounds1=(0, R), bounds2=(0, 2 * np.pi)) params = dict(c1=shift_D * R * R, c2=0, k=0, D=shift_D) @@ -134,8 +134,12 @@ def target(): $\phi(x,y) = (1 - s^8)\sin(k_x(x - 0.5))\cos(k_y y)$. """ + import numpy as np from sympde.topology.analytical_mapping import TargetMapping from sympde.topology.domain import Square + from sympy import Rational, cos, pi, sin + + from psydac.utilities.operators import Laplacian domain_log = Square("Omega", bounds1=(0, 1), bounds2=(0, 2 * np.pi)) params = dict(c1=0, c2=0, k=Rational(3, 10), D=Rational(2, 10)) @@ -167,8 +171,12 @@ def czarny(): $\phi(x,y) = (1 - s^8)\sin(\pi x)\cos(\pi y)$. """ + import numpy as np from sympde.topology.analytical_mapping import CzarnyMapping from sympde.topology.domain import Square + from sympy import Rational, cos, pi, sin + + from psydac.utilities.operators import Laplacian domain_log = Square("Omega", bounds1=(0, 1), bounds2=(0, 2 * np.pi)) params = dict(c1=0, c2=0, eps=Rational(1, 5), b=Rational(7, 5)) @@ -240,6 +248,7 @@ def __init__(self, S, M, P, alpha): C0PolarProjection_V0, C1PolarProjection_U0, ) + from psydac.linalg.stencil import StencilMatrix assert isinstance(S, StencilMatrix) assert isinstance(M, StencilMatrix) @@ -257,6 +266,8 @@ def __init__(self, S, M, P, alpha): self.W0 = W0 def dot(self, x, out=None): + from psydac.linalg.stencil import StencilVector + if out is None: y = self.M._domain.zeros() else: @@ -339,7 +350,7 @@ def compute_errors(phi, phi_ref, M, S): ErrorDiagnostics Reference L2 and H1 norms and relative L2 and H1 errors. """ - + import numpy as np # L2 and H1 norms ref_l2_2 = M.dot_inner(phi_ref.coeffs, phi_ref.coeffs) @@ -386,6 +397,8 @@ def plot_solution(use_spline_mapping, model, ncells, periodic, V0_h, refine=10): Refinement factor for plotting. """ import matplotlib.pyplot as plt + import numpy as np + from mpi4py import MPI from psydac.cad.geometry import Geometry from psydac.ddm.cart import DomainDecomposition @@ -394,7 +407,7 @@ def plot_solution(use_spline_mapping, model, ncells, periodic, V0_h, refine=10): from psydac.utilities.utils import refine_array_1d if use_spline_mapping: - geometry = Geometry(filename="geo.h5", comm=MPI.COMM_SELF) + geometry = Geometry.from_file("geo.h5", comm=MPI.COMM_SELF) map_discrete = [*geometry.mappings.values()].pop() Vnew = map_discrete.space map_plot = map_discrete @@ -488,7 +501,10 @@ def run_poisson_2d( mpi_comm, # given by function 'parallel_run_from_cli' ): import os + from time import sleep, time + import numpy as np + from mpi4py import MPI from sympde.calculus import dot, grad from sympde.expr import BilinearForm, LinearForm, integral from sympde.topology import ScalarFunctionSpace, elements_of diff --git a/psydac/feec/polar/examples/polar_model_2d.py b/psydac/feec/polar/examples/polar_model_2d.py index 40cd24dae..8e5742ada 100644 --- a/psydac/feec/polar/examples/polar_model_2d.py +++ b/psydac/feec/polar/examples/polar_model_2d.py @@ -138,10 +138,7 @@ def build_geometry( # In order to create a sympde.Domain object from this mapping we have # to create first a HDF5 file and then load as sympde.Domain.fromfile - geometry = Geometry.from_discrete_mapping( - map_discrete, - comm=mpi_comm, - ) + geometry = Geometry.from_discrete_mapping(map_discrete) geometry.export(filename) self._geometry_export_time = time() - t0 diff --git a/psydac/feec/tests/test_feec_conf_projectors_cart_2d.py b/psydac/feec/tests/test_feec_conf_projectors_cart_2d.py index cf64ecf27..b327278dd 100644 --- a/psydac/feec/tests/test_feec_conf_projectors_cart_2d.py +++ b/psydac/feec/tests/test_feec_conf_projectors_cart_2d.py @@ -4,7 +4,6 @@ # for full license details. # #---------------------------------------------------------------------------# import pytest -from collections import OrderedDict import numpy as np from sympy import Tuple, lambdify @@ -43,7 +42,7 @@ def get_polynomial_function(degree, hom_bc_axes, domain): g0_y = (y - 0.75)**degree[1] expr = g0_x * g0_y - callable_function = lambdify(domain.coordinates, expr) + callable_function = lambdify(domain.coordinates, expr, modules=['numpy']) return expr, callable_function diff --git a/psydac/fem/partitioning.py b/psydac/fem/partitioning.py index e9c587916..2ba87642b 100644 --- a/psydac/fem/partitioning.py +++ b/psydac/fem/partitioning.py @@ -67,11 +67,15 @@ def partition_coefficients(domain_decomposition, spaces, min_blocks=None): if min_blocks is None: min_blocks = [None] * ndims - for s, e, V, mb in zip(global_starts, global_ends, spaces, min_blocks): - if V.periodic or mb is None: - assert all(e-s+1 >= V.degree) - else: - assert all(e-s+1 >= mb) + # Each process must own enough coefficients along each axis. The global + # starts and ends are the same on all processes, hence so is the error + for axis, (s, e, V, mb) in enumerate(zip(global_starts, global_ends, spaces, min_blocks)): + min_coeffs = V.degree if V.periodic or mb is None else mb + owned_coeffs = int(min(e - s + 1)) + if owned_coeffs < min_coeffs: + raise ValueError(f'Along axis {axis}, an MPI process owns only {owned_coeffs} coefficients,' + f' but at least {min_coeffs} are needed: use fewer MPI processes along' + ' this axis (e.g. with mpi_dims_mask), or more cells') return global_starts, global_ends diff --git a/psydac/fem/tests/test_partitioning.py b/psydac/fem/tests/test_partitioning.py new file mode 100644 index 000000000..cfcdff92e --- /dev/null +++ b/psydac/fem/tests/test_partitioning.py @@ -0,0 +1,28 @@ +#---------------------------------------------------------------------------# +# This file is part of PSYDAC which is released under MIT License. See the # +# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE # +# for full license details. # +#---------------------------------------------------------------------------# +import numpy as np +import pytest +from mpi4py import MPI + +from psydac.ddm.cart import DomainDecomposition +from psydac.fem.splines import SplineSpace +from psydac.fem.tensor import TensorFemSpace + +#============================================================================== +# A single process owns all the coefficients, which are at least `degree` +@pytest.mark.mpi(min_size=2) +@pytest.mark.parametrize('periodic', [True, False]) +def test_partition_coefficients_too_few(periodic: bool) -> None: + + # Some processes own a single cell, hence less than `degree` coefficients + comm = MPI.COMM_WORLD + degree = 3 + ncells = max(comm.size, degree) + W = SplineSpace(degree=degree, grid=np.linspace(0, 1, ncells + 1), periodic=periodic) + ddm = DomainDecomposition([ncells], [periodic], comm=comm) + + with pytest.raises(ValueError, match=f'at least {degree} are needed'): + TensorFemSpace(ddm, W) diff --git a/psydac/fem/tests/test_plot_field_2d.py b/psydac/fem/tests/test_plot_field_2d.py index a099aee02..bf870bdf7 100644 --- a/psydac/fem/tests/test_plot_field_2d.py +++ b/psydac/fem/tests/test_plot_field_2d.py @@ -46,16 +46,16 @@ def test_plot_field(use_scalar_field, use_multipatch): degree = [2, 2] A = Square('A',bounds1=(0.5, 1.), bounds2=(0, np.pi/2)) - mapping_1 = PolarMapping('M1',2, c1= 0., c2= 0., rmin = 0., rmax=1.) + mapping_1 = PolarMapping('M1', 2, c1= 0., c2= 0., rmin = 0., rmax=1.) D1 = mapping_1(A) if use_multipatch: B = Square('B',bounds1=(0.5, 1.), bounds2=(np.pi/2, np.pi)) - mapping_2 = PolarMapping('M2',2, c1= 0., c2= 0., rmin = 0., rmax=1.) - D2 = mapping_2(B) + mapping_2 = PolarMapping('M2', 2, c1= 0., c2= 0., rmin = 0., rmax=1.) + D2 = mapping_2(B) - connectivity = [((0,1,1),(1,1,-1))] - patches = [D1,D2] + patches = [D1, D2] + connectivity = [((0, 1, 1), (1, 1,-1), 1)] domain = Domain.join(patches, connectivity, 'domain') else: domain = D1 diff --git a/psydac/linalg/kron.py b/psydac/linalg/kron.py index 045c6d015..f67691f41 100644 --- a/psydac/linalg/kron.py +++ b/psydac/linalg/kron.py @@ -654,8 +654,7 @@ def _inslice_to_temp(self, inslice, target): Copies data to an internal, 1-dimensional temporary array. Does not allocate any new array. """ - targetview = target[:self._localsize] - targetview.shape = inslice.shape + targetview = target[:self._localsize].reshape(inslice.shape, copy=False) targetview[:] = inslice @@ -664,11 +663,9 @@ def _reorder_temp_to_temp(self, source, target, i): Reorders the dimensions of the temporary arrays, and copies data from one to another. Does not allocate any new array. """ - sourceview = source[:self._localsize] - sourceview.shape = self._shapes[i] + sourceview = source[:self._localsize].reshape(self._shapes[i], copy=False) - targetview = target[:self._localsize] - targetview.shape = self._shapes[i+1] + targetview = target[:self._localsize].reshape(self._shapes[i+1], copy=False) targetview[:] = sourceview.transpose(self._perm) @@ -677,8 +674,7 @@ def _reorder_temp_to_outslice(self, source, outslice): Reorders the dimensions of the temporary array for a final time, and copies it to the output. Does not allocate any new array. """ - sourceview = source[:self._localsize] - sourceview.shape = self._shapes[-1] + sourceview = source[:self._localsize].reshape(self._shapes[-1], copy=False) outslice[:] = sourceview.transpose(self._perm) @@ -729,8 +725,8 @@ def solve_pass(self, workmem, tempmem): Ignored, it exists for compatibility with the parallel solver. """ # reshape necessary memory in column-major - view = workmem[:self._datasize] - view.shape = (self._numrhs,self._dimrhs) + view = workmem[:self._datasize].reshape( + (self._numrhs, self._dimrhs), copy=False) # call solver in in-place mode self._solver.solve(view, out=view) @@ -879,11 +875,12 @@ def _blocked_to_contiguous(self, blocked, contiguous): Copies from a blocked view to a contiguous view. Equals roughly a partial transpose, if the block sizes in the cartesian grid are the same. """ - blocked_view = blocked[:self._datasize] - blocked_view.shape = (self._mlocal,self._nglobal) + blocked_view = blocked[:self._datasize].reshape( + (self._mlocal, self._nglobal), copy=False) for start, end in zip(self._cartstart, self._cartend): - contiguouspart = contiguous[start*self._mlocal:end*self._mlocal] - contiguouspart.shape = (self._mlocal,end-start) + contiguouspart = np.reshape( + contiguous[start*self._mlocal:end*self._mlocal], + (self._mlocal, end-start), copy=False) blocked_view[:,start:end] = contiguouspart def _contiguous_to_blocked(self, blocked, contiguous): @@ -891,11 +888,12 @@ def _contiguous_to_blocked(self, blocked, contiguous): Copies from a contiguous view to a blocked view. Equals roughly a partial transpose, if the block sizes in the cartesian grid are the same. """ - blocked_view = blocked[:self._datasize] - blocked_view.shape = (self._mlocal,self._nglobal) + blocked_view = blocked[:self._datasize].reshape( + (self._mlocal, self._nglobal), copy=False) for start, end in zip(self._cartstart, self._cartend): - contiguouspart = contiguous[start*self._mlocal:end*self._mlocal] - contiguouspart.shape = (self._mlocal,end-start) + contiguouspart = np.reshape( + contiguous[start*self._mlocal:end*self._mlocal], + (self._mlocal, end-start), copy=False) contiguouspart[:] = blocked_view[:,start:end] def solve_pass(self, workmem, tempmem): diff --git a/psydac/mapping/discrete.py b/psydac/mapping/discrete.py index 61f7e0ac9..d9e128506 100644 --- a/psydac/mapping/discrete.py +++ b/psydac/mapping/discrete.py @@ -4,8 +4,6 @@ # for full license details. # #---------------------------------------------------------------------------# from itertools import product -import string -import random import yaml import numpy as np @@ -819,6 +817,10 @@ def __init__(self, *components, name=None): SplineMapping.__init__(self, *components, name=name) + # The weights must be defined on the same space as the coordinates + assert isinstance(weights, FemField) + assert weights.space is self._space + self._weights = NurbsMapping.Weights(self) self._weights_field = weights @@ -843,14 +845,13 @@ def from_control_points_weights(cls, tensor_space, control_points, weights): fields += [FemField(tensor_space)] # Get spline coefficients for each coordinate X_i - # we store w*x where w is the weight and x is the control point + # we store the Cartesian control points x, not the homogeneous w*x starts = tensor_space.coeff_space.starts ends = tensor_space.coeff_space.ends idx_to = tuple(slice(s, e+1) for s,e in zip(starts, ends)) for i, field in enumerate(fields[:-1]): idx_from = (*idx_to, i) -# idw_from = tuple(idx_to) - field.coeffs[idx_to] = control_points[idx_from] #* weights[idw_from] + field.coeffs[idx_to] = control_points[idx_from] # weights idx_from = tuple(idx_to) diff --git a/psydac/mapping/discrete_gallery.py b/psydac/mapping/discrete_gallery.py index 5075d8c85..c9a472f89 100644 --- a/psydac/mapping/discrete_gallery.py +++ b/psydac/mapping/discrete_gallery.py @@ -3,6 +3,8 @@ # LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE # # for full license details. # #---------------------------------------------------------------------------# +from typing import Iterable + import numpy as np from mpi4py import MPI @@ -33,6 +35,11 @@ 'spherical_shell', ) +__all__ = ( + 'get_available_mappings', + 'discrete_mapping', +) + class Collela3D( Mapping ): _expressions = {'x':'2.*(x1 + 0.1*sin(2.*pi*x1)*sin(2.*pi*x2)) - 1.', @@ -40,78 +47,149 @@ class Collela3D( Mapping ): 'z':'2.*x3 - 1.'} #============================================================================== -def discrete_mapping(mapping, ncells, degree, **kwargs): +def get_available_mappings(ldim): + """ + Get a list of `mapping` values accepted as argument to `discrete_mapping`. + + Parameters + ---------- + ldim : int + The number of logical dimensions of the topological domain. + + Returns + ------- + tuple + All the accepted values for the `mapping` parameter. + """ + assert isinstance(ldim, int), f'ldim must be int, got {type(ldim).__name__} instead' + assert ldim > 0, f'ldim must be > 0, got {ldim} instead' + + if ldim == 2: + return ('identity', 'collela', 'circle', 'annulus', 'quarter_annulus', + 'target', 'czarny') + elif ldim == 3: + return ('identity', 'collela', 'spherical_shell') + else: + return () - comm = kwargs.pop('comm', MPI.COMM_WORLD) - return_space = kwargs.pop('return_space', False) +#============================================================================== +def discrete_mapping(mapping, ncells, degree, *, + comm = MPI.COMM_WORLD, + mpi_dims_mask = None, + return_space = False): + """ + Create a SplineMapping by interpolating one of the available analytical mappings. + + Parameters + ---------- + mapping : str + The name of the mapping. See `available_mappings` to get the options. + + ncells : Iterable[int] + The number of cells along each logical dimension. + + degree : Iterable[int] + The spline degree along each logical dimension. + + comm : MPI.Intracomm, optional + The MPI intracommunicator. + + mpi_dims_mask : Iterable[bool], optional + True if the dimension is to be used in the domain decomposition + (=default for each dimension). If mpi_dims_mask[i]=False, the i-th + dimension will not be decomposed. + + return_space : bool, optional + Whether this function should also return the discrete space it creates. + + Returns + ------- + map_discrete : SplineMapping + The spline mapping created. + + space : TensorFemSpace + The space of the components of the spline mapping. + Only returned if `return_space` is True. + """ + # Check types + assert isinstance(mapping, str) + assert isinstance(ncells, Iterable) + assert isinstance(degree, Iterable) + assert isinstance(comm, MPI.Intracomm) or comm is None + assert isinstance(return_space, bool) + + # Check consistency of ncells and degree + assert all(isinstance(n, int) and n >= 1 for n in ncells) + assert all(isinstance(d, int) and d >= 0 for d in degree) + assert len(ncells) == len(degree) mapping = mapping.lower() - dim = len(ncells) - if dim not in [2, 3]: + ldim = len(ncells) + if ldim not in [2, 3]: raise NotImplementedError('Only 2D and 3D mappings are available') # ... - if dim == 2: + if ldim == 2: # Input parameters if mapping == 'identity': - map_symbolic = IdentityMapping('M', dim=dim) + map_symbolic = IdentityMapping('M', dim=ldim) limits = ((0, 1), (0, 1)) periodic = (False, False) elif mapping == 'collela': default_params = dict(k1=1.0, k2=1.0, eps=0.1) - map_symbolic = CollelaMapping2D('M', dim=dim, **default_params) + map_symbolic = CollelaMapping2D('M', dim=ldim, **default_params) limits = ((0, 1), (0, 1)) periodic = (False, False) elif mapping == 'circle': default_params = dict(rmin=0.0, rmax=1.0, c1=0.0, c2=0.0) - map_symbolic = PolarMapping('M', dim=dim, **default_params) + map_symbolic = PolarMapping('M', dim=ldim, **default_params) limits = ((0, 1), (0, 2*np.pi)) periodic = (False, True) elif mapping == 'annulus': default_params = dict(rmin=0.0, rmax=1.0, c1=0.0, c2=0.0) - map_symbolic = PolarMapping('M', dim=dim, **default_params) + map_symbolic = PolarMapping('M', dim=ldim, **default_params) limits = ((1, 4), (0, 2*np.pi)) periodic = (False, True) elif mapping == 'quarter_annulus': default_params = dict(rmin=0.0, rmax=1.0, c1=0.0, c2=0.0) - map_symbolic = PolarMapping('M', dim=dim, **default_params) + map_symbolic = PolarMapping('M', dim=ldim, **default_params) limits = ((1, 4), (0, np.pi/2)) periodic = (False, False) elif mapping == 'target': default_params = dict(c1=0, c2=0, k=0.3, D=0.2) - map_symbolic = TargetMapping('M', dim=dim, **default_params) + map_symbolic = TargetMapping('M', dim=ldim, **default_params) limits = ((0, 1), (0, 2*np.pi)) periodic = (False, True) elif mapping == 'czarny': default_params = dict(c2=0, b=1.4, eps=0.3) - map_symbolic = CzarnyMapping('M', dim=dim, **default_params) + map_symbolic = CzarnyMapping('M', dim=ldim, **default_params) limits = ((0, 1), (0, 2*np.pi)) periodic = (False, True) else: raise ValueError("Required 2D mapping not available") - elif dim == 3: + elif ldim == 3: # Input parameters if mapping == 'identity': - map_symbolic = IdentityMapping('M', dim=dim) + map_symbolic = IdentityMapping('M', dim=ldim) limits = ((0, 1), (0, 1), (0, 1)) periodic = ( False, False, False) elif mapping == 'collela': - map_symbolic = Collela3D('M', dim=dim) + map_symbolic = Collela3D('M', dim=ldim) limits = ((0, 1), (0, 1), (0, 1)) periodic = ( False, False, False) elif mapping == 'spherical_shell': - map_symbolic = SphericalMapping('M', dim=dim) + map_symbolic = SphericalMapping('M', dim=ldim) limits = ((1, 4), (0, np.pi), (0, np.pi/2)) periodic = ( False, False, False) @@ -121,7 +199,8 @@ def discrete_mapping(mapping, ncells, degree, **kwargs): # ... # Create the domain decomposition - domain_decomposition = DomainDecomposition(ncells=ncells, periods=periodic, comm=comm) + domain_decomposition = DomainDecomposition(ncells=ncells, periods=periodic, + comm=comm, mpi_dims_mask=mpi_dims_mask) # Create 1D spline spaces, not distributed spaces_1d = [SplineSpace(grid=np.linspace(*lims, num=nc+1), degree=p, periodic=per) diff --git a/psydac/mapping/tests/test_discrete_mapping.py b/psydac/mapping/tests/test_discrete_mapping.py index 3bb74a641..85319d57d 100644 --- a/psydac/mapping/tests/test_discrete_mapping.py +++ b/psydac/mapping/tests/test_discrete_mapping.py @@ -16,6 +16,7 @@ from psydac.api.discretization import discretize from psydac.core.bsplines import cell_index +from psydac.fem.basic import FemField from psydac.fem.tensor import TensorFemSpace from psydac.fem.splines import SplineSpace from psydac.mapping.discrete import NurbsMapping @@ -309,3 +310,26 @@ def test_nurbs_circle(): J_i = disk.gradient(u=x1, v=x2) assert np.allclose(J_i[:2], J_p, atol=ATOL, rtol=RTOL) + +#============================================================================== +def make_tensor_space_2d(degree: tuple[int, int]) -> TensorFemSpace: + """Create a serial 2D spline space on a uniform grid of 3 x 5 cells.""" + ncells = [3, 5] + domain_decomposition = DomainDecomposition(ncells, [False, False]) + spaces = [SplineSpace(degree=p, grid=np.linspace(0.0, 1.0, n + 1)) + for p, n in zip(degree, ncells)] + return TensorFemSpace(domain_decomposition, *spaces) + +@pytest.mark.parametrize('weights_kind', ['other_space', 'array']) +def test_nurbs_mapping_wrong_weights(weights_kind: str) -> None: + + V = make_tensor_space_2d((2, 2)) + x, y = FemField(V), FemField(V) + + if weights_kind == 'other_space': + weights = FemField(make_tensor_space_2d((1, 3))) + else: + weights = np.ones(V.nbasis) + + with pytest.raises(AssertionError): + NurbsMapping(x, y, weights) diff --git a/pyproject.toml b/pyproject.toml index 9492defb3..f192034ef 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,7 +31,7 @@ dependencies = [ 'termcolor', # Our packages from PyPi - 'sympde == 0.19.3', + 'sympde == 0.20.0', 'pyccel >= 2.2.3', 'gelato == 0.12', @@ -76,3 +76,40 @@ psydac = "psydac.cmd.main:psydac_command" setup = ['--default-library=static'] dist = ['--include-subprojects'] install = ['--only=python-modules'] + +[tool.pylint.main] +py-version = '3.10' +ignore-paths = ['psydac/pyccel/.*', '.*/__e?pyccel__/.*', '.*/__psydac__/.*'] +ignored-modules = ['mpi4py'] # C extension, members not visible to Pylint + +[tool.pylint.basic] +no-docstring-rgx = '^(_|test_|Test)' + +[tool.pylint.'messages control'] +disable = [ + # formatting owned by black / isort + 'line-too-long', + 'trailing-whitespace', + 'bad-indentation', + 'multiple-statements', + 'missing-final-newline', + 'trailing-newlines', + 'wrong-import-order', + # style choices for mathematical code + 'invalid-name', + 'too-many-arguments', + 'too-many-positional-arguments', + 'too-many-locals', + 'too-many-branches', + 'too-many-statements', + 'too-many-instance-attributes', + 'too-few-public-methods', + 'too-many-nested-blocks', + # pytest fixtures shadow outer names by design + 'redefined-outer-name', + 'duplicate-code', + 'fixme', +] + +[tool.isort] +profile = 'black' # sort imports as black formats them diff --git a/subprojects/igakit b/subprojects/igakit index 92ee0975d..2a0c9ee8f 160000 --- a/subprojects/igakit +++ b/subprojects/igakit @@ -1 +1 @@ -Subproject commit 92ee0975d7089e813c1171bd7dbf35a179cf31a4 +Subproject commit 2a0c9ee8f529b1bff86567adf2750bc86b0a9dcf