3D initial conditions in Gray-Scott

Author: brownbaerchen

pySDC/implementations/problem_classes/GrayScott_MPIFFT.py

@@ -200,7 +200,6 @@ def u_exact(self, t, seed=10700000):
             Exact solution.
         """
         assert t == 0.0, 'Exact solution only valid as initial condition'
-        assert self.ndim == 2, 'The initial conditions are 2D for now..'
 
         xp = self.xp
 
@@ -222,12 +221,13 @@ def u_exact(self, t, seed=10700000):
             _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
+                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
 
-                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)
+                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]
+                for m in masks[1:]:
+                    mask = xp.logical_and(mask, m)
 
                 _u[mask] = rng.random()
                 _v[mask] = rng.random()
@@ -236,6 +236,7 @@ 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)
 

pySDC/implementations/problem_classes/generic_MPIFFT_Laplacian.py

@@ -11,6 +11,7 @@ class IMEX_Laplacian_MPIFFT(Problem):
     r"""
     Generic base class for IMEX problems using a spectral method to solve the Laplacian implicitly and a possible rest
     explicitly. The FFTs are done with``mpi4py-fft`` [1]_.
+    Works in two and three dimensions.
 
     Parameters
     ----------

pySDC/projects/GPU/configs/GS_configs.py

@@ -27,6 +27,7 @@ class GrayScott(Config):
     num_frames = 200
     sweeper_type = 'IMEX'
     res_per_blob = 2**7
+    ndim = 3
 
     def get_LogToFile(self, ranks=None):
         import numpy as np
@@ -49,16 +50,14 @@ def process_solution(L):
                     't': L.time + L.dt,
                     'u': uend[0].get().view(np.ndarray),
                     'v': uend[1].get().view(np.ndarray),
-                    'X': L.prob.X[0].get().view(np.ndarray),
-                    'Y': L.prob.X[1].get().view(np.ndarray),
+                    'X': L.prob.X.get().view(np.ndarray),
                 }
             else:
                 return {
                     't': L.time + L.dt,
                     'u': uend[0],
                     'v': uend[1],
-                    'X': L.prob.X[0],
-                    'Y': L.prob.X[1],
+                    'X': L.prob.X,
                 }
 
         def logging_condition(L):
@@ -75,7 +74,7 @@ def logging_condition(L):
         LogToFile.logging_condition = logging_condition
         return LogToFile
 
-    def plot(self, P, idx, n_procs_list):  # pragma: no cover
+    def plot(self, P, idx, n_procs_list, projection='xy'):  # pragma: no cover
         import numpy as np
         from matplotlib import ticker as tkr
 
@@ -99,19 +98,49 @@ def plot(self, P, idx, n_procs_list):  # pragma: no cover
             vmax['u'] = max([vmax['u'], buffer[f'u-{rank}']['u'].real.max()])
 
         for rank in range(n_procs_list[2]):
-            im = ax.pcolormesh(
-                buffer[f'u-{rank}']['X'],
-                buffer[f'u-{rank}']['Y'],
-                buffer[f'u-{rank}']['v'].real,
-                vmin=vmin['v'],
-                vmax=vmax['v'],
-                cmap='binary',
-            )
+            if len(buffer[f'u-{rank}']['X']) == 2:
+                ax.set_xlabel('$x$')
+                ax.set_ylabel('$y$')
+                im = ax.pcolormesh(
+                    buffer[f'u-{rank}']['X'][0],
+                    buffer[f'u-{rank}']['X'][1],
+                    buffer[f'u-{rank}']['v'].real,
+                    vmin=vmin['v'],
+                    vmax=vmax['v'],
+                    cmap='binary',
+                )
+            else:
+                v3d = buffer[f'u-{rank}']['v'].real
+
+                if projection == 'xy':
+                    slices = [slice(None), slice(None), v3d.shape[2] // 2]
+                    x = buffer[f'u-{rank}']['X'][0][*slices]
+                    y = buffer[f'u-{rank}']['X'][1][*slices]
+                    ax.set_xlabel('$x$')
+                    ax.set_ylabel('$y$')
+                elif projection == 'xz':
+                    slices = [slice(None), v3d.shape[1] // 2, slice(None)]
+                    x = buffer[f'u-{rank}']['X'][0][*slices]
+                    y = buffer[f'u-{rank}']['X'][2][*slices]
+                    ax.set_xlabel('$x$')
+                    ax.set_ylabel('$z$')
+                elif projection == 'yz':
+                    slices = [v3d.shape[0] // 2, slice(None), slice(None)]
+                    x = buffer[f'u-{rank}']['X'][1][*slices]
+                    y = buffer[f'u-{rank}']['X'][2][*slices]
+                    ax.set_xlabel('$y$')
+                    ax.set_ylabel('$z$')
+
+                im = ax.pcolormesh(
+                    x,
+                    y,
+                    v3d[*slices],
+                    vmin=vmin['v'],
+                    vmax=vmax['v'],
+                    cmap='binary',
+                )
             fig.colorbar(im, cax, format=tkr.FormatStrFormatter('%.1f'))
             ax.set_title(f't={buffer[f"u-{rank}"]["t"]:.2f}')
-            ax.set_xlabel('$x$')
-            ax.set_ylabel('$y$')
-            ax.set_aspect(1.0)
             ax.set_aspect(1.0)
         return fig
 
@@ -130,7 +159,7 @@ def get_description(self, *args, res=-1, **kwargs):
         desc['sweeper_params']['QI'] = 'MIN-SR-S'
         desc['sweeper_params']['QE'] = 'PIC'
 
-        desc['problem_params']['nvars'] = (2**8 if res == -1 else res,) * 2
+        desc['problem_params']['nvars'] = (2**8 if res == -1 else res,) * self.ndim
         desc['problem_params']['Du'] = 0.00002
         desc['problem_params']['Dv'] = 0.00001
         desc['problem_params']['A'] = 0.04