New hook, used in boussinesq_2d_imex example to set the GMRES tolerance

Author: pancetta

examples/boussinesq_2d_imex/HookClass.py

@@ -24,7 +24,25 @@ def __init__(self):
         #self.counter = 0
 
     def dump_sweep(self,status):
-        pass
+        """
+        Set new GMRES tolerance depending on the previous SDC residual
+
+        Args:
+            status: status object per step
+        """
+        super(plot_solution,self).dump_sweep(status)
+        self.level.prob.gmres_tol = max(self.level.status.residual*self.level.prob.gmres_tol_factor,self.level.prob.gmres_tol_limit)
+
+    def dump_pre_iteration(self,status):
+        """
+        Set new GMRES tolerance depending on the initial SDC residual
+
+        Args:
+            status: status object per step
+        """
+        super(plot_solution,self).dump_pre_iteration(status)
+        self.level.sweep.compute_residual()
+        self.level.prob.gmres_tol = max(self.level.status.residual*self.level.prob.gmres_tol_factor,self.level.prob.gmres_tol_limit)
 
     def dump_step(self,status):
         """

examples/boussinesq_2d_imex/ProblemClass.py

@@ -66,7 +66,7 @@ def __init__(self, cparams, dtype_u, dtype_f):
         assert 'order' in cparams
         assert 'gmres_maxiter' in cparams
         assert 'gmres_restart' in cparams
-        assert 'gmres_tol' in cparams
+        assert 'gmres_tol_limit' in cparams
 
         # add parameters as attributes for further reference
         for k,v in cparams.items():
@@ -86,6 +86,7 @@ def __init__(self, cparams, dtype_u, dtype_f):
         self.D_upwind   = getBoussinesq2DUpwindMatrix( self.N, self.h[0], self.u_adv , self.order_upw)
 
         self.logger = logging()
+        self.gmres_tol = None
 
     def solve_system(self,rhs,factor,u0,t):
         """
@@ -103,6 +104,8 @@ def solve_system(self,rhs,factor,u0,t):
 
         b         = rhs.values.flatten()
         cb        = Callback()
+        print(self.gmres_tol)
+
         sol, info = LA.gmres( self.Id - factor*self.M, b, x0=u0.values.flatten(), tol=self.gmres_tol, restart=self.gmres_restart, maxiter=self.gmres_maxiter, callback=cb)
         # If this is a dummy call with factor==0.0, do not log because it should not be counted as a solver call
         if factor!=0.0:

examples/boussinesq_2d_imex/plotgmrescounter.py

@@ -15,9 +15,9 @@
   uimex = np.load('rkimex.npy')
   uref  = np.load('uref.npy')
 
-  print "Estimated discretisation error of DIRK: %5.3e" % ( np.linalg.norm(udirk.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) )
-  print "Estimated discretisation error of SDC:  %5.3e" % ( np.linalg.norm(uend.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) )
-  print "Estimated discretisation error of RK-IMEX:  %5.3e" % ( np.linalg.norm(uimex.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) )
+  print("Estimated discretisation error of DIRK: %5.3e" % ( np.linalg.norm(udirk.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) ))
+  print("Estimated discretisation error of SDC:  %5.3e" % ( np.linalg.norm(uend.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) ))
+  print("Estimated discretisation error of RK-IMEX:  %5.3e" % ( np.linalg.norm(uimex.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) ))
 
   fs = 8
   rcParams['figure.figsize'] = 5.0, 2.5

examples/boussinesq_2d_imex/rungmrescounter.py

@@ -37,7 +37,7 @@
     swparams = {}
     swparams['collocation_class'] = collclass.CollGaussLegendre
     swparams['num_nodes'] = 3
-    swparams['do_LU'] = False
+    swparams['do_LU'] = True
 
     sparams = {}
     sparams['maxiter'] = 4
@@ -47,7 +47,7 @@
     # setup parameters "in time"
     t0     = 0
     Tend   = 3000
-    Nsteps =  500
+    Nsteps =  100
     dt = Tend/float(Nsteps)
 
     # This comes as read-in for the problem class
@@ -63,7 +63,8 @@
     pparams['order_upw'] = [5]
     pparams['gmres_maxiter'] = [500]
     pparams['gmres_restart'] = [10]
-    pparams['gmres_tol']     = [1e-9]
+    pparams['gmres_tol_limit'] = [1e-3]
+    pparams['gmres_tol_factor'] = [0.1]
 
     # This comes as read-in for the transfer operations
     tparams = {}
