diff --git a/CHANGELOG.md b/CHANGELOG.md index ad4d4aba0..f9dc0e1c3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,8 @@ All notable changes to this project will be documented in this file. - #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] Add `CLAUDE.md` with developer rules ### Fixed diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 000000000..d25c1bb03 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,74 @@ +# 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 everywhere +- meson-python + mesonpy 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) + +## Code standards +- Follow PEP 8 style guide whenever possible +- Docstrings of public functions and classes follow Numpydoc conventions +- pylint + black + isort for linting/typing +- 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" + +## 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`) +- Test file names follow `test_.py` +- Test function names follow `test_` or `test__` +- Parametrize unit tests with `@pytest.parametrize` to minimize code duplication +- Keep run time of tests at a minimum +- Aim for 100% coverage on newly committed code + +## 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/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 + +## Always do +- Before claiming done: + . Run `python -m black --check $(git ls-files "*.py")` + . Run `python -m isort --check $(git ls-files "*.py" | grep -v "FOUND_DUPLICATED_IMPORT.py")` + . Run `pylint --disable=all --enable=unused-import ./psydac --ignore-paths="__pyccel*" --ignore-paths="__epyccel*"` +- 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/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/tests/test_geometry.py b/psydac/cad/tests/test_geometry.py index 18ad343b1..1811fde14 100644 --- a/psydac/cad/tests/test_geometry.py +++ b/psydac/cad/tests/test_geometry.py @@ -365,14 +365,6 @@ 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 #============================================================================== diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index fc3cce3cc..d0e6746b2 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -7,7 +7,7 @@ This module provides iterative solvers and preconditioners. """ -from math import sqrt +from math import inf, sqrt import numpy as np import warnings @@ -239,7 +239,8 @@ def solve_without_pc(self, b, out=None): template = "| {:7d} | {:19.2e} |" print(template.format(1, sqrt(am))) - # Iterate to convergence + # Iterate to convergence (iteration 1 is the initial residual) + m = 1 for m in range(2, maxiter+1): if am < tol_sqr: m -= 1 @@ -339,7 +340,8 @@ def solve_with_pc(self, b, out=None): template = "| {:7d} | {:19.2e} |" print( template.format(1, sqrt(nrmr_sqr))) - # Iterate to convergence + # Iterate to convergence (iteration 1 is the initial residual) + k = 1 for k in range(2, maxiter+1): if nrmr_sqr < tol_sqr: @@ -906,10 +908,10 @@ def solve_with_pc(self, b, out=None): while res_sqr > tol_sqr and niter < maxiter: - # v = A @ pp, vp = PC @ v, alphap = rhop/(vp.rp0) + # v = A @ pp, vp = PC @ v, alphap = rhop/(rp0.vp) A.dot(pp, out=v) pc.dot(v, out=vp) - alphap = rhop / vp.inner(rp0) + alphap = rhop / rp0.inner(vp) # s = r - alphap*v, sp = PC @ s r.copy(out=s) @@ -942,8 +944,8 @@ def solve_with_pc(self, b, out=None): tp *= omegap rp -= tp - # rhop_new = rp.rp0, betap = (alphap*rhop_new)/(omegap*rhop) - rhop_new = rp.inner(rp0) + # rhop_new = rp0.rp, betap = (alphap*rhop_new)/(omegap*rhop) + rhop_new = rp0.inner(rp) betap = (alphap*rhop_new) / (omegap*rhop) rhop = 1*rhop_new @@ -1218,12 +1220,17 @@ def solve(self, b, out=None): ynorm = sqrt(x.inner(x)) rnorm = phibar - if ynorm == 0 or Anorm == 0:test1 = inf - #else:test1 = rnorm / (Anorm*ynorm) # ||r|| / (||A|| ||x||) - else:test1 = rnorm # ||r|| + if ynorm == 0 or Anorm == 0: + test1 = inf + #else: + # test1 = rnorm / (Anorm*ynorm) # ||r|| / (||A|| ||x||) + else: + test1 = rnorm # ||r|| - if Anorm == 0:test2 = inf - else:test2 = root / Anorm # ||Ar|| / (||A|| ||r||) + if Anorm == 0: + test2 = inf + else: + test2 = root / Anorm # ||Ar|| / (||A|| ||r||) # Estimate cond(A). # In this version we look at the diagonals of R in the @@ -1485,6 +1492,15 @@ def solve(self, b, out=None): if conlim > 0:ctol = 1 / conlim normr = beta + # Early exit as in SciPy: the solution of A x = 0 is x = 0, and + # A^H (b - A x) = 0 means that x is already a least-squares solution + if normb == 0: + x *= 0.0 + normr = 0.0 + istop = 1 + elif alpha * beta == 0: + istop = 1 if beta == 0 else 2 + # Reverse the order here from the original matlab code because if verbose: @@ -1495,7 +1511,8 @@ def solve(self, b, out=None): template = "| {:7d} | {:19.2e} |" # Main iteration loop. - for itn in range(1, maxiter + 1): + while istop == 0 and itn < maxiter: + itn += 1 # Perform the next step of the bidiagonalization to obtain the # next beta, u, alpha, v. These satisfy the relations @@ -1596,7 +1613,7 @@ def solve(self, b, out=None): test1 = normr / normb if (normA * normr) != 0:test2 = normar / (normA * normr) - else:test2 = np.infty + else:test2 = inf test3 = 1 / condA t1 = test1 / (1 + normA * normx / normb) rtol = btol + atol * normA * normx / normb @@ -1777,8 +1794,6 @@ def solve(self, b, out=None): # Iterate to convergence for k in range(maxiter): - if am < tol: - break # run Arnoldi self.arnoldi(k, p) @@ -1788,22 +1803,25 @@ def solve(self, b, out=None): # update the residual vector beta.append(- sn[k] * beta[k]) - beta[k] *= cn[k] + beta[k] *= cn[k].conjugate() am = abs(beta[k+1]) if verbose: print( template.format( k+2, am ) ) + if am < tol: + break + if verbose: - print( "+---------+---------------------+") - # calculate result - y = self.solve_triangular(self._H[:k, :k], beta[:k]) # system of upper triangular matrix + print( "+---------+---------------------+") + # calculate result from all k+1 Arnoldi vectors + y = self.solve_triangular(self._H[:k+1, :k+1], beta[:k+1]) # system of upper triangular matrix - for i in range(k): + for i in range(k+1): x.mul_iadd(y[i], self._Q[i]) # Convergence information - self._info = {'niter': k+1, 'success': bool(am < tol), 'res_norm': am} + self._info = {'niter': k+2, 'success': bool(am < tol), 'res_norm': am} if recycle: x.copy(out=self._options["x0"]) @@ -1828,7 +1846,7 @@ def arnoldi(self, k, p): self._A.dot( self._Q[k] , out=p) # Krylov vector for i in range(k + 1): # Modified Gram-Schmidt, keeping Hessenberg matrix - h[i] = p.inner(self._Q[i]) + h[i] = self._Q[i].inner(p) p.mul_iadd(-h[i], self._Q[i]) h[k+1] = sqrt(p.inner(p).real) @@ -1841,23 +1859,24 @@ def arnoldi(self, k, p): def apply_givens_rotation(self, k, sn, cn): # Apply Givens rotation to last column of H + # Rotation [[conj(c), conj(s)], [-s, c]] is unitary also for complex c, s h = self._H[:k+2, k] for i in range(k): h_i_prev = h[i] - h[i] *= cn[i] - h[i] += sn[i] * h[i+1] + h[i] *= cn[i].conjugate() + h[i] += sn[i].conjugate() * h[i+1] h[i+1] *= cn[i] h[i+1] -= sn[i] * h_i_prev - - mod = (h[k]**2 + h[k+1]**2)**0.5 + + mod = (abs(h[k])**2 + abs(h[k+1])**2)**0.5 cn.append( h[k] / mod ) sn.append( h[k+1] / mod ) - h[k] *= cn[k] - h[k] += sn[k] * h[k+1] + h[k] *= cn[k].conjugate() + h[k] += sn[k].conjugate() * h[k+1] h[k+1] = 0. # becomes triangular def dot(self, b, out=None): diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index ee122d876..1eeab1fa4 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -85,9 +85,7 @@ def define_data(n, p, matrix_data, dtype=float): def test_solver_tridiagonal(n, p, dtype, solver, use_jacobi_pc, verbose=False): # Quickly skip tests that are not relevant - if solver == 'BiCGSTAB' and use_jacobi_pc and dtype == complex: - pytest.skip("Preconditioned BiCGSTAB only works for real matrices") - elif solver == 'MINRES' and dtype == complex: + if solver == 'MINRES' and dtype == complex: pytest.skip("MINRES only works for real matrices") # Also skip some problematic tests for now -- see Issue #557 @@ -252,6 +250,74 @@ def test_solver_tridiagonal(n, p, dtype, solver, use_jacobi_pc, verbose=False): assert errh_norm < tol assert (solver == 'CG' and use_jacobi_pc) or errc_norm < tol +#=============================================================================== +# Diagonal operators make the stopping tests divide by a vanishing norm +@pytest.mark.parametrize( + ('solver', 'diagonal', 'expected'), [('MINRES', 0.0, 0.0), ('LSMR', 2.0, 0.5)] +) +def test_solver_diagonal(solver: str, diagonal: float, expected: float) -> None: + + V, A, _ = define_data(6, 1, [0.0, diagonal, 0.0]) + b = V.zeros() + b[:] = 1.0 + + x = inverse(A, solver, tol=1e-10) @ b + + assert np.array_equal(x.toarray(), np.full(6, expected)) + +#=============================================================================== +# GMRES converges in at most n iterations, and must report the true residual +@pytest.mark.parametrize('maxiter', [1, 6, 12]) +@pytest.mark.parametrize('dtype', [float, complex]) +def test_GMRES_solve(maxiter: int, dtype: type) -> None: + + n = 12 + diagonals = [-7-2j, -6-2j, -1-10j] if dtype == complex else [-7, -6, -1] + _, A, xe = define_data(n, 1, diagonals, dtype=dtype) + b = A @ xe + + solver = inverse(A, 'GMRES', tol=1e-10, maxiter=maxiter) + x = solver @ b + info = solver.get_info() + + r = b - A @ x + assert np.isclose(info['res_norm'], np.sqrt(r.inner(r).real), rtol=1e-8, atol=1e-12) + assert info['success'] == (maxiter == n) + +#=============================================================================== +# With maxiter=1, CG only evaluates the residual of the initial guess +@pytest.mark.parametrize('use_jacobi_pc', [False, True]) +def test_ConjugateGradient_solve_maxiter_1(use_jacobi_pc: bool) -> None: + + V, A, xe = define_data_hermitian(6, 1) + pc = A.diagonal(inverse=True) if use_jacobi_pc else None + + solver = inverse(A, 'CG', pc=pc, tol=1e-10, maxiter=1) + x = solver @ (A @ xe) + info = solver.get_info() + + assert np.array_equal(x.toarray(), V.zeros().toarray()) + assert info['niter'] == 1 + assert not info['success'] + +#=============================================================================== +# LSMR returns before iterating if b = 0 (x = 0) or if A^H (b - A x0) = 0 (x = x0) +@pytest.mark.parametrize(('diagonal', 'rhs'), [(2.0, 0.0), (0.0, 1.0)]) +def test_LSMR_solve_early_exit(diagonal: float, rhs: float) -> None: + + V, A, xe = define_data(6, 1, [0.0, diagonal, 0.0]) + b = V.zeros() + b[:] = rhs + + solver = inverse(A, 'LSMR', x0=xe, tol=1e-10) + x = solver @ b + info = solver.get_info() + + expected = xe if rhs else V.zeros() + assert np.array_equal(x.toarray(), expected.toarray()) + assert info['niter'] == 0 + assert info['success'] + #=============================================================================== def test_LST_preconditioner(comm=None): diff --git a/pyproject.toml b/pyproject.toml index ad867dc80..36c290a3e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -76,3 +76,37 @@ 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', +]