Merge remote-tracking branch 'upstream/master' into HEAD

Author: brownbaerchen

.github/workflows/ci_pipeline.yml

@@ -173,7 +173,7 @@ jobs:
   user_firedrake_tests:
     runs-on: ubuntu-latest
     container:
-      image: firedrakeproject/firedrake-vanilla:2025-01
+      image: firedrakeproject/firedrake-vanilla:latest
       options: --user root
       volumes:
         - ${{ github.workspace }}:/repositories
@@ -205,7 +205,7 @@ jobs:
           firedrake-clean
           cd ./pySDC
           coverage run -m pytest --continue-on-collection-errors -v --durations=0 /repositories/pySDC/pySDC/tests -m firedrake
-        timeout-minutes: 120
+        timeout-minutes: 45
       - name: Make coverage report
         run: |
           . /home/firedrake/firedrake/bin/activate

README.md

@@ -19,8 +19,8 @@ implemented.
 
 -   Variants of SDC: explicit, implicit, IMEX, multi-implicit, Verlet,
     multi-level, diagonal, multi-step
--   Variants of PFASST: virtual parallel or MPI-based parallel,
-    classical of multigrid perspective
+-   Variants of PFASST: virtually parallel or MPI-based parallel,
+    classical or multigrid perspective
 -   8 tutorials: from setting up a first collocation problem to SDC,
     PFASST and advanced topics
 -   Projects: many documented projects with defined and tested outcomes
@@ -43,12 +43,11 @@ The code is hosted on GitHub, see
 will give you a core version of `pySDC` to work with,
 working with the developer version is most often the better choice. We
 thus recommend to checkout the code from GitHub and install the
