added test equation

Author: pancetta

pySDC/implementations/problem_classes/TestEquation_0D.py

@@ -0,0 +1,110 @@
+from __future__ import division
+import numpy as np
+import scipy.sparse as sp
+from scipy.sparse.linalg import splu
+
+from pySDC.core.Problem import ptype
+from pySDC.core.Errors import ParameterError, ProblemError
+
+
+# noinspection PyUnusedLocal
+class testequation0d(ptype):
+    """
+    Example implementing a bundle of test equations at once (via diagonal matrix)
+
+    Attributes:
+        A: digonal matrix containing the parameters
+    """
+
+    def __init__(self, problem_params, dtype_u, dtype_f):
+        """
+        Initialization routine
+
+        Args:
+            problem_params (dict): custom parameters for the example
+            dtype_u: mesh data type for solution
+            dtype_f: mesh data type for RHS
+        """
+
+        # these parameters will be used later, so assert their existence
+        essential_keys = ['lambdas', 'u0']
+        for key in essential_keys:
+            if key not in problem_params:
+                msg = 'need %s to instantiate problem, only got %s' % (key, str(problem_params.keys()))
+                raise ParameterError(msg)
+
+        assert not any(isinstance(i, list) for i in problem_params['lambdas']), \
+            'ERROR: expect flat list here, got %s' % problem_params['lambdas']
+        nvars = len(problem_params['lambdas'])
+        assert nvars > 0, 'ERROR: expect at least one lambda parameter here'
+
+        # invoke super init, passing number of dofs, dtype_u and dtype_f
+        super(testequation0d, self).__init__(init=nvars, dtype_u=dtype_u, dtype_f=dtype_f,
+                                             params=problem_params)
+
+        self.A = self.__get_A(self.params.lambdas)
+
+    @staticmethod
+    def __get_A(lambdas):
+        """
+        Helper function to assemble FD matrix A in sparse format
+
+        Args:
+            lambdas (list): list of lambda parameters
+
+        Returns:
+            scipy.sparse.csc_matrix: diagonal matrix A in CSC format
+        """
+
+        A = sp.diags(lambdas)
+        return A
+
+    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)
+        f.values = self.A.dot(u.values)
+        return f
+
+    def solve_system(self, rhs, factor, u0, t):
+        """
+        Simple linear solver for (I-factor*A)u = rhs
+
+        Args:
+            rhs (dtype_f): right-hand side for the linear system
+            factor (float): abbrev. for the local stepsize (or any other factor required)
+            u0 (dtype_u): initial guess for the iterative solver
+            t (float): current time (e.g. for time-dependent BCs)
+
+        Returns:
+            dtype_u: solution as mesh
+        """
+
+        me = self.dtype_u(self.init)
+        L = splu(sp.eye(self.params.nvars, format='csc') - factor * self.A)
+        me.values = L.solve(rhs.values)
+        return me
+
+    def u_exact(self, t):
+        """
+        Routine to compute the exact solution at time t
+
+        Args:
+            t (float): current time
+
+        Returns:
+            dtype_u: exact solution
+        """
+
+        me = self.dtype_u(self.init)
+        me.values = self.params.u0 * np.exp(t * np.array(self.params.lambdas))
+        return me

pySDC/projects/matrixPFASST/allinclusive_matrix_nonMPI.py

@@ -5,6 +5,7 @@
 from pySDC.implementations.controller_classes.allinclusive_multigrid_nonMPI import allinclusive_multigrid_nonMPI
 
 from pySDC.implementations.datatype_classes.mesh import mesh
+from pySDC.implementations.datatype_classes.complex_mesh import mesh as cmesh
 from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
 from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
 
