matrices LHS and RHS in test_imexsweeper are now created by a function of the imex_sweeper class

Author: Daniel Ruprecht

pySDC/sweeper_classes/imex_1st_order.py

@@ -156,12 +156,20 @@ def compute_end_point(self):
         return None
 
     def get_sweeper_mats(self):
+      """
+      Returns the three matrices Q, QI, QE which define the sweeper. The first row and column, corresponding to the left starting value, are removed to correspond to the notation
+      introduced in Ruprecht & Speck, ``Spectral deferred corrections with fast-wave slow-wave splitting'', 2016
+      """
       QE = self.QE[1:,1:]
       QI = self.QI[1:,1:]
       Q  = self.coll.Qmat[1:,1:]
       return QE, QI, Q
 
     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.
+      """
       QE, QI, Q = self.get_sweeper_mats()
       if lambdas==[None,None]:
         pass
@@ -171,4 +179,8 @@ def get_scalar_problems_sweeper_mats(self, lambdas=[None, None]):
         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)"
-      return QE, QI, Q
+      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

tests/test_imexsweeper.py

@@ -32,13 +32,9 @@ def setupLevelStepProblem(self):
     problem = level.prob
     return step, level, problem, nnodes
 
-  def setupSweeperMatrices(self, step, level, problem):
-    nnodes  = step.levels[0].sweep.coll.num_nodes
-    # Build SDC sweep matrix
-    QE, QI, Q = level.sweep.get_sweeper_mats()
-    dt = step.status.dt
-    LHS = np.eye(nnodes) - step.status.dt*( problem.lambda_f[0]*QI + problem.lambda_s[0]*QE )
-    RHS = step.status.dt*( (problem.lambda_f[0]+problem.lambda_s[0])*Q - (problem.lambda_f[0]*QI + problem.lambda_s[0]*QE) )
+  def setupSweeperMatrices(self, level, problem):
+    lambdas = [ problem.lambda_f[0] , problem.lambda_s[0] ]
+    LHS, RHS = level.sweep.get_scalar_problems_sweeper_mats( lambdas = lambdas )
     return LHS, RHS
 
   #
@@ -112,7 +108,7 @@ def test_sweepequalmatrix(self):
       # Perform node-to-node SDC sweep
       level.sweep.update_nodes()
 
-      LHS, RHS = self.setupSweeperMatrices(step, level, problem)
+      LHS, RHS = self.setupSweeperMatrices(level, problem)
 
       unew = np.linalg.inv(LHS).dot( u0full + RHS.dot(u0full) )
       usweep = np.array([ level.u[l].values.flatten() for l in range(1,nnodes+1) ])
@@ -193,7 +189,7 @@ def test_manysweepsequalmatrix(self):
         level.sweep.update_nodes()
       usweep = np.array([ level.u[l].values.flatten() for l in range(1,nnodes+1) ])
 
-      LHS, RHS = self.setupSweeperMatrices(step, level, problem)
+      LHS, RHS = self.setupSweeperMatrices(level, problem)
       unew = u0full
       for i in range(0,K):
         unew = np.linalg.inv(LHS).dot( u0full + RHS.dot(unew) )
@@ -207,6 +203,7 @@ def test_manysweepsequalmatrix(self):
         Mat_sweep = Mat_sweep + np.linalg.matrix_power(Pinv.dot(RHS),i).dot(Pinv)
       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"
+
 
   #
   # Make sure that update function for K sweeps computed from K-sweep matrix gives same result as K sweeps in node-to-node form plus compute_end_point
@@ -226,7 +223,7 @@ def test_manysweepupdate(self):
       level.sweep.compute_end_point()
       uend_sweep = level.uend.values
 
-      LHS, RHS = self.setupSweeperMatrices(step, level, problem)
+      LHS, RHS = self.setupSweeperMatrices(level, problem)
       # Build single matrix representing K sweeps    
       Pinv = np.linalg.inv(LHS)
       Mat_sweep = np.linalg.matrix_power(Pinv.dot(RHS), K)
@@ -238,6 +235,7 @@ def test_manysweepupdate(self):
       uend_matrix = update*self.pparams['u0']
       assert abs(uend_matrix - uend_sweep)<1e-14, "Node-to-node sweep plus update yields different result than update function computed through K-sweep matrix"
 
+
   #
   # Make sure that creating a sweeper object with a collocation object with right_is_node=False and do_coll_update=False throws an exception
   #