beautified problem_classes

Author: pancetta

projects/parallelSDC/nonlinear_playground.py

@@ -33,7 +33,7 @@ def main():
     problem_params['nu'] = 1
     problem_params['nvars'] = 255
     problem_params['lambda0'] = 5.0
-    problem_params['maxiter'] = 50
+    problem_params['newton_maxiter'] = 50
     problem_params['newton_tol'] = 1E-12
     problem_params['interval'] = (-5, 5)
 

projects/parallelSDC/preconditioner_playground.py

@@ -119,7 +119,7 @@ def main():
 
                     problem_params['nu'] = 1
                     problem_params['lambda0'] = param
-                    problem_params['maxiter'] = 20
+                    problem_params['newton_maxiter'] = 20
                     problem_params['newton_tol'] = 1E-09
                     problem_params['interval'] = (-5, 5)
 
@@ -215,5 +215,5 @@ def plot_iterations():
 
 
 if __name__ == "__main__":
-    # main()
+    main()
     plot_iterations()

pySDC/core/Errors.py

@@ -46,3 +46,10 @@ class ControllerError(Exception):
     Error class handling/indicating problems with the controller
     """
     pass
+
+
+class ProblemError(Exception):
+    """
+    Error class handling/indicating problems with the problem classes
+    """
+    pass

pySDC/implementations/problem_classes/AdvectionEquation_1D_FD.py

@@ -1,15 +1,17 @@
 from __future__ import division
 import numpy as np
 import scipy.sparse as sp
-import scipy.sparse.linalg as sLA
-import scipy.linalg as LA
+from scipy.sparse.linalg import splu
 
 from pySDC.core.Problem import ptype
+from pySDC.core.Errors import ParameterError, ProblemError
 
 
+# noinspection PyUnusedLocal
 class advection1d(ptype):
     """
-    Example implementing the unforced 1D advection equation with periodic BC in [0,1], discretized using upwinding finite differences
+    Example implementing the unforced 1D advection equation with periodic BC in [0,1],
+    discretized using upwinding finite differences
 
     Attributes:
         A: FD discretization of the gradient operator using upwinding
@@ -21,29 +23,33 @@ def __init__(self, problem_params, dtype_u, dtype_f):
         Initialization routine
 
         Args:
-            problem_params: custom parameters for the example
+            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
-        assert 'nvars' in problem_params, 'ERROR: need number of nvars for the problem class'
-        assert 'c' in problem_params, 'ERROR: need advection parameter for the problem class'
-        assert 'freq' in problem_params, 'ERROR: need frequency parameter for the problem class'
+        essential_keys = ['nvars', 'c', 'freq']
+        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
-        assert (problem_params['nvars'])%2 == 0, 'ERROR: the setup requires nvars = 2^p'
+        if (problem_params['nvars']) % 2 != 0:
+            raise ProblemError('setup requires nvars = 2^p')
 
-        if not 'order' in problem_params:
+        if 'order' not in problem_params:
             problem_params['order'] = 1
-        if not 'type' in problem_params:
+        if 'type' not in problem_params:
             problem_params['type'] = 'upwind'
 
         # invoke super init, passing number of dofs, dtype_u and dtype_f
-        super(advection1d,self).__init__(init=problem_params['nvars'], dtype_u=dtype_u, dtype_f=dtype_f, params=problem_params)
+        super(advection1d, self).__init__(init=problem_params['nvars'], dtype_u=dtype_u, dtype_f=dtype_f,
+                                          params=problem_params)
 
         # compute dx and get discretization matrix A
-        self.dx = 1.0/self.params.nvars
+        self.dx = 1.0 / self.params.nvars
         self.A = self.__get_A(self.params.nvars, self.params.c, self.dx, self.params.order, self.params.type)
 
     @staticmethod
@@ -52,14 +58,14 @@ def __get_A(N, c, dx, order, type):
         Helper function to assemble FD matrix A in sparse format
 
         Args:
-            N: number of dofs
-            c: diffusion coefficient
-            dx: distance between two spatial nodes
-            order: specifies order of discretization
-            type: upwind or centered differences
+            N (int): number of dofs
+            c (float): diffusion coefficient
+            dx (float): distance between two spatial nodes
+            order (int): specifies order of discretization
+            type (string): upwind or centered differences
 
         Returns:
-         matrix A in CSC format
+            scipy.sparse.csc_matrix: matrix A in CSC format
         """
 
         coeff = None
@@ -81,8 +87,7 @@ def __get_A(N, c, dx, order, type):
                 zero_pos = 4
                 coeff = 1.0 / 60.0
             else:
-                print("Order " + order + " not implemented.")
-                exit()
+                raise ProblemError("Order " + str(order) + " not implemented.")
 
         else:
 
@@ -111,66 +116,65 @@ def __get_A(N, c, dx, order, type):
                 coeff = 1.0 / 60.0
                 zero_pos = 5
             else:
-                print("Order " + order + " not implemented.")
-                exit()
+                raise ProblemError("Order " + str(order) + " not implemented.")
 
-        dstencil = np.concatenate((stencil, np.delete(stencil,zero_pos-1)))
-        offsets = np.concatenate(([N-i-1 for i in reversed(range(zero_pos-1))],[i-zero_pos+1 for i in range(zero_pos-1,len(stencil))]))
-        doffsets = np.concatenate((offsets, np.delete(offsets,zero_pos-1)-N))
+        dstencil = np.concatenate((stencil, np.delete(stencil, zero_pos - 1)))
+        offsets = np.concatenate(([N - i - 1 for i in reversed(range(zero_pos - 1))],
+                                  [i - zero_pos + 1 for i in range(zero_pos - 1, len(stencil))]))
+        doffsets = np.concatenate((offsets, np.delete(offsets, zero_pos - 1) - N))
 
-        A = c * coeff * (1.0/dx) * sp.diags(dstencil,doffsets,shape=(N,N),format='csc')
+        A = sp.diags(dstencil, doffsets, shape=(N, N), format='csc')
+        A *= c * coeff * (1.0 / dx)
 
         return A
 
-
-    def eval_f(self,u,t):
+    def eval_f(self, u, t):
         """
         Routine to evaluate the RHS
 
         Args:
