simple nonlinear ODE test case

pancetta · pancetta · commit e7c55450abf0 · 2021-09-16T10:37:22.000+02:00
diff --git a/pySDC/implementations/problem_classes/nonlinear_ODE_1.py b/pySDC/implementations/problem_classes/nonlinear_ODE_1.py
@@ -0,0 +1,114 @@
+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 nonlinear_ODE_1(ptype):
+    """
+    Example implementing some simple nonlinear ODE with a singularity in the derivative, taken from
+    https://www.osti.gov/servlets/purl/6111421 (Problem E-4)
+    """
+
+    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', 'newton_maxiter', 'newton_tol']
+        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(nonlinear_ODE_1, 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[:] = t - t ** 2 / 4
+        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[:] = np.sqrt(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)
+
+        # start newton iteration
+        n = 0
+        res = 99
+        while n < self.params.newton_maxiter:
+
+            # form the function g with g(u) = 0
+            g = u - dt * np.sqrt(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) / (2 * np.sqrt(1 - 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
diff --git a/pySDC/playgrounds/Gander/matrix_based.py b/pySDC/playgrounds/Gander/matrix_based.py
@@ -0,0 +1,122 @@
+from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
+from pySDC.core.Sweeper import sweeper
+import numpy as np
+import matplotlib.pyplot as plt
+
+def iteration_vs_estimate():
+
+    M = 5
+    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)
+    exit()
+    I = np.eye(M)
+
+    lam = -0.7
+    print(1/np.linalg.norm(Qd, np.inf), np.linalg.norm(np.linalg.inv(Qd), np.inf))
+    C = I - lam * Q
+    P = I - lam * Qd
+    Pinv = np.linalg.inv(P)
+
+    R = Pinv.dot(lam * (Q - Qd))
+    rho = max(abs(np.linalg.eigvals(R)))
+    infnorm = np.linalg.norm(R, np.inf)
+    twonorm = np.linalg.norm(R, 2)
+
+    # uex = np.exp(lam)
+    u0 = np.ones(M, dtype=np.complex128)
+    uex = np.linalg.inv(C).dot(u0)[-1]
+    u = u0.copy()
+    # u = np.random.rand(M)
+    res = u0 - C.dot(u)
+    err = []
+    resnorm = []
+    for k in range(K):
+        u += Pinv.dot(res)
+        res = u0 - C.dot(u)
+        err.append(abs(u[-1] - uex))
+        resnorm.append(np.linalg.norm(res, np.inf))
+        print(k, resnorm[-1], err[-1])
+
+    plt.figure()
+    plt.semilogy(range(K), err, 'o-', color='red', label='error')
+    plt.semilogy(range(K), resnorm, 'd-', color='orange', label='residual')
+    plt.semilogy(range(K), [err[0] * rho ** k for k in range(K)], '--', color='green', label='spectral')
+    plt.semilogy(range(K), [err[0] * infnorm ** k for k in range(K)], '--', color='blue', label='infnorm')
+    plt.semilogy(range(K), [err[0] * twonorm ** k for k in range(K)], '--', color='cyan', label='twonorm')
+    plt.semilogy(range(K), [err[0] * ((-1/lam) ** (1/M)) ** k for k in range(K)], 'x--', color='black', label='est')
+    plt.semilogy(range(K), [err[0] * (abs(lam) * np.linalg.norm(Q-Qd)) ** k for k in range(K)], 'x-.', color='black', label='est')
+    # plt.semilogy(range(K), [err[0] * (1/abs(lam) + np.linalg.norm((I-np.linalg.inv(Qd).dot(Q)) ** k)) for k in range(K)], 'x-.', color='black', label='est')
+    plt.grid()
+    plt.legend()
+    plt.show()
+
+
+def estimates_over_lambda():
+    M = 5
+    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:]
+
+    Qdinv = np.linalg.inv(Qd)
+    I = np.eye(M)
+    # lam = 1/np.linalg.eigvals(Q)[0]
+    # print(np.linalg.inv(I - lam * Q))
+    # print(np.linalg.norm(1/lam * Qdinv, np.inf))
+    # exit()
+
+    # print(np.linalg.norm((I-Qdinv.dot(Q)) ** K))
+    # lam_list = np.linspace(start=20j, stop=0j, num=1000, endpoint=False)
+    lam_list = np.linspace(start=-1000, stop=0, num=100, endpoint=True)
+
+    rho = []
+    est = []
+    infnorm = []
+    twonorm = []
+    for lam in lam_list:
+        # C = I - lam * Q
+        P = I - lam * Qd
+        Pinv = np.linalg.inv(P)
+
+        R = Pinv.dot(lam * (Q - Qd))
+        # w, V = np.linalg.eig(R)
+        # Vinv = np.linalg.inv(V)
+        # assert np.linalg.norm(V.dot(np.diag(w)).dot(Vinv) - R) < 1E-14, np.linalg.norm(V.dot(np.diag(w)).dot(Vinv) - R)
+        rho.append(max(abs(np.linalg.eigvals(R))))
+        # est.append(0.62*(-1/lam) ** (1/(M-1)))  # M = 3
+        # est.append(0.57*(-1/lam) ** (1/(M-1)))  # M = 5
+        # est.append(0.71*(-1/lam) ** (1/(M-1.5)))  # M = 7
+        # est.append(0.92*(-1/lam) ** (1/(M-2.5)))  # M = 9
+        # est.append((-1/lam) ** (1/(M)))
+        est.append(1000*(1/abs(lam)) ** (1/(1)))
+        # est.append(abs(lam))
+        # est.append(np.linalg.norm((I-Qdinv.dot(Q)) ** M) + 1/abs(lam))
+        # est.append(np.linalg.norm(np.linalg.inv(I - lam * Qd), np.inf))
+        # est.append(np.linalg.norm(np.linalg.inv(I - np.linalg.inv(I - lam * np.diag(Qd)).dot(lam*np.tril(Qd, -1))), np.inf))
+
+        infnorm.append(np.linalg.norm(np.linalg.matrix_power(R, M), np.inf))
+        # infnorm.append(np.linalg.norm(Vinv.dot(R).dot(V), np.inf))
+        twonorm.append(np.linalg.norm(np.linalg.matrix_power(R, M), 2))
+        # twonorm.append(np.linalg.norm(Vinv.dot(R).dot(V), 2))
+    plt.figure()
+    plt.semilogy(lam_list, rho, '--', color='green', label='spectral')
+    plt.semilogy(lam_list, est, '--', color='red', label='est')
+    # plt.semilogy(lam_list, infnorm, 'x--', color='blue', label='infnorm')
+    # plt.semilogy(lam_list, twonorm, 'd--', color='cyan', label='twonorm')
+    # plt.semilogy(np.imag(lam_list), rho, '--', color='green', label='spectral')
+    # plt.semilogy(np.imag(lam_list), est, '--', color='red', label='est')
+    # plt.semilogy(np.imag(lam_list), infnorm, '--', color='blue', label='infnorm')
+    # plt.semilogy(np.imag(lam_list), twonorm, '--', color='cyan', label='twonorm')
+    plt.grid()
+    plt.legend()
+    plt.show()
+
+if __name__ == '__main__':
+    iteration_vs_estimate()
+    # estimates_over_lambda()
diff --git a/pySDC/playgrounds/ODEs/ODE1_playground.py b/pySDC/playgrounds/ODEs/ODE1_playground.py
@@ -0,0 +1,66 @@
+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.nonlinear_ODE_1 import nonlinear_ODE_1
+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'] = 0.01
+
+    # 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'] = 5E-11
+    problem_params['newton_maxiter'] = 500
+    problem_params['u0'] = 0.0
+
+    # 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'] = nonlinear_ODE_1
+    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 = 5
+
+    # 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()
