Bayreuth 1D SDC code is up and running, very stable

Author: pancetta

pySDC/implementations/problem_classes/AllenCahn_1D_FD.py

@@ -5,11 +5,10 @@
 
 from pySDC.core.Errors import ParameterError, ProblemError
 from pySDC.core.Problem import ptype
-from pySDC.implementations.datatype_classes.mesh import mesh
+from pySDC.implementations.datatype_classes.mesh import mesh, rhs_imex_mesh
 
 
-# noinspection PyUnusedLocal
-class allencahn_front(ptype):
+class allencahn_front_fullyimplicit(ptype):
     """
     Example implementing the Allen-Cahn equation in 1D with finite differences and inhomogeneous Dirichlet-BC,
     with driving force, 0-1 formulation (Bayreuth example)
@@ -44,7 +43,7 @@ def __init__(self, problem_params, dtype_u=mesh, dtype_f=mesh):
             problem_params['stop_at_nan'] = True
 
         # invoke super init, passing number of dofs, dtype_u and dtype_f
-        super(allencahn_front, self).__init__(problem_params['nvars'], dtype_u, dtype_f, problem_params)
+        super(allencahn_front_fullyimplicit, self).__init__(problem_params['nvars'], dtype_u, dtype_f, problem_params)
 
         # compute dx and get discretization matrix A
         self.dx = (self.params.interval[1] - self.params.interval[0]) / (self.params.nvars + 1)
@@ -92,7 +91,6 @@ def solve_system(self, rhs, factor, u0, t):
         """
 
         u = self.dtype_u(u0).values
-        z = self.dtype_u(self.init, val=0.0).values
         eps2 = self.params.eps ** 2
         dw = self.params.dw
 
@@ -188,17 +186,96 @@ def u_exact(self, t):
         return me
 
 
-# noinspection PyUnusedLocal
-class allencahn_front_finel(allencahn_front):
+class allencahn_front_semiimplicit(allencahn_front_fullyimplicit):
     """
     Example implementing the Allen-Cahn equation in 1D with finite differences and inhomogeneous Dirichlet-BC,
-    with driving force, 0-1 formulation (Bayreuth example), Finel's trick/parametrization
+    with driving force, 0-1 formulation (Bayreuth example), semi-implicit time-stepping
 
     Attributes:
         A: second-order FD discretization of the 1D laplace operator
         dx: distance between two spatial nodes
     """
 
