Merge branch 'pyccel:devel' into codecov-test

.github/workflows/continuous-integration.yml

@@ -58,6 +58,7 @@ jobs:
       - name: Set default MPI and HDF5 C compilers on Ubuntu
         if: matrix.os == 'ubuntu-latest'
         run: |
+          sudo apt-get update
           sudo apt-get install --reinstall openmpi-bin libhdf5-openmpi-dev
 
       - name: Install non-Python dependencies on macOS
@@ -159,9 +160,9 @@ jobs:
           python -m pip install .
           python -m pip freeze
 
-      - name: Test pyccel optimization flags
+      - name: Test Pyccel optimization flags
         run: |
-            pyccel --verbose ./psydac/linalg/kernels/axpy_kernels.py
+          pytest --pyargs psydac -m pyccel --capture=no
 
       - name: Initialize test directory
         run: |

.github/workflows/documentation.yml

@@ -18,9 +18,9 @@ jobs:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN}}
     steps:
     - name: Checkout
-      uses: actions/checkout@v3
+      uses: actions/checkout@v4
     - name: Set up Python
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: 3.9
     - name: Install non-Python dependencies on Ubuntu
@@ -37,9 +37,9 @@ jobs:
         make -C docs html
         python docs/update_links.py
     - name: Setup Pages
-      uses: actions/configure-pages@v3
+      uses: actions/configure-pages@v5
     - name: Upload artifact
-      uses: actions/upload-pages-artifact@v1
+      uses: actions/upload-pages-artifact@v3
       with:
         path: 'docs/build/html'
 
@@ -53,4 +53,4 @@ jobs:
     steps:
     - name: Deploy to GitHub Pages
       id: deployment
-      uses: actions/deploy-pages@v2
+      uses: actions/deploy-pages@v4

README.md

@@ -186,7 +186,7 @@ python3 <PSYDAC-PATH>/mpi_tester.py --pyargs psydac -m "parallel and petsc"
 
 Many of Psydac's low-level Python functions can be translated to a compiled language using the [Pyccel](https://github.com/pyccel/pyccel) transpiler. Currently, all of those functions are collected in modules which follow the name pattern `[module]_kernels.py`.
 
