improved fft init

Author: pancetta

pySDC/implementations/problem_classes/AllenCahn_2D_FFT.py

@@ -58,17 +58,14 @@ def __init__(self, problem_params, dtype_u=mesh, dtype_f=rhs_imex_mesh):
 
         kx = np.zeros(self.init[0])
         ky = np.zeros(self.init[1] // 2 + 1)
-        for i in range(0, int(self.init[0] / 2) + 1):
-            kx[i] = 2 * np.pi / self.params.L * i
-        for i in range(0, int(self.init[1] // 2) + 1):
-            ky[i] = 2 * np.pi / self.params.L * i
-        for i in range(int(self.init[0] / 2) + 1, self.init[0]):
-            kx[i] = 2 * np.pi / self.params.L * (-self.init[0] + i)
-
-        self.lap = np.zeros((self.init[0], self.init[1] // 2 + 1))
-        for i in range(self.init[0]):
-            for j in range(self.init[1] // 2 + 1):
-                self.lap[i, j] = -kx[i] ** 2 - ky[j] ** 2
+
+        kx[:int(self.init[0] / 2) + 1] = 2 * np.pi / self.params.L * np.arange(0, int(self.init[0] / 2) + 1)
+        kx[int(self.init[0] / 2) + 1:] = 2 * np.pi / self.params.L * \
+            np.arange(int(self.init[0] / 2) + 1 - self.init[0], 0)
+        ky[:] = 2 * np.pi / self.params.L * np.arange(0, self.init[1] // 2 + 1)
+
+        xv, yv = np.meshgrid(kx, ky, indexing='ij')
+        self.lap = -xv ** 2 - yv ** 2
 
         rfft_in = pyfftw.empty_aligned(self.init, dtype='float64')
         fft_out = pyfftw.empty_aligned((self.init[0], self.init[1] // 2 + 1), dtype='complex128')
@@ -133,10 +130,9 @@ def u_exact(self, t):
         assert t == 0, 'ERROR: u_exact only valid for t=0'
         me = self.dtype_u(self.init, val=0.0)
         if self.params.init_type == 'circle':
-            for i in range(self.params.nvars[0]):
-                for j in range(self.params.nvars[1]):
-                    r2 = self.xvalues[i] ** 2 + self.xvalues[j] ** 2
-                    me.values[i, j] = np.tanh((self.params.radius - np.sqrt(r2)) / (np.sqrt(2) * self.params.eps))
+            xv, yv = np.meshgrid(self.xvalues, self.xvalues, indexing='ij')
+            me.values[:, :] = np.tanh((self.params.radius - np.sqrt(xv ** 2 + yv ** 2)) /
+                                      (np.sqrt(2) * self.params.eps))
         elif self.params.init_type == 'checkerboard':
             xv, yv = np.meshgrid(self.xvalues, self.xvalues)
             me.values[:, :] = np.sin(2.0 * np.pi * xv) * np.sin(2.0 * np.pi * yv)
@@ -172,19 +168,7 @@ def __init__(self, problem_params, dtype_u=mesh, dtype_f=rhs_imex_mesh):
         """
         super(allencahn2d_imex_stab, self).__init__(problem_params=problem_params, dtype_u=dtype_u, dtype_f=dtype_f)
 
-        kx = np.zeros(self.init[0])
-        ky = np.zeros(self.init[1] // 2 + 1)
-        for i in range(0, int(self.init[0] / 2) + 1):
-            kx[i] = 2 * np.pi / self.params.L * i
-        for i in range(0, int(self.init[1] // 2) + 1):
-            ky[i] = 2 * np.pi / self.params.L * i
-        for i in range(int(self.init[0] / 2) + 1, self.init[0]):
-            kx[i] = 2 * np.pi / self.params.L * (-self.init[0] + i)
-
-        self.lap = np.zeros((self.init[0], self.init[1] // 2 + 1))
-        for i in range(self.init[0]):
-            for j in range(self.init[1] // 2 + 1):
-                self.lap[i, j] = -kx[i] ** 2 - ky[j] ** 2 - 2.0 / self.params.eps ** 2
+        self.lap -= 2.0 / self.params.eps ** 2
 
     def eval_f(self, u, t):
         """

pySDC/implementations/problem_classes/AllenCahn_2D_parFFT.py

@@ -53,31 +53,32 @@ def __init__(self, problem_params, dtype_u=mesh, dtype_f=rhs_imex_mesh):
         if problem_params['nvars'][0] % 2 != 0:
             raise ProblemError('the setup requires nvars = 2^p per dimension')
 
-        self.FFT = R2C(np.asarray(problem_params['nvars']), np.asarray((1, 1)), problem_params['comm'], "double")
+        self.FFT = R2C(np.asarray(problem_params['nvars']),
+                       np.asarray((problem_params['L'], problem_params['L'])),
+                       problem_params['comm'], "double")
 
         # invoke super init, passing number of dofs, dtype_u and dtype_f
         super(allencahn2d_imex, self).__init__(init=self.FFT.real_shape(), dtype_u=dtype_u, dtype_f=dtype_f,
                                                params=problem_params)
 
+        self.rank = self.params.comm.Get_rank()
+
         self.dx = self.params.L / self.params.nvars[0]  # could be useful for hooks, too.
-        self.xvalues = np.array([i * self.dx - self.params.L / 2.0 for i in range(self.init[0])])
+        self.xvalues = np.array([i * self.dx - self.params.L / 2.0
+                                 for i in range(self.rank * self.init[0], (self.rank + 1) * self.init[0])])
         self.yvalues = np.array([i * self.dx - self.params.L / 2.0 for i in range(self.init[1])])
 
         cinit = self.FFT.complex_shape()
 
         kx = np.zeros(cinit[0])
         ky = np.zeros(cinit[1])
-        for i in range(0, int(cinit[0] / 2) + 1):
-            kx[i] = 2 * np.pi / self.params.L * i
-        for i in range(0, cinit[1]):
-            ky[i] = 2 * np.pi / self.params.L * i
-        for i in range(int(cinit[0] / 2) + 1, cinit[0]):
-            kx[i] = 2 * np.pi / self.params.L * (-cinit[0] + i)
 
-        self.lap = np.zeros((cinit[0], cinit[1]))
-        for i in range(cinit[0]):
-            for j in range(cinit[1]):
-                self.lap[i, j] = -kx[i] ** 2 - ky[j] ** 2
+        kx[:int(cinit[0] / 2) + 1] = 2 * np.pi / self.params.L * np.arange(0, int(cinit[0] / 2) + 1)
+        kx[int(cinit[0] / 2) + 1:] = 2 * np.pi / self.params.L * np.arange(int(cinit[0] / 2) + 1 - cinit[0], 0)
+        ky[:] = 2 * np.pi / self.params.L * np.arange(self.rank * cinit[1], (self.rank + 1) * cinit[1])
+
+        xv, yv = np.meshgrid(kx, ky, indexing='ij')
+        self.lap = -xv ** 2 - yv ** 2
 
         print(self.init, cinit)
 
@@ -143,12 +144,11 @@ def u_exact(self, t):
         assert t == 0, 'ERROR: u_exact only valid for t=0'
         me = self.dtype_u(self.init, val=0.0)
         if self.params.init_type == 'circle':
-            for i in range(self.init[0]):
-                for j in range(self.init[1]):
-                    r2 = self.xvalues[i] ** 2 + self.yvalues[j] ** 2
-                    me.values[i, j] = np.tanh((self.params.radius - np.sqrt(r2)) / (np.sqrt(2) * self.params.eps))
+            xv, yv = np.meshgrid(self.xvalues, self.yvalues, indexing='ij')
+            me.values[:, :] = np.tanh((self.params.radius - np.sqrt(xv ** 2 + yv ** 2)) /
+                                      (np.sqrt(2) * self.params.eps))
         elif self.params.init_type == 'checkerboard':
-            xv, yv = np.meshgrid(self.xvalues, self.xvalues)
+            xv, yv = np.meshgrid(self.xvalues, self.yvalues, indexing='ij')
             me.values[:, :] = np.sin(2.0 * np.pi * xv) * np.sin(2.0 * np.pi * yv)
         elif self.params.init_type == 'random':
             me.values[:, :] = np.random.uniform(-1, 1, self.init)

pySDC/playgrounds/parallel/AllenCahn_contracting_circle_FFT.py

@@ -10,7 +10,9 @@
 from pySDC.implementations.problem_classes.AllenCahn_2D_parFFT import allencahn2d_imex, allencahn2d_imex_stab
 from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
 from pySDC.implementations.transfer_classes.TransferMesh_FFT2D import mesh_to_mesh_fft2d
-from pySDC.projects.TOMS.AllenCahn_monitor import monitor
+from pySDC.playgrounds.parallel.AllenCahn_parallel_monitor import monitor
+
+import pySDC.helpers.plot_helper as plt_helper
 
 
 # http://www.personal.psu.edu/qud2/Res/Pre/dz09sisc.pdf
@@ -98,7 +100,7 @@ def run_SDC_variant(variant=None):
 
     # setup parameters "in time"
     t0 = 0.0
-    Tend = 0.032
+    Tend = 0.02
 
     # set MPI communicator
     comm = MPI.COMM_WORLD
@@ -128,6 +130,8 @@ def run_SDC_variant(variant=None):
                                                                         time_rank, time_size))
 
     description['problem_params']['comm'] = space_comm
+    # set level depending on rank
+    controller_params['logger_level'] = controller_params['logger_level'] if space_rank == 0 else 99
 
     # instantiate controller
     controller = controller_MPI(controller_params=controller_params, description=description, comm=time_comm)
@@ -136,9 +140,19 @@ def run_SDC_variant(variant=None):
     P = controller.S.levels[0].prob
     uinit = P.u_exact(t0)
 
+    # if time_rank == 0:
+    #     plt_helper.plt.imshow(uinit.values)
+    #     plt_helper.savefig(f'uinit_{space_rank}', save_pdf=False, save_pgf=False, save_png=True)
+    # exit()
+
     # call main function to get things done...
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
+    if time_rank == 0:
+        plt_helper.plt.imshow(uend.values)
+        plt_helper.savefig(f'uend_{space_rank}', save_pdf=False, save_pgf=False, save_png=True)
+    # exit()
+
     rank = comm.Get_rank()
 
     # filter statistics by variant (number of iterations)

pySDC/playgrounds/parallel/AllenCahn_parallel_monitor.py

@@ -0,0 +1,75 @@
+import numpy as np
+
+from pySDC.core.Hooks import hooks
+
+
+class monitor(hooks):
+
+    def __init__(self):
+        """
+        Initialization of Allen-Cahn monitoring
+        """
+        super(monitor, self).__init__()
+
+        self.init_radius = None
+
+    def pre_run(self, step, level_number):
+        """
+        Overwrite standard pre run hook
+
+        Args:
+            step (pySDC.Step.step): the current step
+            level_number (int): the current level number
+        """
+        super(monitor, self).pre_run(step, level_number)
+        L = step.levels[0]
+
+        c = np.count_nonzero(L.u[0].values > 0.0)
+        radius = np.sqrt(c / np.pi) * L.prob.dx
+
+        radius1 = 0
+        rows, cols = np.where(L.u[0].values > 0.0)
+        for r in rows:
+            radius1 = max(radius1, abs(L.prob.xvalues[r]))
+
+        # rows1 = np.where(L.u[0].values[int((L.prob.init[0]) / 2), :int((L.prob.init[0]) / 2)] > -0.99)
+        # rows2 = np.where(L.u[0].values[int((L.prob.init[0]) / 2), :int((L.prob.init[0]) / 2)] < 0.99)
+        # interface_width = (rows2[0][-1] - rows1[0][0]) * L.prob.dx / L.prob.params.eps
+
+        self.init_radius = L.prob.params.radius
+
+        if L.time == 0.0:
+            self.add_to_stats(process=step.status.slot, time=L.time, level=-1, iter=step.status.iter,
+                              sweep=L.status.sweep, type='computed_radius', value=radius)
+            self.add_to_stats(process=step.status.slot, time=L.time, level=-1, iter=step.status.iter,
+                              sweep=L.status.sweep, type='exact_radius', value=self.init_radius)
+            # self.add_to_stats(process=step.status.slot, time=L.time, level=-1, iter=step.status.iter,
+            #                   sweep=L.status.sweep, type='interface_width', value=interface_width)
+
+    def post_step(self, step, level_number):
+        """
+        Overwrite standard post step hook
+
+        Args:
+            step (pySDC.Step.step): the current step
+            level_number (int): the current level number
+        """
+        super(monitor, self).post_step(step, level_number)
+
+        # some abbreviations
+        L = step.levels[0]
+
+        c = np.count_nonzero(L.uend.values >= 0.0)
+        radius = np.sqrt(c / np.pi) * L.prob.dx
+
+        exact_radius = np.sqrt(max(self.init_radius ** 2 - 2.0 * (L.time + L.dt), 0))
+        # rows1 = np.where(L.uend.values[int((L.prob.init[0]) / 2), :int((L.prob.init[0]) / 2)] > -0.99)
+        # rows2 = np.where(L.uend.values[int((L.prob.init[0]) / 2), :int((L.prob.init[0]) / 2)] < 0.99)
+        # interface_width = (rows2[0][-1] - rows1[0][0]) * L.prob.dx / L.prob.params.eps
+
+        self.add_to_stats(process=step.status.slot, time=L.time + L.dt, level=-1, iter=step.status.iter,
+                          sweep=L.status.sweep, type='computed_radius', value=radius)
+        self.add_to_stats(process=step.status.slot, time=L.time + L.dt, level=-1, iter=step.status.iter,
+                          sweep=L.status.sweep, type='exact_radius', value=exact_radius)
+        # self.add_to_stats(process=step.status.slot, time=L.time + L.dt, level=-1, iter=step.status.iter,
+        #                   sweep=L.status.sweep, type='interface_width', value=interface_width)