Merge remote-tracking branch 'upstream/master' into Lorenz_resilience…

…_plots

docs/source/index.rst

@@ -53,6 +53,7 @@ Projects
    projects/compression.rst
    projects/second_order.rst
    projects/monodomain.rst
+   projects/GPU.rst
 
 
 API documentation

docs/source/projects/GPU.rst

@@ -0,0 +1 @@
+.. include:: /../../pySDC/projects/GPU/README.rst

pySDC/helpers/plot_helper.py

@@ -43,11 +43,13 @@ def figsize_by_journal(journal, scale, ratio):  # pragma: no cover
         'JSC_beamer': 426.79135,
         'Springer_Numerical_Algorithms': 338.58778,
         'JSC_thesis': 434.26027,
+        'TUHH_thesis': 426.79135,
     }
     # 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,
     }
     assert (
         journal in textwidths.keys()

pySDC/helpers/pysdc_helper.py

@@ -15,7 +15,7 @@ class FrozenClass(object):
 
     def __setattr__(self, key, value):
         """
-        Function called when setting arttributes
+        Function called when setting attributes
 
         Args:
             key: the attribute
@@ -35,7 +35,7 @@ def __getattr__(self, key):
         if key in self.attrs:
             return None
         else:
-            super().__getattr__(key)
+            super().__getattribute__(key)
 
     @classmethod
     def add_attr(cls, key, raise_error_if_exists=False):

pySDC/implementations/convergence_controller_classes/step_size_limiter.py

@@ -146,11 +146,75 @@ def get_new_step_size(self, controller, S, **kwargs):
                         S,
                     )
                     L.status.dt_new = dt_new
-                elif abs(L.status.dt_new / L.params.dt - 1) < self.params.dt_rel_min_slope:
+                elif abs(L.status.dt_new / L.params.dt - 1) < self.params.dt_rel_min_slope and not S.status.restart:
                     L.status.dt_new = L.params.dt
                     self.log(
                         f"Step size did not change sufficiently to warrant step size change, keeping {L.status.dt_new:.2e}",
                         S,
                     )
 
         return None
+
+
+class StepSizeRounding(ConvergenceController):
+    """
+    Class to round step size when using adaptive step size selection.
+    """
+
+    def setup(self, controller, params, description, **kwargs):
+        """
+        Define parameters here
+
+        Args:
+            controller (pySDC.Controller): The controller
+            params (dict): The params passed for this specific convergence controller
+            description (dict): The description object used to instantiate the controller
+
+        Returns:
+            (dict): The updated params dictionary
+        """
+        defaults = {
+            "control_order": +93,
+            "digits": 1,
+            "fac": 5,
+        }
+        return {**defaults, **super().setup(controller, params, description, **kwargs)}
+
+    @staticmethod
+    def _round_step_size(dt, fac, digits):
+        dt_rounded = None
+        exponent = np.log10(dt) // 1
+
+        dt_norm = dt / 10 ** (exponent - digits)
+        dt_norm_round = (dt_norm // fac) * fac
+        dt_rounded = dt_norm_round * 10 ** (exponent - digits)
+        return dt_rounded
+
+    def get_new_step_size(self, controller, S, **kwargs):
+        """
+        Enforce an upper and lower limit to the step size here.
+        Be aware that this is only tested when a new step size has been determined. That means if you set an initial
+        value for the step size outside of the limits, and you don't do any further step size control, that value will
+        go through.
+        Also, the final step is adjusted such that we reach Tend as best as possible, which might give step sizes below
+        the lower limit set here.
+
+        Args:
+            controller (pySDC.Controller): The controller
+            S (pySDC.Step): The current step
+
+        Returns:
+            None
+        """
+        for L in S.levels:
+            if L.status.dt_new is not None:
+                dt_rounded = self._round_step_size(L.status.dt_new, self.params.fac, self.params.digits)
+
+                if L.status.dt_new != dt_rounded:
+                    self.log(
+                        f"Step size rounded from {L.status.dt_new:.6e} to {dt_rounded:.6e}",
+                        S,
+                    )
+                    L.status.dt_new = dt_rounded
+
+        return None

pySDC/implementations/datatype_classes/cupy_mesh.py

@@ -5,6 +5,11 @@
 except ImportError:
     MPI = None
 
+try:
+    from pySDC.helpers.NCCL_communicator import NCCLComm
+except ImportError:
+    NCCLComm = None
+
 
 class cupy_mesh(cp.ndarray):
     """
@@ -31,7 +36,7 @@ def __new__(cls, init, val=0.0, **kwargs):
             obj[:] = init[:]
         elif (
             isinstance(init, tuple)
-            and (init[1] is None or isinstance(init[1], MPI.Intracomm))
+            and (init[1] is None or isinstance(init[1], MPI.Intracomm) or isinstance(init[1], NCCLComm))
             and isinstance(init[2], cp.dtype)
         ):
             obj = cp.ndarray.__new__(cls, init[0], dtype=init[2], **kwargs)
@@ -62,12 +67,15 @@ def __abs__(self):
             float: absolute maximum of all mesh values
         """
         # take absolute values of the mesh values
-        local_absval = float(cp.amax(cp.ndarray.__abs__(self)))
+        local_absval = cp.max(cp.ndarray.__abs__(self))
 
         if self.comm is not None:
             if self.comm.Get_size() > 1:
-                global_absval = 0.0
-                global_absval = max(self.comm.allreduce(sendobj=local_absval, op=MPI.MAX), global_absval)
+                global_absval = local_absval * 0
+                if isinstance(self.comm, NCCLComm):
+                    self.comm.Allreduce(sendbuf=local_absval, recvbuf=global_absval, op=MPI.MAX)
+                else:
+                    global_absval = self.comm.allreduce(sendobj=float(local_absval), op=MPI.MAX)
             else:
                 global_absval = local_absval
         else:

pySDC/implementations/problem_classes/GrayScott_MPIFFT.py

@@ -1,10 +1,8 @@
 import scipy.sparse as sp
-from mpi4py import MPI
-from mpi4py_fft import PFFT
 
 from pySDC.core.errors import ProblemError
-from pySDC.core.problem import Problem, WorkCounter
-from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh, comp2_mesh
+from pySDC.core.problem import WorkCounter
+from pySDC.implementations.datatype_classes.mesh import comp2_mesh
 from pySDC.implementations.problem_classes.generic_MPIFFT_Laplacian import IMEX_Laplacian_MPIFFT
 
 from mpi4py_fft import newDistArray
@@ -48,6 +46,8 @@ class grayscott_imex_diffusion(IMEX_Laplacian_MPIFFT):
         Denotes the period of the function to be approximated for the Fourier transform.
     comm : COMM_WORLD, optional
         Communicator for ``mpi4py-fft``.
+    num_blobs : int, optional
+        Number of blobs in the initial conditions. Negative values give rectangles.
 
     Attributes
     ----------
@@ -71,18 +71,35 @@ class grayscott_imex_diffusion(IMEX_Laplacian_MPIFFT):
     .. [3] https://www.chebfun.org/examples/pde/GrayScott.html
     """
 
-    def __init__(self, Du=1.0, Dv=0.01, A=0.09, B=0.086, **kwargs):
-        kwargs['L'] = 2.0
-        super().__init__(dtype='d', alpha=1.0, x0=-kwargs['L'] / 2.0, **kwargs)
+    def __init__(
+        self,
+        Du=1.0,
+        Dv=0.01,
+        A=0.09,
+        B=0.086,
+        L=2.0,
+        num_blobs=1,
+        **kwargs,
+    ):
+        super().__init__(dtype='d', alpha=1.0, x0=-L / 2.0, L=L, **kwargs)
 
         # prepare the array with two components
         shape = (2,) + (self.init[0])
         self.iU = 0
         self.iV = 1
         self.ncomp = 2  # needed for transfer class
-        self.init = (shape, self.comm, self.xp.dtype('float'))
 
-        self._makeAttributeAndRegister('Du', 'Dv', 'A', 'B', localVars=locals(), readOnly=True)
+        self.init = (shape, self.comm, self.xp.dtype('complex') if self.spectral else self.xp.dtype('float'))
+
+        self._makeAttributeAndRegister(
+            'Du',
+            'Dv',
+            'A',
+            'B',
+            'num_blobs',
+            localVars=locals(),
+            readOnly=True,
+        )
 
         # prepare "Laplacians"
         self.Ku = -self.Du * self.K2
@@ -168,7 +185,7 @@ def solve_system(self, rhs, factor, u0, t):
 
         return me
 
-    def u_exact(self, t):
+    def u_exact(self, t, seed=10700000):
         r"""
         Routine to compute the exact solution at time :math:`t = 0`, see [3]_.
 
@@ -185,19 +202,135 @@ def u_exact(self, t):
         assert t == 0.0, 'Exact solution only valid as initial condition'
         assert self.ndim == 2, 'The initial conditions are 2D for now..'
 
-        me = self.dtype_u(self.init, val=0.0)
+        xp = self.xp
+
+        _u = xp.zeros_like(self.X[0])
+        _v = xp.zeros_like(self.X[0])
+
+        rng = xp.random.default_rng(seed)
+
+        if self.num_blobs < 0:
+            """
+            Rectangles with stationary background, see arXiv:1501.01990
+            """
+            F, k = self.A, self.B - self.A
+            A = xp.sqrt(F) / (F + k)
+
+            # set stable background state from Equation 2
+            assert 2 * k < xp.sqrt(F) - 2 * F, 'Kill rate is too large to facilitate stable background'
+            _u[...] = (A - xp.sqrt(A**2 - 4)) / (2 * A)
+            _v[...] = xp.sqrt(F) * (A + xp.sqrt(A**2 - 4)) / 2
+
+            for _ in range(-self.num_blobs):
+                x0, y0 = rng.random(size=2) * self.L[0] - self.L[0] / 2
+                lx, ly = rng.random(size=2) * self.L[0] / self.nvars[0] * 30
+
+                mask_x = xp.logical_and(self.X[0] > x0, self.X[0] < x0 + lx)
+                mask_y = xp.logical_and(self.X[1] > y0, self.X[1] < y0 + ly)
+                mask = xp.logical_and(mask_x, mask_y)
+
+                _u[mask] = rng.random()
+                _v[mask] = rng.random()
+
+        elif self.num_blobs > 0:
+            """
+            Blobs as in https://www.chebfun.org/examples/pde/GrayScott.html
+            """
+
+            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
+                        )
+                    )
+                    _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
+                        )
+                    )
+
+            _u += 1
+        else:
+            raise NotImplementedError
 
-        # This assumes that the box is [-L/2, L/2]^2
+        u = self.u_init
         if self.spectral:
-            tmp = 1.0 - self.xp.exp(-80.0 * ((self.X[0] + 0.05) ** 2 + (self.X[1] + 0.02) ** 2))
-            me[0, ...] = self.fft.forward(tmp)
-            tmp = self.xp.exp(-80.0 * ((self.X[0] - 0.05) ** 2 + (self.X[1] - 0.02) ** 2))
-            me[1, ...] = self.fft.forward(tmp)
+            u[0, ...] = self.fft.forward(_u)
+            u[1, ...] = self.fft.forward(_v)
         else:
-            me[0, ...] = 1.0 - self.xp.exp(-80.0 * ((self.X[0] + 0.05) ** 2 + (self.X[1] + 0.02) ** 2))
-            me[1, ...] = self.xp.exp(-80.0 * ((self.X[0] - 0.05) ** 2 + (self.X[1] - 0.02) ** 2))
+            u[0, ...] = _u
+            u[1, ...] = _v
 
-        return me
+        return u
+
+    def get_fig(self, n_comps=2):  # pragma: no cover
+        """
+        Get a figure suitable to plot the solution of this problem
+
+        Args:
+        n_comps (int): Number of components that fit in the solution
+
+        Returns
+        -------
+        self.fig : matplotlib.pyplot.figure.Figure
+        """
+        import matplotlib.pyplot as plt
+        from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+        plt.rcParams['figure.constrained_layout.use'] = True
+
+        if n_comps == 2:
+            self.fig, axs = plt.subplots(1, 2, sharex=True, sharey=True, figsize=((6, 3)))
+            divider = make_axes_locatable(axs[1])
+            self.cax = divider.append_axes('right', size='3%', pad=0.03)
+        else:
+            self.fig, ax = plt.subplots(1, 1, figsize=((6, 5)))
+            divider = make_axes_locatable(ax)
+            self.cax = divider.append_axes('right', size='3%', pad=0.03)
+        return self.fig
+
+    def plot(self, u, t=None, fig=None):  # pragma: no cover
+        r"""
+        Plot the solution. Please supply a figure with the same structure as returned by ``self.get_fig``.
+
+        Parameters
+        ----------
+        u : dtype_u
+            Solution to be plotted
+        t : float
+            Time to display at the top of the figure
+        fig : matplotlib.pyplot.figure.Figure
+            Figure with the correct structure
+
+        Returns
+        -------
+        None
+        """
+        fig = self.get_fig(n_comps=2) if fig is None else fig
+        axs = fig.axes
+
+        vmin = u.min()
+        vmax = u.max()
+        for i, label in zip([self.iU, self.iV], [r'$u$', r'$v$']):
+            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)
+
+        if t is not None:
+            fig.suptitle(f't = {t:.2e}')
+        axs[0].set_xlabel(r'$x$')
+        axs[0].set_ylabel(r'$y$')
+        fig.colorbar(im, self.cax)
 
 
 class grayscott_imex_linear(grayscott_imex_diffusion):

pySDC/implementations/problem_classes/generic_MPIFFT_Laplacian.py

@@ -66,8 +66,6 @@ def setup_GPU(cls):
     def __init__(
         self, nvars=None, spectral=False, L=2 * np.pi, alpha=1.0, comm=MPI.COMM_WORLD, dtype='d', useGPU=False, x0=0.0
     ):
-        """Initialization routine"""
-
         if useGPU:
             self.setup_GPU()
 
@@ -98,7 +96,7 @@ def __init__(
         # invoke super init, passing the communicator and the local dimensions as init
         super().__init__(init=(tmp_u.shape, comm, tmp_u.dtype))
         self._makeAttributeAndRegister(
-            'nvars', 'spectral', 'L', 'alpha', 'comm', 'x0', localVars=locals(), readOnly=True
+            'nvars', 'spectral', 'L', 'alpha', 'comm', 'x0', 'useGPU', localVars=locals(), readOnly=True
         )
 
         # get local mesh