some improvements for AC FFT example + parallel playground

Author: pancetta

pySDC/implementations/problem_classes/AllenCahn_2D_FFT.py

@@ -15,8 +15,8 @@ class allencahn2d_imex(ptype):
         xvalues: grid points in space
         dx: mesh width
         lap: spectral operator for Laplacian
-        fft_object: planned FFT for forward transformation
-        ifft_object: planned IFFT for backward transformation
+        rfft_object: planned real FFT for forward transformation
+        irfft_object: planned IFFT for backward transformation
     """
 
     def __init__(self, problem_params, dtype_u=mesh, dtype_f=rhs_imex_mesh):
@@ -57,28 +57,25 @@ def __init__(self, problem_params, dtype_u=mesh, dtype_f=rhs_imex_mesh):
         self.xvalues = np.array([i * self.dx - self.params.L / 2.0 for i in range(self.params.nvars[0])])
 
         kx = np.zeros(self.init[0])
-        ky = np.zeros(self.init[1])
+        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):
+        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)
-        for i in range(int(self.init[1] / 2) + 1, self.init[1]):
-            ky[i] = 2 * np.pi / self.params.L * (-self.init[1] + i)
 
-        self.lap = np.zeros(self.init)
+        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]):
+            for j in range(self.init[1] // 2 + 1):
                 self.lap[i, j] = -kx[i] ** 2 - ky[j] ** 2
 
-        # TODO: cleanup and move to real-valued FFT
-        fft_in = pyfftw.empty_aligned(self.init, dtype='complex128')
-        fft_out = pyfftw.empty_aligned(self.init, dtype='complex128')
-        ifft_in = pyfftw.empty_aligned(self.init, dtype='complex128')
-        ifft_out = pyfftw.empty_aligned(self.init, dtype='complex128')
-        self.fft_object = pyfftw.FFTW(fft_in, fft_out, direction='FFTW_FORWARD', axes=(0, 1))
-        self.ifft_object = pyfftw.FFTW(ifft_in, ifft_out, direction='FFTW_BACKWARD', axes=(0, 1))
+        rfft_in = pyfftw.empty_aligned(self.init, dtype='float64')
+        fft_out = pyfftw.empty_aligned((self.init[0], self.init[1] // 2 + 1), dtype='complex128')
+        ifft_in = pyfftw.empty_aligned((self.init[0], self.init[1] // 2 + 1), dtype='complex128')
+        irfft_out = pyfftw.empty_aligned(self.init, dtype='float64')
+        self.rfft_object = pyfftw.FFTW(rfft_in, fft_out, direction='FFTW_FORWARD', axes=(0, 1))
+        self.irfft_object = pyfftw.FFTW(ifft_in, irfft_out, direction='FFTW_BACKWARD', axes=(0, 1))
 
     def eval_f(self, u, t):
         """
@@ -94,8 +91,8 @@ def eval_f(self, u, t):
 
         f = self.dtype_f(self.init)
         v = u.values.flatten()
-        tmp = self.lap * self.fft_object(u.values)
-        f.impl.values[:] = np.real(self.ifft_object(tmp))
+        tmp = self.lap * self.rfft_object(u.values)
+        f.impl.values[:] = np.real(self.irfft_object(tmp))
         if self.params.eps > 0:
             f.expl.values = 1.0 / self.params.eps ** 2 * v * (1.0 - v ** self.params.nu)
             f.expl.values = f.expl.values.reshape(self.params.nvars)
@@ -117,8 +114,8 @@ def solve_system(self, rhs, factor, u0, t):
 
         me = self.dtype_u(self.init)
 
-        tmp = self.fft_object(rhs.values) / (1.0 - factor * self.lap)
-        me.values[:] = np.real(self.ifft_object(tmp))
+        tmp = self.rfft_object(rhs.values) / (1.0 - factor * self.lap)
+        me.values[:] = np.real(self.irfft_object(tmp))
 
         return me
 
@@ -149,3 +146,85 @@ def u_exact(self, t):
             raise NotImplementedError('type of initial value not implemented, got %s' % self.params.init_type)
 
         return me
+
+
+class allencahn2d_imex_stab(allencahn2d_imex):
+    """
+    Example implementing Allen-Cahn equation in 2D using FFTs for solving linear parts, IMEX time-stepping with
+    stabilized splitting
+
+    Attributes:
+        xvalues: grid points in space
+        dx: mesh width
+        lap: spectral operator for Laplacian
+        rfft_object: planned real FFT for forward transformation
+        irfft_object: planned IFFT for backward transformation
+    """
+
+    def __init__(self, problem_params, dtype_u=mesh, dtype_f=rhs_imex_mesh):
+        """
+        Initialization routine
+
+        Args:
+            problem_params (dict): custom parameters for the example
+            dtype_u: mesh data type (will be passed to parent class)
+            dtype_f: mesh data type wuth implicit and explicit parts (will be passed to parent class)
+        """
+        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
+
+    def eval_f(self, u, t):
+        """
+        Routine to evaluate the RHS
+
+        Args:
+            u (dtype_u): current values
+            t (float): current time
+
+        Returns:
+            dtype_f: the RHS
+        """
+
+        f = self.dtype_f(self.init)
+        v = u.values.flatten()
+        tmp = self.lap * self.rfft_object(u.values)
+        f.impl.values[:] = np.real(self.irfft_object(tmp))
+        if self.params.eps > 0:
+            f.expl.values = 1.0 / self.params.eps ** 2 * v * (1.0 - v ** self.params.nu) + \
+                            2.0 / self.params.eps ** 2 * v
+            f.expl.values = f.expl.values.reshape(self.params.nvars)
+        return f
+
+    def solve_system(self, rhs, factor, u0, t):
+        """
+        Simple FFT solver for the diffusion part
+
+        Args:
+            rhs (dtype_f): right-hand side for the linear system
+            factor (float) : abbrev. for the node-to-node stepsize (or any other factor required)
+            u0 (dtype_u): initial guess for the iterative solver (not used here so far)
+            t (float): current time (e.g. for time-dependent BCs)
+
+        Returns:
+            dtype_u: solution as mesh
+        """
+
+        me = self.dtype_u(self.init)
+
+        tmp = self.rfft_object(rhs.values) / (1.0 - factor * self.lap)
+        me.values[:] = np.real(self.irfft_object(tmp))
+
+        return me

