adds function get_scalar_problems_manysweep_mat to imex_sweeper which generates a matrix representing K many sweeps of IMEX-SDC; modified test_imexsweeper to test that this formulation is equivalent to K sweeps using update_nodes

Author: Daniel Ruprecht

examples/fwsw/plot_stability.py

@@ -18,8 +18,8 @@
     N_s = 100
     N_f = 400
 
-    lam_s_max = 12.0
-    lam_f_max = 24.0
+    lam_s_max = 3.0
+    lam_f_max = 12.0
     lambda_s = 1j*np.linspace(0.0, lam_s_max, N_s)
     lambda_f = 1j*np.linspace(0.0, lam_f_max, N_f)
 
@@ -31,7 +31,7 @@
     swparams = {}
     swparams['collocation_class'] = collclass.CollGaussLegendre
     swparams['num_nodes'] = 3
-    K = 0
+    K = 3
     do_coll_update = True
 
     #
@@ -47,10 +47,6 @@
     nnodes  = step.levels[0].sweep.coll.num_nodes
     level   = step.levels[0]
     problem = level.prob
-  
-    QE = level.sweep.QE[1:,1:]
-    QI = level.sweep.QI[1:,1:]
-    Q  = level.sweep.coll.Qmat[1:,1:]
 
     stab = np.zeros((N_f, N_s), dtype='complex')
 
@@ -59,15 +55,9 @@
         lambda_fast = lambda_f[j]
         lambda_slow = lambda_s[i]
         if K is not 0:
-
-          LHS = np.eye(nnodes) - step.status.dt*( lambda_fast*QI + lambda_slow*QE )
-          RHS = step.status.dt*( (lambda_fast+lambda_slow)*Q - (lambda_fast*QI + lambda_slow*QE) )
-
-          Pinv = np.linalg.inv(LHS)
-          Mat_sweep = np.linalg.matrix_power(Pinv.dot(RHS), K)
-          for k in range(0,K):
-            Mat_sweep = Mat_sweep + np.linalg.matrix_power(Pinv.dot(RHS),k).dot(Pinv)
-
+          lambdas = [ lambda_fast, lambda_slow ]
+          LHS, RHS = level.sweep.get_scalar_problems_sweeper_mats( lambdas = lambdas )
+          Mat_sweep = level.sweep.get_scalar_problems_manysweep_mat( nsweeps = K, lambdas = lambdas )
         else:
           # Compute stability function of collocation solution
           Mat_sweep = np.linalg.inv(np.eye(nnodes)-step.status.dt*(lambda_fast + lambda_slow)*Q)

pySDC/sweeper_classes/imex_1st_order.py

@@ -169,6 +169,8 @@ def get_scalar_problems_sweeper_mats(self, lambdas=[None, None]):
       """
       For a scalar problem, an IMEX-SDC sweep can be written in matrix formulation. This function returns the corresponding matrices. 
       See Ruprecht & Speck, ``Spectral deferred corrections with fast-wave slow-wave splitting'', 2016 for the derivation.
+
+      The first entry in lambdas is lambda_fast, the second is lambda_slow.
       """
       QE, QI, Q = self.get_sweeper_mats()
       if lambdas==[None,None]:
@@ -178,9 +180,19 @@ def get_scalar_problems_sweeper_mats(self, lambdas=[None, None]):
       else:
         lambda_fast = lambdas[0]
         lambda_slow = lambdas[1]
-      assert abs(lambda_slow)<=abs(lambda_fast), "First entry in parameter lambdas (lambda_fast) has to be greater than second entry (lambda_slow)"
       nnodes = self.coll.num_nodes
       dt = self.level.dt
       LHS = np.eye(nnodes) - dt*( lambda_fast*QI + lambda_slow*QE )
       RHS = dt*( (lambda_fast+lambda_slow)*Q - (lambda_fast*QI + lambda_slow*QE) )
       return LHS, RHS
