Merge remote-tracking branch 'origin/master' into neuralpint

Author: tlunet

.github/workflows/ci_pipeline.yml

@@ -39,7 +39,7 @@ jobs:
       fail-fast: false
       matrix:
         env: ['base', 'fenics', 'mpi4py', 'petsc', 'pytorch']
-        python: ['3.9', '3.10', '3.11', '3.12', '3.13']
+        python: ['3.10', '3.11', '3.12', '3.13']
     defaults:
       run:
         shell: bash -l {0}

README.md

@@ -5,6 +5,8 @@
 [![fair-software.eu](https://img.shields.io/badge/fair--software.eu-%E2%97%8F%20%20%E2%97%8F%20%20%E2%97%8F%20%20%E2%97%8F%20%20%E2%97%8F-green)](https://fair-software.eu)
 [![SQAaaS badge shields.io](https://img.shields.io/badge/sqaaas%20software-silver-lightgrey)](https://api.eu.badgr.io/public/assertions/aS8J0NDTTjCyYP6iVufviQ "SQAaaS silver badge achieved")
 [![PyPI - Downloads](https://img.shields.io/pypi/dm/pySDC?logo=pypi)](https://pypistats.org/packages/pysdc)
+[![HiRSE Code Promo Badge](https://img.shields.io/badge/Promo-8db427?label=HiRSE&labelColor=005aa0&link=https%3A%2F%2Fgo.fzj.de%2FCodePromo)](https://go.fzj.de/CodePromo)
+[![SWH](https://archive.softwareheritage.org/badge/origin/https://github.com/Parallel-in-Time/pySDC/)](https://archive.softwareheritage.org/browse/origin/?origin_url=https://github.com/Parallel-in-Time/pySDC)
 
 # Welcome to pySDC!
 
@@ -33,7 +35,7 @@ implemented.
 -   Continuous integration via [GitHub
     Actions](https://github.com/Parallel-in-Time/pySDC/actions) and
     [Gitlab CI](https://gitlab.hzdr.de/r.speck/pysdc/-/pipelines) (through the [GitHub2Gitlab Action](https://github.com/jakob-fritz/github2lab_action))
--   Fully compatible with Python 3.9 - 3.12, runs at least on Ubuntu
+-   Fully compatible with Python 3.10 - 3.13, runs at least on Ubuntu
 
 ## Getting started
 

pySDC/helpers/blocks.py

@@ -102,7 +102,7 @@ def ranks(self):
     @property
     def localBounds(self):
         iLocList, nLocList = [], []
-        for rank, nPoints, nBlocks in zip(self.ranks, self.gridSizes, self.nBlocks):
+        for rank, nPoints, nBlocks in zip(self.ranks, self.gridSizes, self.nBlocks, strict=True):
             n0 = nPoints // nBlocks
             nRest = nPoints - nBlocks * n0
             nLoc = n0 + 1 * (rank < nRest)

pySDC/helpers/fieldsIO.py

@@ -159,7 +159,8 @@ def fromFile(cls, fileName):
         fieldsIO : :class:`FieldsIO`
             The specialized `FieldsIO` adapted to the file.
         """
-        assert os.path.isfile(fileName), f"not a file ({fileName})"
+        if not os.path.isfile(fileName):
+            raise FileNotFoundError(f"not a file ({fileName})")
         with open(fileName, "rb") as f:
             STRUCT, DTYPE = np.fromfile(f, dtype=H_DTYPE, count=2)
             fieldsIO: FieldsIO = cls.STRUCTS[STRUCT](DTYPES[DTYPE], fileName)
@@ -719,7 +720,7 @@ def writeFields_MPI(fileName, dtypeIdx, algo, nSteps, nVar, gridSizes):
 
     iLoc, nLoc = blocks.localBounds
     Rectilinear.setupMPI(comm, iLoc, nLoc)
-    s = [slice(i, i + n) for i, n in zip(iLoc, nLoc)]
+    s = [slice(i, i + n) for i, n in zip(iLoc, nLoc, strict=True)]
     u0 = u0[(slice(None), *s)]
 
     f1 = Rectilinear(DTYPES[dtypeIdx], fileName)

pySDC/helpers/plot_helper.py

@@ -45,13 +45,15 @@ def figsize_by_journal(journal, scale, ratio):  # pragma: no cover
         'Springer_proceedings': 347.12354,
         'JSC_thesis': 434.26027,
         'TUHH_thesis': 426.79135,
+        'Nature_CS': 372.0,
     }
     # store text height in points here, get this from LaTeX using \the\textheight
     textheights = {
         'JSC_beamer': 214.43411,
         'JSC_thesis': 635.5,
         'TUHH_thesis': 631.65118,
         'Springer_proceedings': 549.13828,
+        'Nature_CS': 552.69478,
     }
     assert (
         journal in textwidths.keys()

pySDC/helpers/spectral_helper.py

@@ -97,6 +97,7 @@ def __init__(self, N, x0=None, x1=None, useGPU=False, useFFTW=False):
 
         if useGPU:
             self.setup_GPU()
+            self.logger.debug('Set up for GPU')
         else:
             self.setup_CPU(useFFTW=useFFTW)
 
@@ -161,6 +162,10 @@ def get_differentiation_matrix(self):
     def get_integration_matrix(self):
         raise NotImplementedError()
 
+    def get_integration_weights(self):
+        """Weights for integration across entire domain"""
+        raise NotImplementedError()
+
     def get_wavenumbers(self):
         """
         Get the grid in spectral space
@@ -379,6 +384,16 @@ def get_integration_matrix(self, lbnd=0):
             raise NotImplementedError(f'This function allows to integrate only from x=0, you attempted from x={lbnd}.')
         return S
 
+    def get_integration_weights(self):
+        """Weights for integration across entire domain"""
+        n = self.xp.arange(self.N, dtype=float)
+
+        weights = (-1) ** n + 1
+        weights[2:] /= 1 - (n**2)[2:]
+
+        weights /= 2 / self.L
+        return weights
+
     def get_differentiation_matrix(self, p=1):
         '''
         Keep in mind that the T2T differentiation matrix is dense.
@@ -808,6 +823,12 @@ def get_integration_matrix(self, p=1):
         k[0] = 1j * self.L
         return self.linalg.matrix_power(self.sparse_lib.diags(1 / (1j * k)), p)
 
+    def get_integration_weights(self):
+        """Weights for integration across entire domain"""
+        weights = self.xp.zeros(self.N)
+        weights[0] = self.L / self.N
+        return weights
+
     def get_plan(self, u, forward, *args, **kwargs):
         if self.fft_lib.__name__ == 'mpi4py_fft.fftw':
             if 'axes' in kwargs.keys():
@@ -1006,7 +1027,6 @@ def __init__(self, comm=None, useGPU=False, debug=False):
         self.BCs = None
 
         self.fft_cache = {}
-        self.fft_dealias_shape_cache = {}
 
         self.logger = logging.getLogger(name='Spectral Discretization')
         if debug:
@@ -1293,11 +1313,12 @@ def setup_BCs(self):
 
         diags = self.xp.ones(self.BCs.shape[0])
         diags[self.BC_zero_index] = 0
-        self.BC_line_zero_matrix = sp.diags(diags)
+        self.BC_line_zero_matrix = sp.diags(diags).tocsc()
 
         # prepare BCs in spectral space to easily add to the RHS
         rhs_BCs = self.put_BCs_in_rhs(self.u_init)
-        self.rhs_BCs_hat = self.transform(rhs_BCs)
+        self.rhs_BCs_hat = self.transform(rhs_BCs).view(self.xp.ndarray)
+        del self.BC_rhs_mask
 
     def check_BCs(self, u):
         """
@@ -1350,7 +1371,7 @@ def put_BCs_in_rhs_hat(self, rhs_hat):
             Generate a mask where we need to set values in the rhs in spectral space to zero, such that can replace them
             by the boundary conditions. The mask is then cached.
             """
-            self._rhs_hat_zero_mask = self.newDistArray().astype(bool)
+            self._rhs_hat_zero_mask = self.newDistArray(forward_output=True).astype(bool).view(self.xp.ndarray)
 
             for axis in range(self.ndim):
                 for bc in self.full_BCs:
@@ -1498,7 +1519,7 @@ def get_indices(self, forward_output=True):
     def get_pfft(self, axes=None, padding=None, grid=None):
         if self.ndim == 1 or self.comm is None:
             return None
-        from mpi4py_fft import PFFT
+        from mpi4py_fft import PFFT, newDistArray
 
         axes = tuple(i for i in range(self.ndim)) if axes is None else axes
         padding = list(padding if padding else [1.0 for _ in range(self.ndim)])
@@ -1530,6 +1551,10 @@ def no_transform(u, *args, **kwargs):
             transforms=transforms,
             grid=grid,
         )
+
+        # do a transform to do the planning
+        _u = newDistArray(pfft, forward_output=False)
+        pfft.backward(pfft.forward(_u))
         return pfft
 
     def get_fft(self, axes=None, direction='object', padding=None, shape=None):
@@ -1885,6 +1910,7 @@ def get_differentiation_matrix(self, axes, **kwargs):
             _D = self.axes[axis].get_differentiation_matrix(**kwargs)
             D = D @ self.expand_matrix_ND(_D, axis)
 
+        self.logger.debug(f'Set up differentiation matrix along axes {axes} with kwargs {kwargs}')
         return D
 
     def get_integration_matrix(self, axes):
@@ -1948,4 +1974,5 @@ def get_basis_change_matrix(self, axes=None, **kwargs):
             _C = self.axes[axis].get_basis_change_matrix(**kwargs)
             C = C @ self.expand_matrix_ND(_C, axis)
 
+        self.logger.debug(f'Set up basis change matrix along axes {axes} with kwargs {kwargs}')
         return C

pySDC/helpers/transfer_helper.py

@@ -27,7 +27,7 @@ def next_neighbors_periodic(p, ps, k):
 
     # zip it
     value_index = []
-    for d, i in zip(distance_to_p, range(distance_to_p.size)):
+    for d, i in zip(distance_to_p, range(distance_to_p.size), strict=True):
         value_index.append((d, i))
     # sort by distance
     value_index_sorted = sorted(value_index, key=lambda s: s[0])
@@ -53,7 +53,7 @@ def next_neighbors(p, ps, k):
     distance_to_p = np.abs(ps - p)
     # zip it
     value_index = []
-    for d, i in zip(distance_to_p, range(distance_to_p.size)):
+    for d, i in zip(distance_to_p, range(distance_to_p.size), strict=True):
         value_index.append((d, i))
     # sort by distance
     value_index_sorted = sorted(value_index, key=lambda s: s[0])
@@ -80,7 +80,7 @@ def continue_periodic_array(arr, nn):
     else:
         cont_arr = [arr[nn[0]]]
         shift = 0.0
-        for n, d in zip(nn[1:], d_nn):
+        for n, d in zip(nn[1:], d_nn, strict=True):
             if d != 1:
                 shift = -1
             cont_arr.append(arr[n] + shift)
@@ -107,7 +107,7 @@ def restriction_matrix_1d(fine_grid, coarse_grid, k=2, periodic=False, pad=1):
 
     if periodic:
         M = np.zeros((coarse_grid.size, fine_grid.size))
-        for i, p in zip(range(n_g), coarse_grid):
+        for i, p in zip(range(n_g), coarse_grid, strict=True):
             nn = next_neighbors_periodic(p, fine_grid, k)
             circulating_one = np.asarray([1.0] + [0.0] * (k - 1))
             cont_arr = continue_periodic_array(fine_grid, nn)
@@ -120,7 +120,7 @@ def restriction_matrix_1d(fine_grid, coarse_grid, k=2, periodic=False, pad=1):
                 M[i, nn] = np.asarray(list(map(lambda x: x(p), bary_pol)))
     else:
         M = np.zeros((coarse_grid.size, fine_grid.size + 2 * pad))
-        for i, p in zip(range(n_g), coarse_grid):
+        for i, p in zip(range(n_g), coarse_grid, strict=True):
             padded_f_grid = border_padding(fine_grid, pad, pad)
             nn = next_neighbors(p, padded_f_grid, k)
             # construct the lagrange polynomials for the k neighbors
@@ -158,7 +158,7 @@ def interpolation_matrix_1d(fine_grid, coarse_grid, k=2, periodic=False, pad=1,
         M = np.zeros((fine_grid.size, coarse_grid.size))
 
         if equidist_nested:
-            for i, p in zip(range(n_f), fine_grid):
+            for i, p in zip(range(n_f), fine_grid, strict=True):
                 if i % 2 == 0:
                     M[i, int(i / 2)] = 1.0
                 else:
@@ -189,7 +189,7 @@ def interpolation_matrix_1d(fine_grid, coarse_grid, k=2, periodic=False, pad=1,
                         M[i, nn] = np.asarray(list(map(lambda x: x(p), bary_pol)))
 
         else:
-            for i, p in zip(range(n_f), fine_grid):
+            for i, p in zip(range(n_f), fine_grid, strict=True):
                 nn = next_neighbors_periodic(p, coarse_grid, k)
                 circulating_one = np.asarray([1.0] + [0.0] * (k - 1))
                 cont_arr = continue_periodic_array(coarse_grid, nn)
@@ -208,7 +208,7 @@ def interpolation_matrix_1d(fine_grid, coarse_grid, k=2, periodic=False, pad=1,
         padded_c_grid = border_padding(coarse_grid, pad, pad)
 
         if equidist_nested:
-            for i, p in zip(range(n_f), fine_grid):
+            for i, p in zip(range(n_f), fine_grid, strict=True):
                 if i % 2 != 0:
                     M[i, int((i - 1) / 2) + 1] = 1.0
                 else:
@@ -231,7 +231,7 @@ def interpolation_matrix_1d(fine_grid, coarse_grid, k=2, periodic=False, pad=1,
                         M[i, nn] = np.asarray(list(map(lambda x: x(p), bary_pol)))
 
         else:
-            for i, p in zip(range(n_f), fine_grid):
+            for i, p in zip(range(n_f), fine_grid, strict=True):
                 nn = next_neighbors(p, padded_c_grid, k)
                 # construct the lagrange polynomials for the k neighbors
                 circulating_one = np.asarray([1.0] + [0.0] * (k - 1))

pySDC/helpers/vtkIO.py

@@ -52,7 +52,7 @@ def field(u):
         return numpy_support.numpy_to_vtk(num_array=u.ravel(order='F'), deep=True, array_type=vtk.VTK_FLOAT)
 
     pointData = vtr.GetPointData()
-    for name, u in zip(varNames, data):
+    for name, u in zip(varNames, data, strict=True):
         uVTK = field(u)
         uVTK.SetName(name)
         pointData.AddArray(uVTK)

pySDC/implementations/problem_classes/AllenCahn_Temp_MPIFFT.py

@@ -134,7 +134,7 @@ def __init__(
         N = self.fft.global_shape()
         k = [np.fft.fftfreq(n, 1.0 / n).astype(int) for n in N[:-1]]
         k.append(np.fft.rfftfreq(N[-1], 1.0 / N[-1]).astype(int))
-        K = [ki[si] for ki, si in zip(k, s)]
+        K = [ki[si] for ki, si in zip(k, s, strict=True)]
         Ks = list(np.meshgrid(*K, indexing='ij', sparse=True))
         Lp = 2 * np.pi / L
         for i in range(ndim):

pySDC/implementations/problem_classes/Brusselator.py

@@ -178,7 +178,7 @@ def plot(self, u, t=None, fig=None):  # pragma: no cover
 
         vmin = u.min()
         vmax = u.max()
-        for i, label in zip([self.iU, self.iV], [r'$u$', r'$v$']):
+        for i, label in zip([self.iU, self.iV], [r'$u$', r'$v$'], strict=True):
             im = axs[i].pcolormesh(self.X[0], self.X[1], u[i], vmin=vmin, vmax=vmax)
             axs[i].set_aspect(1)
             axs[i].set_title(label)