pySDC/playgrounds/fft/AllenCahn_contracting_circle_FFT.py

@@ -8,7 +8,7 @@
 from pySDC.helpers.stats_helper import filter_stats, sort_stats
 from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
-from pySDC.implementations.problem_classes.AllenCahn_2D_FFT import allencahn2d_imex
+from pySDC.implementations.problem_classes.AllenCahn_2D_FFT 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
@@ -32,7 +32,7 @@ def setup_parameters():
     level_params = dict()
     level_params['restol'] = 1E-08
     level_params['dt'] = 1E-03
-    level_params['nsweeps'] = [2, 1]
+    level_params['nsweeps'] = [3, 1]
 
     # initialize sweeper parameters
     sweeper_params = dict()
@@ -46,13 +46,13 @@ def setup_parameters():
     problem_params = dict()
     problem_params['nu'] = 2
     problem_params['L'] = 1.0
-    problem_params['nvars'] = [(128, 128), (32, 32)]
-    problem_params['eps'] = [0.04, 0.04]
+    problem_params['nvars'] = [(256, 256), (64, 64)]
+    problem_params['eps'] = [0.04, 0.16]
     problem_params['radius'] = 0.25
 
     # initialize step parameters
     step_params = dict()
-    step_params['maxiter'] = 20
+    step_params['maxiter'] = 50
 
     # initialize controller parameters
     controller_params = dict()
@@ -72,13 +72,12 @@ def setup_parameters():
     return description, controller_params
 
 
-def run_SDC_variant(variant=None, inexact=False):
+def run_SDC_variant(variant=None):
     """
     Routine to run particular SDC variant
 
     Args:
         variant (str): string describing the variant
-        inexact (bool): flag to use inexact nonlinear solve (or nor)
 
     Returns:
         timing (float)
@@ -92,6 +91,9 @@ def run_SDC_variant(variant=None, inexact=False):
     if variant == 'semi-implicit':
         description['problem_class'] = allencahn2d_imex
         description['sweeper_class'] = imex_1st_order
+    elif variant == 'semi-implicit-stab':
+        description['problem_class'] = allencahn2d_imex_stab
+        description['sweeper_class'] = imex_1st_order
     else:
         raise NotImplemented('Wrong variant specified, got %s' % variant)
 
@@ -100,8 +102,7 @@ def run_SDC_variant(variant=None, inexact=False):
     Tend = 0.032
 
     # instantiate controller
-    controller = controller_nonMPI(num_procs=8, controller_params=controller_params,
-                                               description=description)
+    controller = controller_nonMPI(num_procs=8, controller_params=controller_params, description=description)
 
     # get initial values on finest level
     P = controller.MS[0].levels[0].prob
@@ -207,13 +208,13 @@ def show_results(fname, cwd=''):
 
         xcoords = [item0[0] for item0 in computed_radii]
         radii = [item0[1] for item0 in computed_radii]
-        if key[0] + ' ' + key[1] == 'semi-implicit exact':
+        if key[0] + ' ' + key[1] == 'semi-implicit-stab exact':
             ax.plot(xcoords, radii, label=(key[0] + ' ' + key[1]).replace('_v2', ' mod.'))
 
         exact_radii = sort_stats(filter_stats(item, type='exact_radius'), sortby='time')
 
-        diff = np.array([abs(item0[1] - item1[1]) for item0, item1 in zip(exact_radii, computed_radii)])
-        max_pos = int(np.argmax(diff))
+        # diff = np.array([abs(item0[1] - item1[1]) for item0, item1 in zip(exact_radii, computed_radii)])
+        # max_pos = int(np.argmax(diff))
         # assert max(diff) < 0.07, 'ERROR: computed radius is too far away from exact radius, got %s' % max(diff)
         # assert 0.028 < computed_radii[max_pos][0] < 0.03, \
         #     'ERROR: largest difference is at wrong time, got %s' % computed_radii[max_pos][0]
@@ -278,9 +279,9 @@ def main(cwd=''):
 
     # Loop over variants, exact and inexact solves
     results = {}
-    for variant in ['semi-implicit']:
+    for variant in ['semi-implicit-stab']:
 
-        results[(variant, 'exact')] = run_SDC_variant(variant=variant, inexact=False)
+        results[(variant, 'exact')] = run_SDC_variant(variant=variant)
 
     # dump result
     fname = 'data/results_SDC_variants_AllenCahn_1E-03'