more cleanup

Author: pancetta

pySDC/implementations/problem_classes/HenonHeiles.py

@@ -50,12 +50,12 @@ def eval_f(self, u, t):
 
     def u_exact(self, t):
         """
-        Routine to compute the exact trajectory at time t
+        Routine to compute the exact/initial trajectory at time t
 
         Args:
             t (float): current time
         Returns:
-            dtype_u: exact position and velocity
+            dtype_u: exact/initial position and velocity
         """
         assert t == 0.0, 'error, u_exact only works for the initial time t0=0'
         me = particles(2)

pySDC/implementations/problem_classes/OuterSolarSystem.py

@@ -3,8 +3,8 @@
 import numpy as np
 
 from pySDC.core.Problem import ptype
-from pySDC.implementations.datatype_classes.particles import particles, fields, acceleration
-from pySDC.core.Errors import ParameterError, ProblemError
+from pySDC.implementations.datatype_classes.particles import particles, acceleration
+from pySDC.core.Errors import ParameterError
 
 
 # noinspection PyUnusedLocal

pySDC/implementations/sweeper_classes/verlet.py

@@ -12,10 +12,10 @@ class verlet(sweeper):
     Second-order sweeper using velocity-Verlet as base integrator
 
     Attributes:
-        S: node-to-node collocation matrix (first order)
-        SQ: node-to-node collocation matrix (second order)
-        ST: node-to-node trapezoidal matrix
-        Sx: node-to-node Euler half-step for position update
+        QQ: 0-to-node collocation matrix (second order)
+        QT: 0-to-node trapezoidal matrix
+        Qx: 0-to-node Euler half-step for position update
+        qQ: update rule for final value (if needed)
     """
 
     def __init__(self, params):
@@ -26,17 +26,15 @@ def __init__(self, params):
             params: parameters for the sweeper
         """
 
-        # call parent's initialization routine
-
         if 'QI' not in params:
             params['QI'] = 'IE'
         if 'QE' not in params:
             params['QE'] = 'EE'
 
+        # call parent's initialization routine
         super(verlet, self).__init__(params)
 
-        # S- and SQ-matrices (derived from Q) and Sx- and ST-matrices for the integrator
-        # [self.S, self.ST, self.SQ, self.Sx, self.QQ] = self.__get_Qd()
+        # Trapezoidal rule, Qx and Double-Q as in the Boris-paper
         [self.QT, self.Qx, self.QQ] = self.__get_Qd()
 
         self.qQ = np.dot(self.coll.weights, self.coll.Qmat[1:, 1:])
@@ -55,8 +53,6 @@ def __get_Qd(self):
         # set implicit and explicit Euler matrices
         QI = self.get_Qdelta_implicit(self.coll, self.params.QI)
         QE = self.get_Qdelta_explicit(self.coll, self.params.QE)
-        # QI = np.zeros(self.coll.Qmat.shape)
-        # QE = np.zeros(self.coll.Qmat.shape)
 
         # trapezoidal rule
         QT = 0.5 * (QI + QE)
@@ -66,6 +62,8 @@ def __get_Qd(self):
 
         QQ = np.zeros(np.shape(self.coll.Qmat))
 
+        # if we have Gauss-Lobatto nodes, we can do a magic trick from the Book
+        # this takes Gauss-Lobatto IIIB and create IIIA out of this
         if isinstance(self.coll, CollGaussLobatto):
 
             for m in range(self.coll.num_nodes):
@@ -74,6 +72,7 @@ def __get_Qd(self):
                                                                self.coll.weights[m])
             QQ = np.dot(self.coll.Qmat, QQ)
 
+        # if we do not have Gauss-Lobatto, just multiply Q (will not get a symplectic method, they say)
         else:
 
             QQ = np.dot(self.coll.Qmat, self.coll.Qmat)
@@ -154,6 +153,7 @@ def integrate(self):
             for j in range(1, self.coll.num_nodes + 1):
                 p[-1].pos += L.dt * (L.dt * self.QQ[m, j] * L.f[j]) + L.dt * self.coll.Qmat[m, j] * L.u[0].vel
                 p[-1].vel += L.dt * self.coll.Qmat[m, j] * L.f[j]
+                # we need to set mass and charge here, too, since the code uses the integral to create new particles
                 p[-1].m = L.u[0].m
                 p[-1].q = L.u[0].q
 
@@ -182,6 +182,7 @@ def compute_end_point(self):
             for m in range(self.coll.num_nodes):
                 L.uend.pos += L.dt * (L.dt * self.qQ[m] * L.f[m + 1]) + L.dt * self.coll.weights[m] * L.u[0].vel
                 L.uend.vel += L.dt * self.coll.weights[m] * L.f[m + 1]
