Adding logistic equation

Author: pancetta

pySDC/implementations/problem_classes/LogisticEquation.py

@@ -0,0 +1,123 @@
+import numpy as np
+
+from pySDC.core.Errors import ParameterError, ProblemError
+from pySDC.core.Problem import ptype
+from pySDC.implementations.datatype_classes.mesh import mesh
+
+
+# noinspection PyUnusedLocal
+class logistics_equation(ptype):
+    """
+    Example implementing the logistic equation, taken from
+    https://www-users.cse.umn.edu/~olver/ln_/odq.pdf (Example 2.2)
+    """
+
+    def __init__(self, problem_params, dtype_u=mesh, dtype_f=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 = ['u0', 'lam', 'newton_maxiter', 'newton_tol', 'direct']
+        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)
+        problem_params['nvars'] = 1
+
+        if 'stop_at_nan' not in problem_params:
+            problem_params['stop_at_nan'] = True
+
+        # invoke super init, passing dtype_u and dtype_f, plus setting number of elements to 2
+        super(logistics_equation, self).__init__((problem_params['nvars'], None, np.dtype('float64')),
+                                                 dtype_u, dtype_f, problem_params)
+
+    def u_exact(self, t):
+        """
+        Exact solution
+
+        Args:
+            t (float): current time
+        Returns:
+            dtype_u: mesh type containing the values
+        """
+
+        me = self.dtype_u(self.init)
+        me[:] = self.params.u0 * np.exp(self.params.lam * t) / (1 - self.params.u0 +
+                                                                self.params.u0 * np.exp(self.params.lam * t))
+        return me
+
+    def eval_f(self, u, t):
+        """
+        Routine to compute the RHS
+
+        Args:
+            u (dtype_u): the current values
+            t (float): current time (not used here)
+        Returns:
+            dtype_f: RHS, 1 component
+        """
+
+        f = self.dtype_f(self.init)
+        f[:] = self.params.lam * u * (1 - u)
+        return f
+
+    def solve_system(self, rhs, dt, u0, t):
+        """
+        Simple Newton solver for the nonlinear equation
+
+        Args:
+            rhs (dtype_f): right-hand side for the nonlinear system
+            dt (float): abbrev. for the node-to-node 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 u
+        """
+        # create new mesh object from u0 and set initial values for iteration
+        u = self.dtype_u(u0)
+
+        if self.params.direct:
+
+            d = (1 - dt * self.params.lam) ** 2 + 4 * dt * self.params.lam * rhs
+            u = (- (1 - dt * self.params.lam) + np.sqrt(d)) / (2 * dt * self.params.lam)
+            return u
+
+        else:
+
+            # start newton iteration
+            n = 0
+            res = 99
+            while n < self.params.newton_maxiter:
+
+                # form the function g with g(u) = 0
+                g = u - dt * self.params.lam * u * (1 - u) - rhs
+
+                # if g is close to 0, then we are done
+                res = np.linalg.norm(g, np.inf)
+                if res < self.params.newton_tol or np.isnan(res):
+                    break
+
+                # assemble dg/du
+                dg = 1 - dt * self.params.lam * (1 - 2 * u)
+                # newton update: u1 = u0 - g/dg
+                u -= 1.0 / dg * g
+
+                # 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:
+                raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
+
+            return u

pySDC/playgrounds/Gander/matrix_based.py

@@ -9,10 +9,10 @@ def iteration_vs_estimate():
     K = 10
     swee = sweeper({'collocation_class': CollGaussRadau_Right, 'num_nodes': M})
     Q = swee.coll.Qmat[1:, 1:]
-    # Qd = swee.get_Qdelta_implicit(swee.coll, 'IE')[1:, 1:]
-    Qd = swee.get_Qdelta_implicit(swee.coll, 'LU')[1:, 1:]
-    Qd = swee.get_Qdelta_explicit(swee.coll, 'EE')[1:, 1:]
-    print(swee.coll.nodes)
+    Qd = swee.get_Qdelta_implicit(swee.coll, 'IE')[1:, 1:]
+    # Qd = swee.get_Qdelta_implicit(swee.coll, 'LU')[1:, 1:]
+    # Qd = swee.get_Qdelta_explicit(swee.coll, 'EE')[1:, 1:]
+    print(np.linalg.norm(Q-Qd, np.inf))
     exit()
     I = np.eye(M)
 