@@ -90,11 +91,11 @@
     cfl_advection    = pparams['u_adv']*dt/P.h[0]
     cfl_acoustic_hor = pparams['c_s']*dt/P.h[0]
     cfl_acoustic_ver = pparams['c_s']*dt/P.h[1]
-    print "Horizontal resolution: %4.2f" % P.h[0]
-    print "Vertical resolution:   %4.2f" % P.h[1]
-    print ("CFL number of advection: %4.2f" % cfl_advection)
-    print ("CFL number of acoustics (horizontal): %4.2f" % cfl_acoustic_hor)
-    print ("CFL number of acoustics (vertical):   %4.2f" % cfl_acoustic_ver)
+    print("Horizontal resolution: %4.2f" % P.h[0])
+    print("Vertical resolution:   %4.2f" % P.h[1])
+    print("CFL number of advection: %4.2f" % cfl_advection)
+    print("CFL number of acoustics (horizontal): %4.2f" % cfl_acoustic_hor)
+    print("CFL number of acoustics (vertical):   %4.2f" % cfl_acoustic_ver)
 
     dirkp = dirk(P, dirk_order)
     u0 = uinit.values.flatten()
@@ -129,18 +130,18 @@
     np.save('rkimex', uimex)
     np.save('uref', uref)
 
-    print " #### Logging report for DIRK-%1i #### " % dirkp.order
-    print "Number of calls to implicit solver: %5i" % dirkp.logger.solver_calls
-    print "Total number of GMRES iterations: %5i" % dirkp.logger.iterations
-    print "Average number of iterations per call: %6.3f" % (float(dirkp.logger.iterations)/float(dirkp.logger.solver_calls))
-    print " "
-    print " #### Logging report for RK-IMEX-%1i #### " % rkimex.order
-    print "Number of calls to implicit solver: %5i" % rkimex.logger.solver_calls
-    print "Total number of GMRES iterations: %5i" % rkimex.logger.iterations
-    print "Average number of iterations per call: %6.3f" % (float(rkimex.logger.iterations)/float(rkimex.logger.solver_calls))
-    print " "
-    print " #### Logging report for SDC-(%1i,%1i) #### " % (swparams['num_nodes'], sparams['maxiter'])
-    print "Number of calls to implicit solver: %5i" % P.logger.solver_calls
-    print "Total number of GMRES iterations: %5i" % P.logger.iterations
-    print "Average number of iterations per call: %6.3f" % (float(P.logger.iterations)/float(P.logger.solver_calls))  
+    print(" #### Logging report for DIRK-%1i #### " % dirkp.order)
+    print("Number of calls to implicit solver: %5i" % dirkp.logger.solver_calls)
+    print("Total number of GMRES iterations: %5i" % dirkp.logger.iterations)
+    print("Average number of iterations per call: %6.3f" % (float(dirkp.logger.iterations)/float(dirkp.logger.solver_calls)))
+    print(" ")
+    print(" #### Logging report for RK-IMEX-%1i #### " % rkimex.order)
+    print("Number of calls to implicit solver: %5i" % rkimex.logger.solver_calls)
+    print("Total number of GMRES iterations: %5i" % rkimex.logger.iterations)
+    print("Average number of iterations per call: %6.3f" % (float(rkimex.logger.iterations)/float(rkimex.logger.solver_calls)))
+    print(" ")
+    print(" #### Logging report for SDC-(%1i,%1i) #### " % (swparams['num_nodes'], sparams['maxiter']))
+    print("Number of calls to implicit solver: %5i" % P.logger.solver_calls)
+    print("Total number of GMRES iterations: %5i" % P.logger.iterations)
+    print("Average number of iterations per call: %6.3f" % (float(P.logger.iterations)/float(P.logger.solver_calls)))
 

examples/boussinesq_2d_imex/standard_integrators.py

@@ -93,7 +93,7 @@ def f_fast(self, u):
 
   def f_fast_solve(self, rhs, alpha, u0):
     cb = Callback()
-    sol, info = LA.gmres( self.problem.Id - alpha*self.problem.M, rhs, x0=u0, tol=self.problem.gmres_tol, restart=self.problem.gmres_restart, maxiter=self.problem.gmres_maxiter, callback=cb)
+    sol, info = LA.gmres( self.problem.Id - alpha*self.problem.M, rhs, x0=u0, tol=self.problem.gmres_tol_limit, restart=self.problem.gmres_restart, maxiter=self.problem.gmres_maxiter, callback=cb)
     if alpha!=0.0:
       #print "RK-IMEX-%1i: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( self.order, cb.getcounter(), cb.getresidual() )
       self.logger.add(cb.getcounter())    
