more playing around

Author: pancetta

pySDC/implementations/problem_classes/AllenCahn_1D_FD.py

@@ -5,7 +5,7 @@
 
 from pySDC.core.Errors import ParameterError, ProblemError
 from pySDC.core.Problem import ptype
-from pySDC.implementations.datatype_classes.mesh import mesh, rhs_imex_mesh
+from pySDC.implementations.datatype_classes.mesh import mesh, rhs_imex_mesh, rhs_comp2_mesh
 
 
 class allencahn_front_fullyimplicit(ptype):
@@ -580,7 +580,7 @@ def __init__(self, problem_params, dtype_u=mesh, dtype_f=rhs_imex_mesh):
         # invoke super init, passing number of dofs, dtype_u and dtype_f
         super(allencahn_periodic_semiimplicit, self).__init__(problem_params, dtype_u, dtype_f)
 
-        self.A -= sp.eye(self.init) * 2.0 / self.params.eps ** 2
+        self.A -= sp.eye(self.init) * 0.0 / self.params.eps ** 2
 
     def solve_system(self, rhs, factor, u0, t):
         """
@@ -613,6 +613,139 @@ def eval_f(self, u, t):
         """
         f = self.dtype_f(self.init)
         f.impl.values = self.A.dot(u.values)
-        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) + 2.0 / self.params.eps ** 2 * u.values
+        f.expl.values = -2.0 / self.params.eps ** 2 * u.values * (1.0 - u.values) * (1.0 - 2.0 * u.values) - \
+            6.0 * self.params.dw * u.values * (1.0 - u.values) + 0.0 / self.params.eps ** 2 * u.values
         return f
+
+
+class allencahn_periodic_multiimplicit(allencahn_periodic_fullyimplicit):
+    """
+    Example implementing the Allen-Cahn equation in 1D with finite differences and periodic BC,
+    with driving force, 0-1 formulation (Bayreuth example)
+    """
+
+    def __init__(self, problem_params, dtype_u=mesh, dtype_f=rhs_comp2_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 (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', 'radius']
+        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']) % 2 != 0:
+            raise ProblemError('setup requires nvars = 2^p')
+
+        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_periodic_multiimplicit, self).__init__(problem_params, dtype_u, dtype_f)
+
+        self.A -= sp.eye(self.init) * 0.0 / self.params.eps ** 2
+
+    def solve_system_1(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(u0)
+        me.values = spsolve(sp.eye(self.params.nvars, format='csc') - factor * self.A, rhs.values)
+        return me
+
+    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.comp1.values = self.A.dot(u.values)
+        f.comp2.values = -2.0 / self.params.eps ** 2 * u.values * (1.0 - u.values) * (1.0 - 2.0 * u.values) - \
+            6.0 * self.params.dw * u.values * (1.0 - u.values) + 0.0 / self.params.eps ** 2 * u.values
+        return f
+
+    def solve_system_2(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
+        """
+
+        u = self.dtype_u(u0).values
+        eps2 = self.params.eps ** 2
+        dw = self.params.dw
+
+        Id = sp.eye(self.params.nvars)
+
+        # start newton iteration
+        n = 0
+        res = 99
+        while n < self.params.newton_maxiter:
+            # print(n)
+            # form the function g with g(u) = 0
+            g = u - rhs.values - factor * (- 2.0 / eps2 * u * (1.0 - u) * (1.0 - 2.0 * u) -
+                                           6.0 * dw * u * (1.0 - u) + 0.0 / self.params.eps ** 2 * u)
+
+            # if g is close to 0, then we are done
+            res = np.linalg.norm(g, np.inf)
+
+            if res < self.params.newton_tol:
+                break
+
+            # assemble dg
+            dg = Id - factor * (- 2.0 / eps2 * sp.diags(
+                (1.0 - u) * (1.0 - 2.0 * u) - u * ((1.0 - 2.0 * u) + 2.0 * (1.0 - u)), offsets=0) - 6.0 * dw * sp.diags(
+                (1.0 - u) - u, offsets=0) + 0.0 / self.params.eps ** 2 * Id)
+
+            # newton update: u1 = u0 - g/dg
+            u -= spsolve(dg, g)
+            # u -= gmres(dg, g, x0=z, tol=self.params.lin_tol)[0]
+            # increase iteration count
+            n += 1
+
+        if np.isnan(res) and self.params.stop_at_nan:
+            raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
+        elif np.isnan(res):
+            self.logger.warning('Newton got nan after %i iterations...' % n)
+
+        if n == self.params.newton_maxiter:
+            self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res))
+
+        self.newton_ncalls += 1
+        self.newton_itercount += n
+
+        me = self.dtype_u(self.init)
+        me.values = u
+
+        return me

pySDC/playgrounds/Allen_Cahn/AllenCahn_Bayreuth_periodic_1D.py

@@ -8,9 +8,12 @@
 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_periodic_fullyimplicit, allencahn_periodic_semiimplicit
+from pySDC.implementations.problem_classes.AllenCahn_1D_FD import allencahn_periodic_fullyimplicit, allencahn_periodic_semiimplicit, allencahn_periodic_multiimplicit
 from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
 from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
+from pySDC.implementations.sweeper_classes.multi_implicit import multi_implicit
+from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
+from pySDC.implementations.transfer_classes.TransferMesh_FFT import mesh_to_mesh_fft
 from pySDC.playgrounds.Allen_Cahn.AllenCahn_monitor_Bayreuth import monitor
 
 
@@ -29,24 +32,24 @@ def setup_parameters():
     level_params = dict()
     level_params['restol'] = 1E-08
     level_params['dt'] = 2E-02
-    level_params['nsweeps'] = 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']
     sweeper_params['QE'] = ['EE']
-    sweeper_params['spread'] = False
+    sweeper_params['spread'] = True
 
     # This comes as read-in for the problem class
     problem_params = dict()
-    problem_params['nvars'] = 128 * 8
-    problem_params['dw'] = -0.04
-    problem_params['eps'] = 0.04
-    problem_params['newton_maxiter'] = 100
+    problem_params['nvars'] = [128 * 8]#, 128 * 4]
+    problem_params['dw'] = [-0.04]
+    problem_params['eps'] = [0.04]
+    problem_params['newton_maxiter'] = 200
     problem_params['newton_tol'] = 1E-08
     problem_params['lin_tol'] = 1E-08
     problem_params['lin_maxiter'] = 100
@@ -55,7 +58,12 @@ def setup_parameters():
 
     # initialize step parameters
     step_params = dict()
-    step_params['maxiter'] = 50
+    step_params['maxiter'] = 20
+
+    # initialize space transfer parameters
+    space_transfer_params = dict()
+    space_transfer_params['rorder'] = 2
+    space_transfer_params['iorder'] = 6
 
     # initialize controller parameters
     controller_params = dict()
@@ -70,6 +78,8 @@ def setup_parameters():
     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_fft  # pass spatial transfer class
+    # description['space_transfer_params'] = space_transfer_params  # pass paramters for spatial transfer
 
     return description, controller_params
 
@@ -101,12 +111,12 @@ def run_SDC_variant(variant=None, inexact=False):
         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
-    #     if inexact:
-    #         description['problem_params']['newton_maxiter'] = 1
-    #         description['problem_params']['lin_maxiter'] = 10
+    elif variant == 'multi-implicit':
+        description['problem_class'] = allencahn_periodic_multiimplicit
+        description['sweeper_class'] = multi_implicit
+        if inexact:
+            description['problem_params']['newton_maxiter'] = 1
+            description['problem_params']['lin_maxiter'] = 10
     # elif variant == 'multi-implicit_v2':
     #     description['problem_class'] = allencahn_multiimplicit_v2
     #     description['sweeper_class'] = multi_implicit
@@ -123,7 +133,7 @@ def run_SDC_variant(variant=None, inexact=False):
 
     # setup parameters "in time"
     t0 = 0
-    Tend = 1.0# * description['level_params']['dt']
+    Tend = description['level_params']['dt']
 
     # instantiate controller
     controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
@@ -315,8 +325,9 @@ def main(cwd=''):
     # Loop over variants, exact and inexact solves
     results = {}
     # for variant in ['multi-implicit', 'semi-implicit', 'fully-implicit', 'semi-implicit_v2', 'multi-implicit_v2']:
-    for variant in ['fully-implicit']:
+    # for variant in ['fully-implicit']:
     # for variant in ['semi-implicit']:
+    for variant in ['multi-implicit']:
 
         results[(variant, 'exact')] = run_SDC_variant(variant=variant, inexact=False)
         # results[(variant, 'inexact')] = run_SDC_variant(variant=variant, inexact=True)

pySDC/playgrounds/parallel/AllenCahn_contracting_circle_FFT.py

@@ -32,7 +32,7 @@ def setup_parameters():
     # initialize level parameters
     level_params = dict()
     level_params['restol'] = 1E-08
-    level_params['dt'] = 1E-03
+    level_params['dt'] = 1E-02
     level_params['nsweeps'] = [3, 1]
 
     # initialize sweeper parameters
@@ -184,7 +184,7 @@ 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)