+
+    def get_scalar_problems_manysweep_mat(self, nsweeps, lambdas=[None,None]):
+      """
+      For a scalar problem, K sweeps of IMEX-SDC can be written in matrix form.
+      """
+      LHS, RHS = self.get_scalar_problems_sweeper_mats(lambdas = lambdas)
+      Pinv = np.linalg.inv(LHS)
+      Mat_sweep = np.linalg.matrix_power(Pinv.dot(RHS), nsweeps)
+      for k in range(0,nsweeps):
+        Mat_sweep += np.linalg.matrix_power(Pinv.dot(RHS), k).dot(Pinv)
+      return Mat_sweep

tests/test_imexsweeper.py

@@ -116,21 +116,22 @@ def test_sweepequalmatrix(self):
   def test_updateformula(self):
     for type in classes:
       self.swparams['collocation_class'] = getattr(pySDC.CollocationClasses, type)
+
       step, level, problem, nnodes = self.setupLevelStepProblem()
       level.sweep.predict()
       u0full = np.array([ level.u[l].values.flatten() for l in range(1,nnodes+1) ])
 
       # Perform update step in sweeper
       level.sweep.update_nodes()
       ustages = np.array([ level.u[l].values.flatten() for l in range(1,nnodes+1) ])
-
       # Compute end value through provided function
       level.sweep.compute_end_point()
       uend_sweep = level.uend.values
       # Compute end value from matrix formulation
       uend_mat   = self.pparams['u0'] + step.status.dt*level.sweep.coll.weights.dot(ustages*(problem.lambda_s[0] + problem.lambda_f[0]))
       assert np.linalg.norm(uend_sweep - uend_mat, np.infty)<1e-14, "Update formula in sweeper gives different result than matrix update formula"
 
+
   #
   # Compute the exact collocation solution by matrix inversion and make sure it is a fixed point
   #
@@ -193,11 +194,7 @@ def test_manysweepsequalmatrix(self):
 
       assert np.linalg.norm(unew - usweep, np.infty)<1e-14, "Doing multiple node-to-node sweeps yields different result than same number of matrix-form sweeps"   
 
-      # Build single matrix representing K sweeps    
-      Pinv = np.linalg.inv(LHS)
-      Mat_sweep = np.linalg.matrix_power(Pinv.dot(RHS), K)
-      for i in range(0,K):
-        Mat_sweep = Mat_sweep + np.linalg.matrix_power(Pinv.dot(RHS),i).dot(Pinv)
+      Mat_sweep = level.sweep.get_scalar_problems_manysweep_mat( nsweeps = K, lambdas = lambdas )
       usweep_onematrix = Mat_sweep.dot(u0full)
       assert np.linalg.norm( usweep_onematrix - usweep, np.infty )<1e-14, "Single-matrix multiple sweep formulation yields different result than multiple sweeps in node-to-node or matrix form form"
 
@@ -208,6 +205,7 @@ def test_manysweepsequalmatrix(self):
   def test_manysweepupdate(self):
     for type in classes:
       self.swparams['collocation_class'] = getattr(pySDC.CollocationClasses, type)
+
       step, level, problem, nnodes = self.setupLevelStepProblem()
       step.levels[0].sweep.predict()
       u0full = np.array([ level.u[l].values.flatten() for l in range(1,nnodes+1) ])
@@ -224,10 +222,7 @@ def test_manysweepupdate(self):
       LHS, RHS = level.sweep.get_scalar_problems_sweeper_mats( lambdas = lambdas )
 
       # Build single matrix representing K sweeps    
-      Pinv = np.linalg.inv(LHS)
-      Mat_sweep = np.linalg.matrix_power(Pinv.dot(RHS), K)
-      for i in range(0,K):
-        Mat_sweep = Mat_sweep + np.linalg.matrix_power(Pinv.dot(RHS),i).dot(Pinv)
+      Mat_sweep = level.sweep.get_scalar_problems_manysweep_mat( nsweeps = K, lambdas = lambdas )
       # Now build update function
       update = 1.0 + (problem.lambda_s[0] + problem.lambda_f[0])*level.sweep.coll.weights.dot(Mat_sweep.dot(np.ones(nnodes)))
       # Multiply u0 by value of update function to get end value directly