-The classical installation translates all kernel files to Fortran without user intervention. This does not happen in the case of an editable install, but the command `psydac-accelerate` is made available to the user instead. This command applies Pyccel to all the kernel files in the source directory, and the C language may be selected instead of Fortran (which is the default).
+The classical installation translates all kernel files to Fortran without user intervention. This does not happen in the case of an editable install, but the command `psydac-accelerate` is made available to the user instead. This command applies Pyccel to all the kernel files in the source directory. The default language is currently Fortran, C should also be supported in a near future.
 
 -   **Only in development mode**:
     ```bash

psydac/api/settings.py

@@ -1,5 +1,8 @@
 # coding: utf-8
+import subprocess # nosec B404
 import platform
+import re
+from packaging.version import Version
 
 
 __all__ = ('PSYDAC_DEFAULT_FOLDER', 'PSYDAC_BACKENDS')
@@ -40,18 +43,26 @@
                        'openmp'  : False}
 # ...
 
+# Get gfortran version
+gfortran_version_output = subprocess.check_output(['gfortran', '--version']).decode('utf-8') # nosec B603, B607
+gfortran_version_string = re.search("(\d+\.\d+\.\d+)", gfortran_version_output).group()
+gfortran_version = Version(gfortran_version_string)
+
 # Platform-dependent flags
-if platform.system() == "Darwin" and platform.machine() == 'arm64':
+if platform.system() == "Darwin" and platform.machine() == 'arm64' and gfortran_version >= Version("14"):
+
     # Apple silicon requires architecture-specific flags (see https://github.com/pyccel/psydac/pull/411)
-    import subprocess # nosec B404
+    # which are only available on GCC version >= 14
     cpu_brand = subprocess.check_output(['sysctl','-n','machdep.cpu.brand_string']).decode('utf-8') # nosec B603, B607
-    if "Apple M1" in cpu_brand:
-        PSYDAC_BACKEND_GPYCCEL['flags'] += ' -mcpu=apple-m1'
+    if   "Apple M1" in cpu_brand: PSYDAC_BACKEND_GPYCCEL['flags'] += ' -mcpu=apple-m1'
+    elif "Apple M2" in cpu_brand: PSYDAC_BACKEND_GPYCCEL['flags'] += ' -mcpu=apple-m2'
+    elif "Apple M3" in cpu_brand: PSYDAC_BACKEND_GPYCCEL['flags'] += ' -mcpu=apple-m3'
     else:
         # TODO: Support later Apple CPU models. Perhaps the CPU naming scheme could be easily guessed
         # based on the output of 'sysctl -n machdep.cpu.brand_string', but I wouldn't rely on this
         # guess unless it has been manually verified. Loud errors are better than silent failures!
         raise SystemError(f"Unsupported Apple CPU '{cpu_brand}'.")
+
 else:
     # Default architecture flags
     PSYDAC_BACKEND_GPYCCEL['flags'] += ' -march=native -mtune=native'

psydac/api/tests/test_epyccel_flags.py

@@ -0,0 +1,35 @@
+import pytest
+
+
+@pytest.mark.pyccel
+def test_epyccel_flags():
+
+    from pyccel.epyccel import epyccel
+    from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL as backend
+
+    kwargs = {'language'     : 'fortran',
+              'compiler'     : backend['compiler'],
+              'fflags'       : backend['flags'],
+              'accelerators' : ['openmp'] if backend['openmp'] else [],
+              'verbose'      : True,
+              }
+
+    # Function to be Pyccel-ized
+    def f(x : float):
+        return 3 * x
+
+    # Pyccel magic
+    # ------------
+    # Fortran code is generated and then compiled with the selected compiler.
+    # The compiled function is callable from Python through the C Python API.
+    # The necessary C wrapper functions are also generated by Pyccel.
+    fast_f = epyccel(f, **kwargs)
+
+    # Check output of pyccelized function
+    assert fast_f(3.5) == f(3.5)
+
+
+# Interactive usage
+if __name__ == '__main__':
+    test_epyccel_flags()
+    print('PASSED')

psydac/fem/tensor.py

@@ -456,7 +456,7 @@ def eval_fields(self, grid, *fields, weights=None, npts_per_cell=None, overlap=0
                 npts_per_cell = (npts_per_cell,) * self.ldim
             for i in range(self.ldim):
                 ncells_i = len(self.breaks[i]) - 1
-                grid[i] = np.reshape(grid[i], newshape=(ncells_i, npts_per_cell[i]))
+                grid[i] = np.reshape(grid[i], (ncells_i, npts_per_cell[i]))
             out_fields = self.eval_fields_regular_tensor_grid(grid, *fields, weights=weights, overlap=overlap)
             # return a list
             return [np.ascontiguousarray(out_fields[..., i]) for i in range(len(fields))]

psydac/fem/vector.py

@@ -333,21 +333,6 @@ def ncells(self):
     def spaces( self ):
         return self._spaces
 
-    @property
-    def is_block(self):
-        """Returns True if all components are identical spaces."""
-        # TODO - improve this tests. for the moment, we only check the degree,
-        #      - shall we check the bc too?
-
-        degree = [V.degree for V in self.spaces]
-        if self.pdim == 1:
-            return len(np.unique(degree)) == 1
-        else:
-            ns = np.asarray(degree[0])
-            for ms in degree[1:]:
-                if not( np.allclose(ns, np.asarray(ms)) ): return False
-            return True
-
     # ...
     def get_refined_space(self, ncells):
         return self._refined_space[tuple(ncells)]

psydac/linalg/block.py

@@ -777,13 +777,18 @@ def __sub__(self, M):
     #--------------------------------------
     # New properties/methods
     #--------------------------------------
-    def diagonal(self, *, inverse = False, out = None):
+    def diagonal(self, *, inverse = False, sqrt = False, out = None):
         """Get the coefficients on the main diagonal as another BlockLinearOperator object.
 
         Parameters
         ----------
         inverse : bool
             If True, get the inverse of the diagonal. (Default: False).
+            Can be combined with sqrt to get the inverse square root.
+
+        sqrt : bool
+            If True, get the square root of the diagonal. (Default: False).
+            Can be combined with inverse to get the inverse square root.
 
         out : BlockLinearOperator
             If provided, write the diagonal entries into this matrix. (Default: None).
@@ -815,7 +820,7 @@ def diagonal(self, *, inverse = False, out = None):
         # Store the diagonal (or its inverse) into `out`
         for i, j in self.nonzero_block_indices:
             if i == j:
-                out[i, i] = self[i, i].diagonal(inverse = inverse, out = out[i, i])
+                out[i, i] = self[i, i].diagonal(inverse = inverse, sqrt = sqrt, out = out[i, i])
 
         return out
 

psydac/linalg/kernels/stencil2coo_kernels.py

@@ -19,7 +19,7 @@ def stencil2coo_1d_C(A:'T[:,:]', data:'T[:]', rows:'int64[:]', cols:'int64[:]',
         for j1 in range(ncl1):
             value = A[i1+pp1,j1]
             if abs(value) == 0.0:continue
-            J = ((I//cm1)*dm1+j1-dp1)%nc1
+            J = ((I*dm1//cm1)+j1-dp1)%nc1
             rows[nnz] = I
             cols[nnz] = J
             data[nnz] = value
@@ -38,7 +38,7 @@ def stencil2coo_1d_F(A:'T[:,:]', data:'T[:]', rows:'int64[:]', cols:'int64[:]',
         for j1 in range(ncl1):
             value = A[i1+pp1,j1]
             if abs(value) == 0.0:continue
-            J = ((I//cm1)*dm1+j1-dp1)%nc1
+            J = ((I*dm1//cm1)+j1-dp1)%nc1
             rows[nnz] = I
             cols[nnz] = J
             data[nnz] = value
@@ -66,8 +66,8 @@ def stencil2coo_2d_C(A:'T[:,:,:,:]', data:'T[:]', rows:'int64[:]', cols:'int64[:
                 for j2 in range(ncl2):
                     value = A[i1+pp1,i2+pp2,j1,j2]
                     if abs(value) == 0.0:continue
-                    jj1 = ((ii1//cm1)*dm1+j1-dp1)%nc1
-                    jj2 = ((ii2//cm2)*dm2+j2-dp2)%nc2
+                    jj1 = ((ii1*dm1//cm1)+j1-dp1)%nc1
+                    jj2 = ((ii2*dm2//cm2)+j2-dp2)%nc2
 
                     J   = jj1*nc2 + jj2
 
@@ -97,8 +97,8 @@ def stencil2coo_2d_F(A:'T[:,:,:,:]', data:'T[:]', rows:'int64[:]', cols:'int64[:
                 for j2 in range(ncl2):
                     value = A[i1+pp1,i2+pp2,j1,j2]
                     if abs(value) == 0.0:continue
-                    jj1 = ((ii1//cm1)*dm1+j1-dp1)%nc1
-                    jj2 = ((ii2//cm2)*dm2+j2-dp2)%nc2
+                    jj1 = ((ii1*dm1//cm1)+j1-dp1)%nc1
+                    jj2 = ((ii2*dm2//cm2)+j2-dp2)%nc2
 
                     J   = jj2*nc1 + jj1
 
@@ -132,9 +132,9 @@ def stencil2coo_3d_C(A:'T[:,:,:,:,:,:]', data:'T[:]', rows:'int64[:]', cols:'int
                         for j3 in range(ncl3):
                             value = A[i1+pp1,i2+pp2,i3+pp3,j1,j2,j3]
                             if abs(value) == 0.0:continue
-                            jj1 = ((ii1//cm1)*dm1+j1-dp1)%nc1
-                            jj2 = ((ii2//cm2)*dm2+j2-dp2)%nc2
-                            jj3 = ((ii3//cm3)*dm3+j3-dp3)%nc3
+                            jj1 = ((ii1*dm1//cm1)+j1-dp1)%nc1
+                            jj2 = ((ii2*dm2//cm2)+j2-dp2)%nc2
+                            jj3 = ((ii3*dm3//cm3)+j3-dp3)%nc3
 
                             J   = jj1*nc2*nc3 + jj2*nc3 + jj3
 
@@ -170,9 +170,9 @@ def stencil2coo_3d_F(A:'T[:,:,:,:,:,:]', data:'T[:]', rows:'int64[:]', cols:'int
                         for j3 in range(ncl3):
                             value = A[i1+pp1,i2+pp2,i3+pp3,j1,j2,j3]
                             if abs(value) == 0.0:continue
-                            jj1 = ((ii1//cm1)*dm1+j1-dp1)%nc1
-                            jj2 = ((ii2//cm2)*dm2+j2-dp2)%nc2
-                            jj3 = ((ii3//cm3)*dm3+j3-dp3)%nc3
+                            jj1 = ((ii1*dm1//cm1)+j1-dp1)%nc1
+                            jj2 = ((ii2*dm2//cm2)+j2-dp2)%nc2
+                            jj3 = ((ii3*dm3//cm3)+j3-dp3)%nc3
 
                             J   = jj3*nc1*nc2 + jj2*nc1 + jj1
 

psydac/linalg/stencil.py

@@ -1372,14 +1372,19 @@ def _exchange_assembly_data_serial(self):
                 self._data[idx_to] += self._data[idx_from]
 
     # ...
-    def diagonal(self, *, inverse = False, out = None):
+    def diagonal(self, *, inverse = False, sqrt = False, out = None):
         """
         Get the coefficients on the main diagonal as a StencilDiagonalMatrix object.
 
         Parameters
         ----------
         inverse : bool
             If True, get the inverse of the diagonal. (Default: False).
+            Can be combined with sqrt to get the inverse square root.
+
+        sqrt : bool
+            If True, get the square root of the diagonal. (Default: False).
+            Can be combined with inverse to get the inverse square root.
 
         out : StencilDiagonalMatrix
             If provided, write the diagonal entries into this matrix. (Default: None).
@@ -1418,6 +1423,9 @@ def diagonal(self, *, inverse = False, out = None):
         else:
             data = diag.copy()
 
+        if sqrt:
+            np.sqrt(data, out=data)
+
         # If needed create a new StencilDiagonalMatrix object
         if out is None:
             out = StencilDiagonalMatrix(V, W, data)

psydac/linalg/tests/test_linalg.py

@@ -618,6 +618,12 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False):
     ### -1,T & T,-1 --- -1,T,T --- -1,T,-1 --- T,-1,-1 --- T,-1,T (the combinations I test)
     ###
 
+    # Square root test
+    scaled_matrix = B * np.random.random() # Ensure the diagonal elements != 1
+    diagonal_values = scaled_matrix.diagonal(sqrt=False).toarray()
+    sqrt_diagonal_values = scaled_matrix.diagonal(sqrt=True).toarray()
+    assert np.array_equal(sqrt_diagonal_values, np.sqrt(diagonal_values))
+
     tol = 1e-5
     C = inverse(B, 'cg', tol=tol)
     P = B.diagonal(inverse=True)

psydac/linalg/tests/test_stencil_matrix.py

@@ -433,7 +433,7 @@ def test_stencil_matrix_2d_serial_spurious_entries( dtype, p1, p2, s1, s2, P1, P
 @pytest.mark.parametrize('dtype', [float])
 @pytest.mark.parametrize('n1', [7])
 @pytest.mark.parametrize('p1', [1,2])
-@pytest.mark.parametrize('s1', [1])
+@pytest.mark.parametrize('s1', [1, 2])
 @pytest.mark.parametrize('P1', [True, False])
 def test_stencil_matrix_1d_serial_toarray( dtype, n1, p1, s1, P1):
     # Select non-zero values based on diagonal index
@@ -489,8 +489,8 @@ def test_stencil_matrix_1d_serial_toarray( dtype, n1, p1, s1, P1):
 @pytest.mark.parametrize('n2', [8, 7])
 @pytest.mark.parametrize('p1', [1, 2])
 @pytest.mark.parametrize('p2', [1, 3])
-@pytest.mark.parametrize('s1', [1])
-@pytest.mark.parametrize('s2', [1])
+@pytest.mark.parametrize('s1', [1, 2])
+@pytest.mark.parametrize('s2', [1, 2])
 @pytest.mark.parametrize('P1', [True, False])
 @pytest.mark.parametrize('P2', [True, False])
 def test_stencil_matrix_2d_serial_toarray( dtype, n1, n2, p1, p2, s1, s2, P1, P2):
@@ -2191,7 +2191,7 @@ def test_stencil_matrix_2d_serial_backend_switch(dtype, n1, n2, p1, p2, s1, s2,
 @pytest.mark.parametrize('dtype', [float, complex])
 @pytest.mark.parametrize('n1', [20, 67])
 @pytest.mark.parametrize('p1', [1, 2, 3])
-@pytest.mark.parametrize('sh1', [1])
+@pytest.mark.parametrize('sh1', [1, 2])
 @pytest.mark.parametrize('P1', [True, False])
 @pytest.mark.parallel
 def test_stencil_matrix_1d_parallel_toarray(dtype, n1, p1, sh1, P1):

psydac/mapping/discrete.py

@@ -243,7 +243,7 @@ def jac_mat_grid(self, grid, npts_per_cell=None, overlap=0):
                 npts_per_cell = (npts_per_cell,) * self.ldim
             for i in range(self.ldim):
                 ncells_i = len(self.space.breaks[i]) - 1
-                grid[i] = np.reshape(grid[i], newshape=(ncells_i, npts_per_cell[i]))
+                grid[i] = np.reshape(grid[i], (ncells_i, npts_per_cell[i]))
             jac_mats = self.jac_mat_regular_tensor_grid(grid, overlap=overlap)
             return jac_mats
 
@@ -413,7 +413,7 @@ def inv_jac_mat_grid(self, grid, npts_per_cell=None, overlap=0):
                 npts_per_cell = (npts_per_cell,) * self.ldim
             for i in range(self.ldim):
                 ncells_i = len(self.space.breaks[i]) - 1
-                grid[i] = np.reshape(grid[i], newshape=(ncells_i, npts_per_cell[i]))
+                grid[i] = np.reshape(grid[i], (ncells_i, npts_per_cell[i]))
             inv_jac_mats = self.inv_jac_mat_regular_tensor_grid(grid, overlap=overlap)
             return inv_jac_mats
 
@@ -584,7 +584,7 @@ def jac_det_grid(self, grid, npts_per_cell=None, overlap=0):
                 npts_per_cell = (npts_per_cell,) * self.ldim
             for i in range(self.ldim):
                 ncells_i = len(self.space.breaks[i]) - 1
-                grid[i] = np.reshape(grid[i], newshape=(ncells_i, npts_per_cell[i]))
+                grid[i] = np.reshape(grid[i], (ncells_i, npts_per_cell[i]))
             jac_dets = self.jac_det_regular_tensor_grid(grid, overlap=overlap)
             return jac_dets
 

pytest.ini

@@ -9,6 +9,7 @@ markers =
     serial: single-process test,
     parallel: test to be run using 'mpiexec',
     petsc: test requiring a working PETSc installation with petsc4py Python bindings
+    pyccel: test for checking Pyccel setup on machine
 
 python_files = test_*.py
 python_classes =