+    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 parent class)
+            dtype_f: mesh data type with implicit and explicit components (will be passed parent class)
+        """
+
+        # these parameters will be used later, so assert their existence
+        essential_keys = ['nvars', 'dw', 'eps', 'newton_maxiter', 'newton_tol', 'interval']
+        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)
+
+        # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on
+        if (problem_params['nvars'] + 1) % 2 != 0:
+            raise ProblemError('setup requires nvars = 2^p - 1')
+
+        if 'stop_at_nan' not in problem_params:
+            problem_params['stop_at_nan'] = True
+
+        # invoke super init, passing number of dofs, dtype_u and dtype_f
+        super(allencahn_front_semiimplicit, self).__init__(problem_params, dtype_u, dtype_f)
+
+    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
+        """
+        # set up boundary values to embed inner points
+        v = 3.0 * np.sqrt(2) * self.params.eps * self.params.dw
+        self.uext.values[0] = 0.5 * (1 + np.tanh((self.params.interval[0] - v * t) / (np.sqrt(2) * self.params.eps)))
+        self.uext.values[-1] = 0.5 * (1 + np.tanh((self.params.interval[1] - v * t) / (np.sqrt(2) * self.params.eps)))
+
+        self.uext.values[1:-1] = u.values[:]
+
+        f = self.dtype_f(self.init)
+        f.impl.values = self.A.dot(self.uext.values)[1:-1]
+        f.expl.values = - 2.0 / self.params.eps ** 2 * u.values * (1.0 - u.values) * (1.0 - 2 * u.values) - \
+            6.0 * self.params.dw * u.values * (1.0 - 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)
+        self.uext.values[0] = 0.0
+        self.uext.values[-1] = 0.0
+        self.uext.values[1:-1] = rhs.values[:]
+        me.values = spsolve(sp.eye(self.params.nvars + 2, format='csc') - factor * self.A, self.uext.values)[1:-1]
+        # me.values = spsolve(sp.eye(self.params.nvars, format='csc') - factor * self.A[1:-1, 1:-1], rhs.values)
+        return me
+
+
+class allencahn_front_finel(allencahn_front_fullyimplicit):
+    """
+    Example implementing the Allen-Cahn equation in 1D with finite differences and inhomogeneous Dirichlet-BC,
+    with driving force, 0-1 formulation (Bayreuth example), Finel's trick/parametrization
+    """
+
     # noinspection PyTypeChecker
     def solve_system(self, rhs, factor, u0, t):
         """

pySDC/playgrounds/Allen_Cahn/AllenCahn_Bayreuth_1D.py

@@ -8,8 +8,9 @@
 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_1D_FD import allencahn_front, allencahn_front_finel
+from pySDC.implementations.problem_classes.AllenCahn_1D_FD import allencahn_front_fullyimplicit, allencahn_front_semiimplicit
 from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
+from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
 from pySDC.playgrounds.Allen_Cahn.AllenCahn_monitor_Bayreuth import monitor
 
 
@@ -26,14 +27,14 @@ def setup_parameters():
 
     # initialize level parameters
     level_params = dict()
-    level_params['restol'] = 1E-08
-    level_params['dt'] = 1.0 / 4
+    level_params['restol'] = 1E-03
+    level_params['dt'] = 16.0 / 1
     level_params['nsweeps'] = 1
 
     # initialize sweeper parameters
     sweeper_params = dict()
     sweeper_params['collocation_class'] = CollGaussRadau_Right
-    sweeper_params['num_nodes'] = [5]
+    sweeper_params['num_nodes'] = [3]
     sweeper_params['Q1'] = ['LU']
     sweeper_params['Q2'] = ['LU']
     sweeper_params['QI'] = ['LU']
@@ -42,23 +43,23 @@ def setup_parameters():
 
     # This comes as read-in for the problem class
     problem_params = dict()
-    problem_params['nvars'] = 2047
+    problem_params['nvars'] = 127
     problem_params['dw'] = -0.04
     problem_params['eps'] = 0.04
     problem_params['newton_maxiter'] = 100
-    problem_params['newton_tol'] = 1E-08
+    problem_params['newton_tol'] = 1E-04
     problem_params['lin_tol'] = 1E-08
     problem_params['lin_maxiter'] = 100
     problem_params['radius'] = 0.25
     problem_params['interval'] = (-0.5, 0.5)
 
     # initialize step parameters
     step_params = dict()
-    step_params['maxiter'] = 20
+    step_params['maxiter'] = 50
 
     # initialize controller parameters
     controller_params = dict()
-    controller_params['logger_level'] = 30
+    controller_params['logger_level'] = 20
     controller_params['hook_class'] = monitor
 
     # fill description dictionary for easy step instantiation
@@ -70,6 +71,7 @@ def setup_parameters():
     description['level_params'] = level_params  # pass level parameters
     description['step_params'] = step_params  # pass step parameters
 
+
     return description, controller_params
 
 
@@ -90,21 +92,16 @@ def run_SDC_variant(variant=None, inexact=False):
 
     # add stuff based on variant
     if variant == 'fully-implicit':
-        # description['problem_class'] = allencahn_front
-        description['problem_class'] = allencahn_front_finel
+        description['problem_class'] = allencahn_front_fullyimplicit
+        # description['problem_class'] = allencahn_front_finel
         description['sweeper_class'] = generic_implicit
         if inexact:
-            description['problem_params']['newton_maxiter'] = 1
-    # elif variant == 'semi-implicit':
-    #     description['problem_class'] = allencahn_semiimplicit
-    #     description['sweeper_class'] = imex_1st_order
-    #     if inexact:
-    #         description['problem_params']['lin_maxiter'] = 10
-    # elif variant == 'semi-implicit_v2':
-    #     description['problem_class'] = allencahn_semiimplicit_v2
-    #     description['sweeper_class'] = imex_1st_order
-    #     if inexact:
-    #         description['problem_params']['newton_maxiter'] = 1
+            description['problem_params']['newton_maxiter'] = 200
+    elif variant == 'semi-implicit':
+        description['problem_class'] = allencahn_front_semiimplicit
+        description['sweeper_class'] = imex_1st_order
+        if inexact:
+            description['problem_params']['lin_maxiter'] = 10
     # elif variant == 'multi-implicit':
     #     description['problem_class'] = allencahn_multiimplicit
     #     description['sweeper_class'] = multi_implicit
@@ -127,7 +124,7 @@ def run_SDC_variant(variant=None, inexact=False):
 
     # setup parameters "in time"
     t0 = 0
-    Tend = 1.0
+    Tend = 32.0
 
     # instantiate controller
     controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
@@ -147,8 +144,8 @@ def run_SDC_variant(variant=None, inexact=False):
     # call main function to get things done...
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
-    plt_helper.plt.plot(uend.values)
-    plt_helper.savefig('uend', save_pdf=False, save_pgf=False, save_png=True)
+    # plt_helper.plt.plot(uend.values)
+    # plt_helper.savefig('uend', save_pdf=False, save_pgf=False, save_png=True)
     # exit()
 
     # filter statistics by variant (number of iterations)
@@ -320,9 +317,10 @@ def main(cwd=''):
     results = {}
     # for variant in ['multi-implicit', 'semi-implicit', 'fully-implicit', 'semi-implicit_v2', 'multi-implicit_v2']:
     for variant in ['fully-implicit']:
+    # for variant in ['semi-implicit']:
 
-        results[(variant, 'exact')] = run_SDC_variant(variant=variant, inexact=False)
-        # results[(variant, 'inexact')] = run_SDC_variant(variant=variant, inexact=True)
+        # results[(variant, 'exact')] = run_SDC_variant(variant=variant, inexact=False)
+        results[(variant, 'inexact')] = run_SDC_variant(variant=variant, inexact=True)
 
     # dump result
     # fname = 'data/results_SDC_variants_AllenCahn_1E-03'

pySDC/playgrounds/Allen_Cahn/AllenCahn_Bayreuth_1D_reference.py

@@ -0,0 +1,132 @@
+import timeit
+import numpy as np
+
+import pySDC.helpers.plot_helper as plt_helper
+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_1D_FD import allencahn_front_fullyimplicit
+from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
+from pySDC.implementations.sweeper_classes.explicit import explicit
+from pySDC.playgrounds.Allen_Cahn.AllenCahn_monitor_Bayreuth import monitor
+
+
+def setup_parameters():
+    """
+    Helper routine to fill in all relevant parameters
+
+    Note that this file will be used for all versions of SDC, containing more than necessary for each individual run
+
+    Returns:
+        description (dict)
+        controller_params (dict)
+    """
+
+    # initialize level parameters
+    level_params = dict()
+    level_params['restol'] = 1E-08
+    level_params['dt'] = 1.0 / (16.0 * 128.0 ** 2)
+    level_params['nsweeps'] = 1
+
+    # initialize sweeper parameters
+    sweeper_params = dict()
+    sweeper_params['collocation_class'] = CollGaussRadau_Right
+    sweeper_params['num_nodes'] = [1]
+    sweeper_params['Q1'] = ['LU']
+    sweeper_params['Q2'] = ['LU']
+    sweeper_params['QI'] = ['LU']
+    sweeper_params['QE'] = ['EE']
+    sweeper_params['spread'] = False
+
+    # This comes as read-in for the problem class
+    problem_params = dict()
+    problem_params['nvars'] = 127
+    problem_params['dw'] = -0.04
+    problem_params['eps'] = 0.04
+    problem_params['newton_maxiter'] = 100
+    problem_params['newton_tol'] = 1E-08
+    problem_params['lin_tol'] = 1E-08
+    problem_params['lin_maxiter'] = 100
+    problem_params['radius'] = 0.25
+    problem_params['interval'] = (-0.5, 0.5)
+
+    # initialize step parameters
+    step_params = dict()
+    step_params['maxiter'] = 20
+
+    # initialize controller parameters
+    controller_params = dict()
+    controller_params['logger_level'] = 30
+    controller_params['hook_class'] = monitor
+
+    # fill description dictionary for easy step instantiation
+    description = dict()
+    description['problem_class'] = allencahn_front_fullyimplicit
+    description['problem_params'] = problem_params  # pass problem parameters
+    description['sweeper_class'] = explicit
+    description['sweeper_params'] = sweeper_params  # pass sweeper parameters
+    description['level_params'] = level_params  # pass level parameters
+    description['step_params'] = step_params  # pass step parameters
+
+    return description, controller_params
+
+
+def run_SDC_variant(variant=None, inexact=False):
+    """
+    Routine to run particular SDC variant
+
+    Args:
+        variant (str): string describing the variant
+        inexact (bool): flag to use inexact nonlinear solve (or nor)
+
+    Returns:
+        results and statistics of the run
+    """
+
+    # load (incomplete) default parameters
+    description, controller_params = setup_parameters()
+
+    # add stuff based on variant
+    if variant == 'explicit':
+        description['problem_class'] = allencahn_front_fullyimplicit
+        description['sweeper_class'] = explicit
+    else:
+        raise NotImplemented('Wrong variant specified, got %s' % variant)
+
+    # setup parameters "in time"
+    t0 = 0
+    Tend = 1.0 / (16.0 * 128.0 ** 2)
+
+    # instantiate controller
+    controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
+
+    # get initial values on finest level
+    P = controller.MS[0].levels[0].prob
+    uinit = P.u_exact(t0)
+
+    # call main function to get things done...
+    # uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
+    wrapped = wrapper(controller.run, u0=uinit, t0=t0, Tend=Tend)
+    print(timeit.timeit(wrapped, number=10000) / 10000.0)
+
+
+def wrapper(func, *args, **kwargs):
+    def wrapped():
+        return func(*args, **kwargs)
+    return wrapped
+
+
+def main(cwd=''):
+    """
+    Main driver
+
+    Args:
+        cwd (str): current working directory (need this for testing)
+    """
+
+    # Loop over variants, exact and inexact solves
+    run_SDC_variant(variant='explicit')
+
+
+if __name__ == "__main__":
+    main()