+                # remember to set mass and charge here, too
                 L.uend.m = L.u[0].m
                 L.uend.q = L.u[0].q
             # add up tau correction of the full interval (last entry)

pySDC/projects/Hamiltonian/simple_problems.py

@@ -2,6 +2,7 @@
 from collections import defaultdict
 import os
 import dill
+import numpy as np
 
 import pySDC.helpers.plot_helper as plt_helper
 
@@ -13,7 +14,7 @@
 from pySDC.implementations.problem_classes.HarmonicOscillator import harmonic_oscillator
 from pySDC.implementations.transfer_classes.TransferParticles_NoCoarse import particles_to_particles
 
-from pySDC.helpers.stats_helper import filter_stats
+from pySDC.helpers.stats_helper import filter_stats, sort_stats
 from pySDC.projects.Hamiltonian.hamiltonian_output import hamiltonian_output
 
 
@@ -125,18 +126,21 @@ def run_simulation(prob=None):
 
     """
 
+    # check what problem type we have and set up corresponding description and variables
     if prob == 'harmonic':
         description, controller_params = setup_harmonic()
         # set time parameters
         t0 = 0.0
         Tend = 50.0
         num_procs = 100
+        maxmeaniter = 5.6
     elif prob == 'henonheiles':
         description, controller_params = setup_henonheiles()
         # set time parameters
         t0 = 0.0
         Tend = 25.0
         num_procs = 100
+        maxmeaniter = 3.9
     else:
         raise NotImplemented('Problem type not implemented, got %s' % prob)
 
@@ -151,6 +155,30 @@ def run_simulation(prob=None):
     # call main function to get things done...
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
+    # filter statistics by type (number of iterations)
+    filtered_stats = filter_stats(stats, type='niter')
+
+    # convert filtered statistics to list of iterations count, sorted by process
+    iter_counts = sort_stats(filtered_stats, sortby='time')
+
+    # compute and print statistics
+    for item in iter_counts:
+        out = 'Number of iterations for time %4.2f: %2i' % item
+        print(out)
+
+    niters = np.array([item[1] for item in iter_counts])
+    out = '   Mean number of iterations: %4.2f' % np.mean(niters)
+    print(out)
+    out = '   Range of values for number of iterations: %2i ' % np.ptp(niters)
+    print(out)
+    out = '   Position of max/min number of iterations: %2i -- %2i' % \
+          (int(np.argmax(niters)), int(np.argmin(niters)))
+    print(out)
+    out = '   Std and var for number of iterations: %4.2f -- %4.2f' % (float(np.std(niters)), float(np.var(niters)))
+    print(out)
+
+    assert np.mean(niters) <= maxmeaniter, 'Mean number of iterations is too high, got %s' % np.mean(niters)
+
     fname = 'data/' + prob + '.dat'
     f = open(fname, 'wb')
     dill.dump(stats, f)
@@ -167,22 +195,25 @@ def show_results(prob=None, cwd=''):
         prob (str): name of the problem
         cwd (str): current working directory
     """
+
+    # read in the dill data
     f = open(cwd + 'data/' + prob + '.dat', 'rb')
     stats = dill.load(f)
     f.close()
 
+    # extract error in hamiltonian and prepare for plotting
     extract_stats = filter_stats(stats, type='err_hamiltonian')
     result = defaultdict(list)
     for k, v in extract_stats.items():
         result[k.iter].append((k.time, v))
     for k, v in result.items():
-        assert k <= 6, 'Number of iterations is too high for %s, got %s' % (prob, k)
         result[k] = sorted(result[k], key=lambda x: x[0])
 
     plt_helper.mpl.style.use('classic')
     plt_helper.setup_mpl()
     plt_helper.newfig(textwidth=238.96, scale=0.89)
 
+    # Rearrange data for easy plotting
     err_ham = 1
     for k, v in result.items():
         time = [item[0] for item in v]

pySDC/projects/Hamiltonian/solar_system.py

@@ -135,10 +135,13 @@ def show_results(prob=None, cwd=''):
         prob (str): name of the problem
         cwd (str): current working directory
     """
+
+    # read in the dill data
     f = open(cwd + 'data/' + prob + '.dat', 'rb')
     stats = dill.load(f)
     f.close()
 
+    # extract error in hamiltonian and prepare for plotting
     extract_stats = filter_stats(stats, type='err_hamiltonian')
     result = defaultdict(list)
     for k, v in extract_stats.items():
@@ -150,6 +153,7 @@ def show_results(prob=None, cwd=''):
     plt_helper.setup_mpl()
     plt_helper.newfig(textwidth=238.96, scale=0.89)
 
+    # Rearrange data for easy plotting
     err_ham = 1
     for k, v in result.items():
         time = [item[0] for item in v]