@@ -129,7 +129,7 @@ def f(self,u):
   #  
   def f_solve(self, b, alpha, u0):
     cb = Callback()
-    sol, info = LA.gmres( self.problem.Id - alpha*(self.problem.D_upwind + self.problem.M), b, x0=u0, tol=self.problem.gmres_tol, restart=self.problem.gmres_restart, maxiter=self.problem.gmres_maxiter, callback=cb)
+    sol, info = LA.gmres( self.problem.Id - alpha*(self.problem.D_upwind + self.problem.M), b, x0=u0, tol=self.problem.gmres_tol_limit, restart=self.problem.gmres_restart, maxiter=self.problem.gmres_maxiter, callback=cb)
     if alpha!=0.0:
       #print "BDF-2: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( cb.getcounter(), cb.getresidual() )
       self.logger.add(cb.getcounter())    
@@ -166,7 +166,7 @@ def f(self,u):
   #  
   def f_solve(self, b, alpha, u0):
     cb = Callback()
-    sol, info = LA.gmres( self.problem.Id - alpha*(self.problem.D_upwind + self.problem.M), b, x0=u0, tol=self.problem.gmres_tol, restart=self.problem.gmres_restart, maxiter=self.problem.gmres_maxiter, callback=cb)
+    sol, info = LA.gmres( self.problem.Id - alpha*(self.problem.D_upwind + self.problem.M), b, x0=u0, tol=self.problem.gmres_tol_limit, restart=self.problem.gmres_restart, maxiter=self.problem.gmres_maxiter, callback=cb)
     if alpha!=0.0:
       #print "BDF-2: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( cb.getcounter(), cb.getresidual() )
       self.logger.add(cb.getcounter())    
@@ -283,7 +283,7 @@ def f(self,u):
   #  
   def f_solve(self, b, alpha, u0):
     cb = Callback()
-    sol, info = LA.gmres( self.problem.Id - alpha*(self.problem.D_upwind + self.problem.M), b, x0=u0, tol=self.problem.gmres_tol, restart=self.problem.gmres_restart, maxiter=self.problem.gmres_maxiter, callback=cb)
+    sol, info = LA.gmres( self.problem.Id - alpha*(self.problem.D_upwind + self.problem.M), b, x0=u0, tol=self.problem.gmres_tol_limit, restart=self.problem.gmres_restart, maxiter=self.problem.gmres_maxiter, callback=cb)
     if alpha!=0.0:
       #print "DIRK-%1i: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( self.order, cb.getcounter(), cb.getresidual() )
       self.logger.add(cb.getcounter())    

pySDC/Hooks.py

@@ -53,6 +53,12 @@ def dump_pre(self,status):
         """
         pass
 
+    def dump_pre_iteration(self,status):
+        """
+        Default routine called before iteration starts
+        """
+        pass
+
 
     def dump_sweep(self,status):
         """

pySDC/PFASST_blockwise.py

@@ -251,6 +251,7 @@ def pfasst(MS):
                 if len(S.levels) > 1 and S.params.predict:
                     S.status.stage = 'PREDICT'
                 else:
+                    S.levels[0].hooks.dump_pre_iteration(S.status)
                     S.status.stage = 'IT_FINE'
 
             return MS
@@ -263,6 +264,7 @@ def pfasst(MS):
 
             for S in MS:
                 # update stage
+                S.levels[0].hooks.dump_pre_iteration(S.status)
                 S.status.stage = 'IT_FINE'
 
             return MS

pySDC/PFASST_stepwise.py

@@ -192,6 +192,7 @@ def pfasst(S):
             if len(S.levels) > 1 and S.params.predict:
                 S.status.stage = 'PREDICT_RESTRICT'
             else:
+                S.levels[0].hooks.dump_pre_iteration(S.status)
                 S.status.stage = 'IT_FINE_SWEEP'
             return S
 
@@ -254,6 +255,7 @@ def pfasst(S):
                 S.transfer(source=S.levels[l],target=S.levels[l-1])
 
             # uodate stage and return
+            S.levels[0].hooks.dump_pre_iteration(S.status)
             S.status.stage = 'IT_FINE_SWEEP'
             return S
 

pySDC/Sweeper.py

@@ -95,6 +95,7 @@ def predict(self):
 
         # indicate that this level is now ready for sweeps
         L.status.unlocked = True
+        L.status.updated = True
 
         return None