pySDC/playgrounds/Gander/testequation.py

@@ -18,13 +18,13 @@ def testequation_setup(prec_type=None, maxiter=None):
     # initialize level parameters
     level_params = dict()
     level_params['restol'] = 0.0
-    level_params['dt'] = 0.25
+    level_params['dt'] = 1.0
     level_params['nsweeps'] = [1]
 
     # initialize sweeper parameters
     sweeper_params = dict()
     sweeper_params['collocation_class'] = CollGaussRadau_Right
-    sweeper_params['num_nodes'] = [3]
+    sweeper_params['num_nodes'] = [5]
     sweeper_params['QI'] = prec_type
     sweeper_params['initial_guess'] = 'spread'
 
@@ -35,7 +35,7 @@ def testequation_setup(prec_type=None, maxiter=None):
     # problem_params['lambdas'] = [[-1.0]]
     # .. or a list of values like this ...
     # problem_params['lambdas'] = [[-1.0, -2.0, 1j, -1j]]
-    problem_params['lambdas'] = [[-1000]]
+    problem_params['lambdas'] = [[-1.0 + 0j]]
     # 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.
 
@@ -45,7 +45,7 @@ def testequation_setup(prec_type=None, maxiter=None):
 
     # initialize controller parameters
     controller_params = dict()
-    controller_params['logger_level'] = 30
+    controller_params['logger_level'] = 20
     controller_params['hook_class'] = error_output
 
     # fill description dictionary for easy step instantiation
@@ -64,11 +64,11 @@ def compare_preconditioners(f=None, list_of_k=None):
 
     # set time parameters
     t0 = 0.0
-    Tend = 1.0
+    Tend = 2.0
 
     for k in list_of_k:
 
-        description_IE, controller_params_IE = testequation_setup(prec_type='IE', maxiter=k)
+        description_IE, controller_params_IE = testequation_setup(prec_type='MIN3', maxiter=k)
         description_LU, controller_params_LU = testequation_setup(prec_type='LU', maxiter=k)
 
         out = f'\nWorking with maxiter = {k}'
@@ -114,7 +114,7 @@ def compare_preconditioners(f=None, list_of_k=None):
 def main():
 
     f = open('comparison_IE_vs_LU.txt', 'w')
-    compare_preconditioners(f=f, list_of_k=[1, 2, 3, 4])
+    compare_preconditioners(f=f, list_of_k=[1])
     f.close()
 
 

pySDC/playgrounds/ODEs/logistics_playground.py

@@ -0,0 +1,68 @@
+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.LogisticEquation import logistics_equation
+from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
+from pySDC.playgrounds.ODEs.trajectory_HookClass import trajectories
+
+
+def main():
+    """
+    Van der Pol's oscillator inc. visualization
+    """
+    # initialize level parameters
+    level_params = dict()
+    level_params['restol'] = 1E-10
+    level_params['dt'] = 1.0
+
+    # initialize sweeper parameters
+    sweeper_params = dict()
+    sweeper_params['collocation_class'] = CollGaussRadau_Right
+    sweeper_params['num_nodes'] = 3
+    sweeper_params['QI'] = 'LU'
+
+    # initialize problem parameters
+    problem_params = dict()
+    problem_params['newton_tol'] = 1E-12
+    problem_params['newton_maxiter'] = 50
+    problem_params['lam'] = 20
+    problem_params['direct'] = True
+    problem_params['u0'] = 0.5
+
+    # initialize step parameters
+    step_params = dict()
+    step_params['maxiter'] = 20
+
+    # initialize controller parameters
+    controller_params = dict()
+    # controller_params['hook_class'] = trajectories
+    controller_params['logger_level'] = 20
+
+    # Fill description dictionary for easy hierarchy creation
+    description = dict()
+    description['problem_class'] = logistics_equation
+    description['problem_params'] = problem_params
+    description['sweeper_class'] = generic_implicit
+    description['sweeper_params'] = sweeper_params
+    description['level_params'] = level_params
+    description['step_params'] = step_params
+
+    # instantiate the controller
+    controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
+
+    # set time parameters
+    t0 = 0.0
+    Tend = level_params['dt']
+
+    # 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)
+
+    uex = P.u_exact(Tend)
+    print(abs(uend - uex))
+
+
+if __name__ == "__main__":
+    main()