-dependencies e.g. by using a [conda](https://conda.io/en/latest/)
-environment. For this, `pySDC` ships with environment files
-which can be found in the folder `etc/`. Use these as e.g.
+dependencies e.g. by using [micromamba](https://mamba.readthedocs.io/en/latest/user_guide/micromamba.html). For this, `pySDC` ships with environment files
+which can be found in the folder `etc/` or within the projects. Use these as e.g.
 
 ``` bash
-conda env create -f etc/environment-base.yml
+micromamba create -f etc/environment-base.yml
 ```
 
 If you want to install the developer version using `pip` directly from the GitHub repository, use this:

etc/environment-petsc.yml

@@ -9,7 +9,7 @@ dependencies:
   - matplotlib>=3.0
   - dill>=0.2.6
   - mpich
-  - petsc4py
+  - petsc4py<3.22
   - mpi4py>=3.0.0
   - pip
   - pip:

pySDC/helpers/firedrake_ensemble_communicator.py

@@ -67,6 +67,9 @@ def Isend(self, buf, dest, tag=MPI.ANY_TAG):
             return self.ensemble.ensemble_comm.Isend(buf=buf, dest=dest, tag=tag)
         return self.ensemble.isend(buf, dest, tag=tag)[0]
 
+    def Free(self):
+        del self
+
 
 def get_ensemble(comm, space_size):
     return fd.Ensemble(comm, space_size)

pySDC/implementations/problem_classes/GrayScott_MPIFFT.py

@@ -19,7 +19,7 @@ class grayscott_imex_diffusion(IMEX_Laplacian_MPIFFT):
         \frac{\partial u}{\partial t} = D_u \Delta u - u v^2 + A (1 - u),
 
     .. math::
-        \frac{\partial v}{\partial t} = D_v \Delta v + u v^2 - B u
+        \frac{\partial v}{\partial t} = D_v \Delta v + u v^2 - B v
 
     in :math:`x \in \Omega:=[-L/2, L/2]^N` with :math:`N=2,3`. Spatial discretization is done by using
     Fast Fourier transformation for solving the linear parts provided by ``mpi4py-fft`` [2]_, see also
@@ -222,7 +222,7 @@ def u_exact(self, t, seed=10700000):
 
             for _ in range(-self.num_blobs):
                 x0 = rng.random(size=self.ndim) * self.L[0] - self.L[0] / 2
-                l = rng.random(size=self.ndim) * self.L[0] / self.nvars[0] * 30
+                l = rng.random(size=self.ndim) * self.L[0] / self.nvars[0] * 80
 
                 masks = [xp.logical_and(self.X[i] > x0[i], self.X[i] < x0[i] + l[i]) for i in range(self.ndim)]
                 mask = masks[0]
@@ -236,33 +236,54 @@ def u_exact(self, t, seed=10700000):
             """
             Blobs as in https://www.chebfun.org/examples/pde/GrayScott.html
             """
-            assert self.ndim == 2, 'The initial conditions are 2D for now..'
-
             inc = self.L[0] / (self.num_blobs + 1)
 
             for i in range(1, self.num_blobs + 1):
                 for j in range(1, self.num_blobs + 1):
-                    signs = (-1) ** rng.integers(low=0, high=2, size=2)
-
-                    # This assumes that the box is [-L/2, L/2]^2
-                    _u[...] += -xp.exp(
-                        -80.0
-                        * (
-                            (self.X[0] + self.x0 + inc * i + signs[0] * 0.05) ** 2
-                            + (self.X[1] + self.x0 + inc * j + signs[1] * 0.02) ** 2
+                    signs = (-1) ** rng.integers(low=0, high=2, size=self.ndim)
+
+                    if self.ndim == 2:
+                        # This assumes that the box is [-L/2, L/2]^2
+                        _u[...] += -xp.exp(
+                            -80.0
+                            * (
+                                (self.X[0] + self.x0 + inc * i + signs[0] * 0.05) ** 2
+                                + (self.X[1] + self.x0 + inc * j + signs[1] * 0.02) ** 2
+                            )
+                        )
+                        _v[...] += xp.exp(
+                            -80.0
+                            * (
+                                (self.X[0] + self.x0 + inc * i - signs[0] * 0.05) ** 2
+                                + (self.X[1] + self.x0 + inc * j - signs[1] * 0.02) ** 2
+                            )
+                        )
+                    elif self.ndim == 3:
+                        z_pos = self.x0 + rng.random() * self.L[2]
+                        # This assumes that the box is [-L/2, L/2]^3
+                        _u[...] += -xp.exp(
+                            -80.0
+                            * (
+                                (self.X[0] + self.x0 + inc * i + signs[0] * 0.05) ** 2
+                                + (self.X[1] + self.x0 + inc * j + signs[1] * 0.02) ** 2
+                                + (self.X[2] - z_pos + signs[2] * 0.035) ** 2
+                            )
                         )
-                    )
-                    _v[...] += xp.exp(
-                        -80.0
-                        * (
-                            (self.X[0] + self.x0 + inc * i - signs[0] * 0.05) ** 2
-                            + (self.X[1] + self.x0 + inc * j - signs[1] * 0.02) ** 2
+                        _v[...] += xp.exp(
+                            -80.0
+                            * (
+                                (self.X[0] + self.x0 + inc * i - signs[0] * 0.05) ** 2
+                                + (self.X[1] + self.x0 + inc * j - signs[1] * 0.02) ** 2
+                                + (self.X[2] - z_pos - signs[2] * 0.035) ** 2
+                            )
                         )
-                    )
+                    else:
+                        raise NotImplementedError
 
             _u += 1
         else:
-            raise NotImplementedError
+            _u[...] = rng.random(_u.shape)
+            _v[...] = rng.random(_v.shape)
 
         u = self.u_init
         if self.spectral:

pySDC/implementations/problem_classes/TestEquation_0D.py

@@ -2,7 +2,7 @@
 import scipy.sparse as nsp
 
 from pySDC.core.problem import Problem, WorkCounter
-from pySDC.implementations.datatype_classes.mesh import mesh
+from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh
 
 
 class testequation0d(Problem):
@@ -145,3 +145,118 @@ def u_exact(self, t, u_init=None, t_init=None):
         me = self.dtype_u(self.init)
         me[:] = u_init * self.xp.exp((t - t_init) * self.lambdas)
         return me
+
+
+class test_equation_IMEX(Problem):
+    dtype_f = imex_mesh
+    dtype_u = mesh
+    xp = np
+    xsp = nsp
+
+    def __init__(self, lambdas_implicit=None, lambdas_explicit=None, u0=0.0):
+        """Initialization routine"""
+
+        if lambdas_implicit is None:
+            re = self.xp.linspace(-30, 19, 50)
+            im = self.xp.linspace(-50, 49, 50)
+            lambdas_implicit = self.xp.array(
+                [[complex(re[i], im[j]) for i in range(len(re))] for j in range(len(im))]
+            ).reshape((len(re) * len(im)))
+        if lambdas_explicit is None:
+            re = self.xp.linspace(-30, 19, 50)
+            im = self.xp.linspace(-50, 49, 50)
+            lambdas_implicit = self.xp.array(
+                [[complex(re[i], im[j]) for i in range(len(re))] for j in range(len(im))]
+            ).reshape((len(re) * len(im)))
+        lambdas_implicit = self.xp.asarray(lambdas_implicit)
+        lambdas_explicit = self.xp.asarray(lambdas_explicit)
+
+        assert lambdas_implicit.ndim == 1, f'expect flat list here, got {lambdas_implicit}'
+        assert lambdas_explicit.shape == lambdas_implicit.shape
+        nvars = lambdas_implicit.size
+        assert nvars > 0, 'expect at least one lambda parameter here'
+
+        # invoke super init, passing number of dofs, dtype_u and dtype_f
+        super().__init__(init=(nvars, None, self.xp.dtype('complex128')))
+
+        self.A = self.xsp.diags(lambdas_implicit)
+        self._makeAttributeAndRegister(
+            'nvars', 'lambdas_implicit', 'lambdas_explicit', 'u0', localVars=locals(), readOnly=True
+        )
+        self.work_counters['rhs'] = WorkCounter()
+
+    def eval_f(self, u, t):
+        """
+        Routine to evaluate the right-hand side of the problem.
+
+        Parameters
+        ----------
+        u : dtype_u
+            Current values of the numerical solution.
+        t : float
+            Current time of the numerical solution is computed.
+
+        Returns
+        -------
+        f : dtype_f
+            The right-hand side of the problem.
+        """
+
+        f = self.dtype_f(self.init)
+        f.impl[:] = u * self.lambdas_implicit
+        f.expl[:] = u * self.lambdas_explicit
+        self.work_counters['rhs']()
+        return f
+
+    def solve_system(self, rhs, factor, u0, t):
+        r"""
+        Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
+
+        Parameters
+        ----------
+        rhs : dtype_f
+            Right-hand side for the linear system.
+        factor : float
+            Abbrev. for the local stepsize (or any other factor required).
+        u0 : dtype_u
+            Initial guess for the iterative solver.
+        t : float
+            Current time (e.g. for time-dependent BCs).
+
+        Returns
+        -------
+        me : dtype_u
+            The solution as mesh.
+        """
+        me = self.dtype_u(self.init)
+        L = 1 - factor * self.lambdas_implicit
+        L[L == 0] = 1  # to avoid potential divisions by zeros
+        me[:] = rhs
+        me /= L
+        return me
+
+    def u_exact(self, t, u_init=None, t_init=None):
+        """
+        Routine to compute the exact solution at time t.
+
+        Parameters
+        ----------
+        t : float
+            Time of the exact solution.
+        u_init : pySDC.problem.testequation0d.dtype_u
+            Initial solution.
+        t_init : float
+            The initial time.
+
+        Returns
+        -------
+        me : dtype_u
+            The exact solution.
+        """
+
+        u_init = (self.u0 if u_init is None else u_init) * 1.0
+        t_init = 0.0 if t_init is None else t_init * 1.0
+
+        me = self.dtype_u(self.init)
+        me[:] = u_init * self.xp.exp((t - t_init) * (self.lambdas_implicit + self.lambdas_explicit))
+        return me

pySDC/implementations/problem_classes/generic_MPIFFT_Laplacian.py

@@ -85,6 +85,7 @@ def __init__(
             collapse=True,
             backend=self.fft_backend,
             comm_backend=self.fft_comm_backend,
+            grid=(-1,),
         )
 
         # get test data to figure out type and dimensions

pySDC/projects/GPU/README.rst

@@ -45,18 +45,84 @@ For instance, use
 
 .. code-block:: bash
  
-    srun -n 4 python work_precision.py --config=GS_USkate --procs=1/1/4 --useGPU=True --mode=run
-    mpirun -np 8 python work_precision.py --config=GS_USkate --procs=1/1/4 --useGPU=True --mode=plot
-    python work_precision.py --config=GS_USkate --procs=1/1/4 --useGPU=True --mode=video
+    srun -n 4 python run_experiment.pyy --config=GS_USkate --procs=1/1/4 --useGPU=True --mode=run
+    mpirun -np 8 python run_experiment.py --config=GS_USkate --procs=1/1/4 --useGPU=True --mode=plot
+    python run_experiment.py --config=GS_USkate --procs=1/1/4 --useGPU=True --mode=video
 
 to first run the problem, then make plots and then make a video for Gray-Scott with the U-Skate configuration (see arXiv:1501.01990).
 
 To do a parallel scaling test, you can go to JUWELS Booster and use, for instance,
 
 .. code-block:: bash
-   python analysis_scripts/parallel_scaling.py --mode=run --scaling=strong --space_time=True --XPU=GPU --problem=GS
-   srun python analysis_scripts/parallel_scaling.py --mode=plot --scaling=strong --space_time=True --XPU=GPU --problem=GS
+
+   python analysis_scripts/parallel_scaling.py --mode=run --space_time=True --XPU=GPU --problem=GS3D
+   python analysis_scripts/parallel_scaling.py --mode=plot --space_time=True --XPU=GPU --problem=GS3D
 
 This will generate jobscripts and submit the jobs. Notice that you have to wait for the jobs to complete before you can plot them.
 
 To learn more about the options for the scripts, run them with `--help`.
+
+Reproducing plots in Thomas Baumann's thesis
+--------------------------------------------
+Keep in mind that the results of the experiments are specific to the hardware that was used in the experiments.
+To record the data for space-time parallel scaling experiments with Gray-Scott and RBC, run the following commands on the specified machines within the directory that contains this README.
+
+.. code-block:: bash
+
+    # run on JUWELS
+    python analysis_scripts/parallel_scaling.py --mode=run --problem=GS3D --XPU=CPU --space_time=False
+    python analysis_scripts/parallel_scaling.py --mode=run --problem=GS3D --XPU=CPU --space_time=True
+
+    # run on JUWELS booster
+    python analysis_scripts/parallel_scaling.py --mode=run --problem=GS3D --XPU=GPU --space_time=False
+    python analysis_scripts/parallel_scaling.py --mode=run --problem=GS3D --XPU=GPU --space_time=True
+
+    # run on JURECA DC
+    python analysis_scripts/parallel_scaling.py --mode=run --problem=RBC --XPU=CPU --space_time=False
+    python analysis_scripts/parallel_scaling.py --mode=run --problem=RBC --XPU=CPU --space_time=True
+
+    # run on JUWELS booster
+    python analysis_scripts/parallel_scaling.py --mode=run --problem=RBC --XPU=GPU --space_time=False
+    python analysis_scripts/parallel_scaling.py --mode=run --problem=RBC --XPU=GPU --space_time=True
+
+These commands will submit a bunch of jobscripts with the individual runs.
+Keep in mind that these are specific to a compute project and some paths are account-specific.
+Most likely, you will have to change options at the top of the file `./etc/generate_jobscript.py` before you can run anything.
+Also, notice that you may not be allowed to request all resources needed for the largest Gray-Scott GPU run during normal operation of JUWELS booster.
+
+After all jobs have run to completion, you have recorded all scaling data and may plot the results with the following command:
+
+.. code-block:: bash
+
+    python paper_plots.py --target=thesis
+
+In order to run the production runs, modify the `path` class attribute of `LargeSim` in `analysis_scripts/large_simulations.py`.
+Then use the following commands on the specified machines:
+
+.. code-block:: bash
+
+    # run on JUWELS booster
+    python analysis_scripts/large_simulations.py --mode=run --problem=GS --XPU=GPU
+
+    # run on JURECA DC
+    python analysis_scripts/large_simulations.py --mode=run --problem=RBC --XPU=CPU
+
+Plotting the results of the Gray-Scott simulation requires a lot of memory and will take very long.
+Modify the paths in `analysis_scripts/plot_large_simulations.py` and then run:
+
+.. code-block:: bash
+
+    python analysis_scripts/3d_plot_GS_large.py --base_path=<path>
+    python analysis_scripts/plot_large_simulations.py --problem=GS
+
+Plotting the results of the Rayleigh-Benard production run is more easy.
+After modifying the paths as earlier, run the following commands:
+
+.. code-block:: bash
+
+    python analysis_scripts/large_simulations.py --mode=plot --problem=RBC --XPU=CPU
+    python analysis_scripts/large_simulations.py --mode=video --problem=RBC --XPU=CPU
+    python analysis_scripts/plot_large_simulations.py --problem=RBC
+    
+Run scripts with `--help` to learn more about parameters.
+Keep in mind that not all features are supported with all problems.