Keep spatially extended residual in the levels after computation

Author: brownbaerchen

pySDC/core/level.py

@@ -83,6 +83,7 @@ def __init__(self, problem_class, problem_params, sweeper_class, sweeper_params,
         self.u = [None] * (self.sweep.coll.num_nodes + 1)
         self.uold = [None] * (self.sweep.coll.num_nodes + 1)
         self.u_avg = [None] * self.sweep.coll.num_nodes
+        self.residual = [None] * self.sweep.coll.num_nodes
         self.f = [None] * (self.sweep.coll.num_nodes + 1)
         self.fold = [None] * (self.sweep.coll.num_nodes + 1)
 

pySDC/core/sweeper.py

@@ -192,14 +192,14 @@ def compute_residual(self, stage=''):
 
         # build QF(u)
         res_norm = []
-        res = self.integrate()
+        L.residual = self.integrate()
         for m in range(self.coll.num_nodes):
-            res[m] += L.u[0] - L.u[m + 1]
+            L.residual[m] += L.u[0] - L.u[m + 1]
             # add tau if associated
             if L.tau[m] is not None:
-                res[m] += L.tau[m]
+                L.residual[m] += L.tau[m]
             # use abs function from data type here
-            res_norm.append(abs(res[m]))
+            res_norm.append(abs(L.residual[m]))
 
         # find maximal residual over the nodes
         if L.params.residual_type == 'full_abs':

pySDC/implementations/controller_classes/controller_ParaDiag_nonMPI.py

@@ -152,12 +152,11 @@ def swap_solution_for_all_at_once_residual(self, local_MS_running):
                 hook.post_comm(step=S, level_number=0, add_to_stats=True)
 
             # compute residuals locally
-            residual = S.levels[0].sweep.get_residual()
-            S.levels[0].status.residual = max(abs(me) for me in residual)
+            S.levels[0].sweep.compute_residual()
 
             # put residual in the solution variables
             for m in range(S.levels[0].sweep.coll.num_nodes):
-                S.levels[0].u[m + 1] = residual[m]
+                S.levels[0].u[m + 1] = S.levels[0].residual[m]
 
     def swap_increment_for_solution(self, local_MS_running):
         """

pySDC/implementations/sweeper_classes/ParaDiagSweepers.py

@@ -122,7 +122,6 @@ def update_nodes(self):
         for m in range(M):
             L.u[m + 1] = y[m]
             if self.params.update_f_evals:
-                raise
                 L.f[m + 1] = P.eval_f(L.u[m + 1], L.time + L.dt * self.coll.nodes[m])
 
         L.status.updated = True
@@ -154,6 +153,10 @@ def get_residual(self):
 
         return residual
 
+    def compute_residual(self, *args, **kwargs):
+        self.eval_f_at_all_nodes()
+        return super().compute_residual(*args, **kwargs)
+
 
 class QDiagonalizationIMEX(QDiagonalization):
     """

pySDC/tutorial/step_9/B_paradiag_for_nonlinear_problems.py

@@ -168,7 +168,6 @@ def residual(_u, u0):
 # do ParaDiag
 sol_paradiag = u.copy() * 0j
 u0 = u.copy()
-
 niter = 0
 res = residual(sol_paradiag, u0)
 while np.max(np.abs(res)) > restol: