From c275e02121aced4c6cef4159078b0769eade7003 Mon Sep 17 00:00:00 2001 From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com> Date: Sat, 16 Nov 2024 13:58:03 +0100 Subject: [PATCH] 3D initial conditions in Gray-Scott --- .../problem_classes/GrayScott_MPIFFT.py | 13 ++-- .../generic_MPIFFT_Laplacian.py | 1 + pySDC/projects/GPU/configs/GS_configs.py | 63 ++++++++++++++----- 3 files changed, 54 insertions(+), 23 deletions(-) diff --git a/pySDC/implementations/problem_classes/GrayScott_MPIFFT.py b/pySDC/implementations/problem_classes/GrayScott_MPIFFT.py index ba14ea762d..796a359b14 100644 --- a/pySDC/implementations/problem_classes/GrayScott_MPIFFT.py +++ b/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) diff --git a/pySDC/implementations/problem_classes/generic_MPIFFT_Laplacian.py b/pySDC/implementations/problem_classes/generic_MPIFFT_Laplacian.py index 5487f7e641..8f4ee38efb 100644 --- a/pySDC/implementations/problem_classes/generic_MPIFFT_Laplacian.py +++ b/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 ---------- diff --git a/pySDC/projects/GPU/configs/GS_configs.py b/pySDC/projects/GPU/configs/GS_configs.py index 6ee1374d5d..f57cdc7ea2 100644 --- a/pySDC/projects/GPU/configs/GS_configs.py +++ b/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