-            u: current values
-            t: current time
+            u (dtype_u): current values
+            t (float): current time
 
         Returns:
-            the RHS
+            dtype_f: the RHS
         """
 
         f = self.dtype_f(self.init)
-        f.values = -self.A.dot(u.values)
+        f.values = -1.0 * self.A.dot(u.values)
         return f
 
-    def solve_system(self,rhs,factor,u0,t):
+    def solve_system(self, rhs, factor, u0, t):
         """
-        Simple linear solver for (I-factor*A)u = rhs
+        Simple linear solver for (I+factor*A)u = rhs
 
         Args:
-            rhs: right-hand side for the linear system
-            factor: abbrev. for the node-to-node stepsize (or any other factor required)
-            u0: initial guess for the iterative solver (not used here so far)
-            t: current time (e.g. for time-dependent BCs)
+            rhs (dtype_f): right-hand side for the linear system
+            factor (float) : abbrev. for the node-to-node stepsize (or any other factor required)
+            u0 (dtype_u): initial guess for the iterative solver (not used here so far)
+            t (float): current time (e.g. for time-dependent BCs)
 
         Returns:
-            solution as mesh
+            dtype_u: solution as mesh
         """
 
         me = self.dtype_u(self.init)
-        me.values = sLA.spsolve(sp.eye(self.params.nvars)+factor*self.A,rhs.values)
+        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):
+    def u_exact(self, t):
         """
         Routine to compute the exact solution at time t
 
         Args:
-            t: current time
+            t (float): current time
 
         Returns:
-            exact solution
+            dtype_u: exact solution
         """
 
         me = self.dtype_u(self.init)
-        xvalues = np.array([i*self.dx for i in range(self.params.nvars)])
-        me.values = np.sin(np.pi*self.params.freq*(xvalues - self.params.c*t))
+        xvalues = np.array([i * self.dx for i in range(self.params.nvars)])
+        me.values = np.sin(np.pi * self.params.freq * (xvalues - self.params.c * t))
         return me
-
-

pySDC/implementations/problem_classes/Auzinger_implicit.py

@@ -3,41 +3,43 @@
 import numpy as np
 
 from pySDC.core.Problem import ptype
+from pySDC.core.Errors import ParameterError
 
 
+# noinspection PyUnusedLocal
 class auzinger(ptype):
     """
     Example implementing the Auzinger initial value problem
     """
 
-    def __init__(self, cparams, dtype_u, dtype_f):
+    def __init__(self, problem_params, dtype_u, dtype_f):
         """
         Initialization routine
 
         Args:
-            cparams: custom parameters for the example
+            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
-        assert 'newton_maxiter' in cparams
-        assert 'newton_tol' in cparams
+        essential_keys = ['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)
 
-        # add parameters as attributes for further reference
-        for k, v in cparams.items():
-            setattr(self, k, v)
         # invoke super init, passing dtype_u and dtype_f, plus setting number of elements to 2
-        super(auzinger, self).__init__(2, dtype_u, dtype_f, cparams)
+        super(auzinger, self).__init__(2, dtype_u, dtype_f, problem_params)
 
     def u_exact(self, t):
         """
         Routine for the exact solution
 
         Args:
-            t: current time
+            t (float): current time
         Returns:
-            mesh type containing the exact solution
+            dtype_u: mesh type containing the exact solution
         """
 
         me = self.dtype_u(self.init)
@@ -50,8 +52,8 @@ def eval_f(self, u, t):
         Routine to compute the RHS for both components simultaneously
 
         Args:
-            t: current time (not used here)
-            u: the current values
+            u (dtype_u): the current values
+            t (float): current time (not used here)
         Returns:
             RHS, 2 components
         """
@@ -68,13 +70,13 @@ def solve_system(self, rhs, dt, u0, t):
         Simple Newton solver for the nonlinear system
 
         Args:
-            rhs: right-hand side for the nonlinear system
-            dt: abbrev. for the node-to-node stepsize (or any other factor required)
-            u0: initial guess for the iterative solver
-            t: current time (e.g. for time-dependent BCs)
+            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:
-            solution u
+            dtype_u: solution u
         """
 
         # create new mesh object from u0 and set initial values for iteration
@@ -123,19 +125,6 @@ def solve_system(self, rhs, dt, u0, t):
         #
         #
         # def solve_system_jacobian(self, dfdu, rhs, factor, u0, t):
-        #     """
-        #     Simple linear solver for (I-dtA)u = rhs
-        #
-        #     Args:
-        #         dfdu: the Jacobian of the RHS of the ODE
-        #         rhs: right-hand side for the linear system
-        #         factor: abbrev. for the node-to-node stepsize (or any other factor required)
-        #         u0: initial guess for the iterative solver (not used here so far)
-        #         t: current time (e.g. for time-dependent BCs)
-        #
-        #     Returns:
-        #         solution as mesh
-        #     """
         #
         #     me = mesh(2)
         #     me.values = LA.spsolve(sp.eye(2) - factor * dfdu, rhs.values)