bugfix, better coarsening, readme

Author: pancetta

pySDC/implementations/problem_classes/OuterSolarSystem.py

@@ -10,7 +10,7 @@
 # noinspection PyUnusedLocal
 class outer_solar_system(ptype):
     """
-    Example implementing the harmonic oscillator
+    Example implementing the outer solar system problem
     """
 
     def __init__(self, problem_params, dtype_u, dtype_f):
@@ -23,6 +23,9 @@ def __init__(self, problem_params, dtype_u, dtype_f):
             dtype_f: acceleration data type (will be passed parent class)
         """
 
+        if 'sun_only' not in problem_params:
+            problem_params['sun_only'] = False
+
         # these parameters will be used later, so assert their existence
         essential_keys = []
         for key in essential_keys:
@@ -32,6 +35,8 @@ def __init__(self, problem_params, dtype_u, dtype_f):
 
         # invoke super init, passing nparts, dtype_u and dtype_f
         super(outer_solar_system, self).__init__((3, 6), dtype_u, dtype_f, problem_params)
+
+        # gravitational constant
         self.G = 2.95912208286E-4
 
     def eval_f(self, u, t):
@@ -46,24 +51,37 @@ def eval_f(self, u, t):
         """
         me = acceleration(self.init, val=0.0)
 
-        for i in range(self.init[-1]):
-            for j in range(i):
-                dx = u.pos.values[:, i] - u.pos.values[:, j]
+        # compute the acceleration due to gravitational forces
+        # ... only with respect to the sun
+        if self.params.sun_only:
+
+            for i in range(1, self.init[-1]):
+                dx = u.pos.values[:, i] - u.pos.values[:, 0]
                 r = np.sqrt(np.dot(dx, dx))
                 df = self.G * dx / (r ** 3)
-                me.values[:, i] -= u.m[j] * df
-                me.values[:, j] += u.m[i] * df
+                me.values[:, i] -= u.m[0] * df
+
+        # ... or with all planets involved
+        else:
+
+            for i in range(self.init[-1]):
+                for j in range(i):
+                    dx = u.pos.values[:, i] - u.pos.values[:, j]
+                    r = np.sqrt(np.dot(dx, dx))
+                    df = self.G * dx / (r ** 3)
+                    me.values[:, i] -= u.m[j] * df
+                    me.values[:, j] += u.m[i] * df
 
         return me
 
     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(self.init)
@@ -82,12 +100,12 @@ def u_exact(self, t):
         me.vel.values[:, 4] = [0.00288930, 0.00114527, 0.00039677]
         me.vel.values[:, 5] = [0.00276725, -0.0017072, -0.00136504]
 
-        me.m[0] = 1.00000597682
-        me.m[1] = 0.000954786104043
-        me.m[2] = 0.000285583733151
-        me.m[3] = 0.0000437273164546
-        me.m[4] = 0.0000517759138449
-        me.m[5] = 1.0 / 130000000.0
+        me.m[0] = 1.00000597682         # Sun
+        me.m[1] = 0.000954786104043     # Jupiter
+        me.m[2] = 0.000285583733151     # Saturn
+        me.m[3] = 0.0000437273164546    # Uranus
+        me.m[4] = 0.0000517759138449    # Neptune
+        me.m[5] = 1.0 / 130000000.0     # Pluto
 
         return me
 

pySDC/projects/Hamiltonian/README.rst

@@ -9,6 +9,7 @@ Simple problems
 
 We first test the integrator for some rather simple problems, namely the harmonic oscillator and the Henon-Heiles problem.
 For both problems we make use of the hook ``hamiltonian_output`` to monitor the deviation from the exact Hamiltonian.
+PFASST is run with 100 processors (virtually parallel) and the deviation from the initial Hamiltonian is shown.
 
 .. include:: doc_hamiltonian_simple.rst
 
@@ -17,7 +18,8 @@ Solar system problem
 
 In this slightly more complex setup we simulate the movement of planets in the outer solar system.
 For this, only six planets are modeled, namely the Sun (which in its mass contains the inner planets), Jupiter, Saturn, Uranus, Neptune and Pluto.
-The gravitational forces are computed using a simple N^2 solver.
+The acceleration due to gravitational forces are computed using a simple N^2 solver.
 All this is implemented within the ``OuterSolarSystem`` problem class.
+Coarsening can be done using only the sun for computing the acceleration.
 
 .. include:: doc_solar_system.rst

pySDC/projects/Hamiltonian/hamiltonian_output.py

@@ -27,7 +27,7 @@ def post_iteration(self, step, level_number):
             step (pySDC.Step.step): the current step
             level_number (int): the current level number
         """
-        super(hamiltonian_output, self).post_step(step, level_number)
+        super(hamiltonian_output, self).post_iteration(step, level_number)
 
         # some abbreviations
         L = step.levels[0]

pySDC/projects/Hamiltonian/solar_system.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
 
@@ -12,7 +13,7 @@
 from pySDC.implementations.problem_classes.OuterSolarSystem import outer_solar_system
 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
 
 
@@ -38,6 +39,7 @@ def setup_outer_solar_system():
 
     # initialize problem parameters for the Penning trap
     problem_params = dict()
+    problem_params['sun_only'] = [False, True]
 
     # initialize step parameters
     step_params = dict()
@@ -93,6 +95,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) <= 4.0, 'Mean number of iterations is too high, got %s' % np.mean(niters)
+
     fname = 'data/' + prob + '.dat'
     f = open(fname, 'wb')
     dill.dump(stats, f)
@@ -118,7 +144,6 @@ def show_results(prob=None, cwd=''):
     for k, v in extract_stats.items():
         result[k.iter].append((k.time, v))
     for k, v in result.items():
-        assert k <= 4, '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')