@@ -26,9 +27,9 @@ def __init__(self, num_procs, controller_params, description):
            description: all the parameters to set up the rest (levels, problems, transfer, ...)
        """
 
-        assert description['dtype_u'] is mesh, \
+        assert description['dtype_u'] is mesh or cmesh, \
             'ERROR: matrix version will only work with mesh data type for u, got %s' % description['dtype_u']
-        assert description['dtype_f'] is mesh, \
+        assert description['dtype_f'] is mesh or cmesh, \
             'ERROR: matrix version will only work with mesh data type for f, got %s' % description['dtype_f']
         assert description['sweeper_class'] is generic_implicit, \
             'ERROR: matrix version will only work with generic_implicit sweeper, got %s' % description['sweeper_class']

pySDC/projects/matrixPFASST/compare_to_propagator.py

@@ -2,10 +2,13 @@
 
 from pySDC.implementations.problem_classes.HeatEquation_1D_FD import heat1d
 from pySDC.implementations.problem_classes.AdvectionEquation_1D_FD import advection1d
+from pySDC.implementations.problem_classes.TestEquation_0D import testequation0d
 from pySDC.implementations.datatype_classes.mesh import mesh
+from pySDC.implementations.datatype_classes.complex_mesh import mesh as cmesh
 from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
 from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
 from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
+from pySDC.implementations.transfer_classes.TransferMesh_NoCoarse import mesh_to_mesh as mesh_to_mesh_nocoarse
 from pySDC.projects.matrixPFASST.allinclusive_matrix_nonMPI import allinclusive_matrix_nonMPI
 
 from pySDC.helpers.stats_helper import filter_stats, sort_stats
@@ -28,7 +31,7 @@ def diffusion_setup(par=0.0):
     sweeper_params = dict()
     sweeper_params['collocation_class'] = CollGaussRadau_Right
     sweeper_params['num_nodes'] = 3
-    sweeper_params['QI'] = 'LU'  # For the IMEX sweeper, the LU-trick can be activated for the implicit part
+    sweeper_params['QI'] = 'LU'
     sweeper_params['spread'] = True
 
     # initialize problem parameters
@@ -84,7 +87,7 @@ def advection_setup(par=0.0):
     sweeper_params = dict()
     sweeper_params['collocation_class'] = CollGaussRadau_Right
     sweeper_params['num_nodes'] = [3]
-    sweeper_params['QI'] = ['LU']  # For the IMEX sweeper, the LU-trick can be activated for the implicit part
+    sweeper_params['QI'] = ['LU']
     sweeper_params['spread'] = True
 
     # initialize problem parameters
@@ -126,6 +129,73 @@ def advection_setup(par=0.0):
     return description, controller_params
 
 
+def testequation_setup():
+    """
+    Setup routine for the test equation
+
+    Args:
+        par (float): parameter for controlling stiffness
+    """
+    # initialize level parameters
+    level_params = dict()
+    level_params['restol'] = 1E-08
+    level_params['dt'] = 0.25
+    level_params['nsweeps'] = [3]
+
+    # initialize sweeper parameters
+    sweeper_params = dict()
+    sweeper_params['collocation_class'] = CollGaussRadau_Right
+    sweeper_params['num_nodes'] = [3]
+    sweeper_params['QI'] = 'LU'
+    sweeper_params['spread'] = True
+
+    # initialize problem parameters
+    problem_params = dict()
+    problem_params['u0'] = 1.0  # initial value (for all instances)
+    # use single values like this...
+    problem_params['lambdas'] = [[-1.0]]
+    # .. or a list of values like this ...
+    problem_params['lambdas'] = [[-1.0, -2.0, 1j, -1j]]
+    # .. or a whole block of values like this
+    ilim_left = -11
+    ilim_right = 0
+    rlim_left = 0
+    rlim_right = 11
+    ilam = 1j * np.logspace(ilim_left, ilim_right, 11)
+    rlam = -1 * np.logspace(rlim_left, rlim_right, 11)
+    lambdas = []
+    for rl in rlam:
+        for il in ilam:
+            lambdas.append(rl + il)
+    problem_params['lambdas'] = [lambdas]
+    # note: PFASST will do all of those at once, but without interaction (realized via diagonal matrix).
+    # The propagation matrix will be diagonal too, corresponding to the respective lambda value.
+
+    # initialize step parameters
+    step_params = dict()
+    step_params['maxiter'] = 50
+
+    # initialize controller parameters
+    controller_params = dict()
+    controller_params['logger_level'] = 30
+    controller_params['predict'] = False
+
+    # fill description dictionary for easy step instantiation
+    description = dict()
+    description['problem_class'] = testequation0d  # pass problem class
+    description['problem_params'] = problem_params  # pass problem parameters
+    description['dtype_u'] = cmesh  # pass data type for u
+    description['dtype_f'] = cmesh  # pass data type for f
+    description['sweeper_class'] = generic_implicit  # pass sweeper
+    description['sweeper_params'] = sweeper_params  # pass sweeper parameters
+    description['level_params'] = level_params  # pass level parameters
+    description['step_params'] = step_params  # pass step parameters
+    description['space_transfer_class'] = mesh_to_mesh_nocoarse  # pass spatial transfer class
+    description['space_transfer_params'] = dict()  # pass paramters for spatial transfer
+
+    return description, controller_params
+
+
 def compare_controllers(type=None, par=0.0, f=None):
     """
     A simple test program to compare PFASST runs with matrix-based and matrix-free controllers
@@ -144,6 +214,8 @@ def compare_controllers(type=None, par=0.0, f=None):
         description, controller_params = diffusion_setup(par)
     elif type == 'advection':
         description, controller_params = advection_setup(par)
+    elif type == 'testequation':
+        description, controller_params = testequation_setup()
     else:
         raise ValueError('No valis setup type provided, aborting..')
 
@@ -196,12 +268,10 @@ def compare_controllers(type=None, par=0.0, f=None):
 
 def main():
 
-    par_list = [1E-02, 1.0, 1E+02]
-
     f = open('comparison_matrix_vs_propagator_detail.txt', 'a')
-    for par in par_list:
-        compare_controllers(type='diffusion', par=par, f=f)
-        compare_controllers(type='advection', par=par, f=f)
+    # compare_controllers(type='diffusion', par=1E-00, f=f)
+    # compare_controllers(type='advection', par=1E-00, f=f)
+    compare_controllers(type='testequation', par=0.0, f=f)
     f.close()