add full solar system problem and run

Author: pancetta

docs/source/projects/doc_solar_system.rst

@@ -4,7 +4,16 @@ Full code: `pySDC/projects/Hamiltonian/solar_system.py <https://github.com/Paral
 
 Results:
 
+.. literalinclude:: ../../../outer_solar_system_out.txt
+
 .. image:: ../../../data/outer_solar_system_hamiltonian.png
-   :scale: 50 %
+   :scale: 100 %
 .. image:: ../../../data/outer_solar_system_positions.png
-   :scale: 50 %
+   :scale: 100 %
+
+.. literalinclude:: ../../../full_solar_system_out.txt
+
+.. image:: ../../../data/full_solar_system_hamiltonian.png
+   :scale: 100 %
+.. image:: ../../../data/full_solar_system_positions.png
+   :scale: 100 %

pySDC/implementations/problem_classes/FullSolarSystem.py

@@ -0,0 +1,91 @@
+from __future__ import division
+
+import numpy as np
+
+from pySDC.implementations.problem_classes.OuterSolarSystem import outer_solar_system
+from pySDC.implementations.datatype_classes.particles import particles, acceleration
+from pySDC.core.Errors import ParameterError
+
+
+# noinspection PyUnusedLocal
+class full_solar_system(outer_solar_system):
+    """
+    Example implementing the full solar system problem
+    """
+
+    def __init__(self, problem_params, dtype_u, dtype_f):
+        """
+        Initialization routine
+
+        Args:
+            problem_params (dict): custom parameters for the example
+            dtype_u: particle data type (will be passed parent class)
+            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:
+            if key not in problem_params:
+                msg = 'need %s to instantiate problem, only got %s' % (key, str(problem_params.keys()))
+                raise ParameterError(msg)
+
+        # invoke parant's class (!) super init, passing nparts, dtype_u and dtype_f
+        super(outer_solar_system, self).__init__((3, 10), dtype_u, dtype_f, problem_params)
+
+        # gravitational constant
+        self.G = 2.95912208286E-4
+
+    def u_exact(self, t):
+        """
+        Routine to compute the exact/initial trajectory at time t
+
+        Args:
+            t (float): current time
+        Returns:
+            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)
+
+        # initial positions and velocities taken from
+        # https://www.aanda.org/articles/aa/full/2002/08/aa1405/aa1405.right.html
+        me.pos.values[:, 0] = [0.0, 0.0, 0.0]
+        me.pos.values[:, 1] = [-2.503321047836E-01,	+1.873217481656E-01, +1.260230112145E-01]
+        me.pos.values[:, 2] = [+1.747780055994E-02,	-6.624210296743E-01, -2.991203277122E-01]
+        me.pos.values[:, 3] = [-9.091916173950E-01,	+3.592925969244E-01, +1.557729610506E-01]
+        me.pos.values[:, 4] = [+1.203018828754E+00,	+7.270712989688E-01, +3.009561427569E-01]
+        me.pos.values[:, 5] = [+3.733076999471E+00,	+3.052424824299E+00, +1.217426663570E+00]
+        me.pos.values[:, 6] = [+6.164433062913E+00, +6.366775402981E+00, +2.364531109847E+00]
+        me.pos.values[:, 7] = [+1.457964661868E+01,	-1.236891078519E+01, -5.623617280033E+00]
+        me.pos.values[:, 8] = [+1.695491139909E+01,	-2.288713988623E+01, -9.789921035251E+00]
+        me.pos.values[:, 9] = [-9.707098450131E+00, -2.804098175319E+01, -5.823808919246E+00]
+
+        me.vel.values[:, 0] = [0.0, 0.0, 0.0]
+        me.vel.values[:, 1] = [-2.438808424736E-02, -1.850224608274E-02, -7.353811537540E-03]
+        me.vel.values[:, 2] = [+2.008547034175E-02, +8.365454832702E-04, -8.947888514893E-04]
+        me.vel.values[:, 3] = [-7.085843239142E-03,	-1.455634327653E-02, -6.310912842359E-03]
+        me.vel.values[:, 4] = [-7.124453943885E-03,	+1.166307407692E-02, +5.542098698449E-03]
+        me.vel.values[:, 5] = [-5.086540617947E-03,	+5.493643783389E-03, +2.478685100749E-03]
+        me.vel.values[:, 6] = [-4.426823593779E-03, +3.394060157503E-03, +1.592261423092E-03]
+        me.vel.values[:, 7] = [+2.647505630327E-03, +2.487457379099E-03, +1.052000252243E-03]
+        me.vel.values[:, 8] = [+2.568651772461E-03,	+1.681832388267E-03, +6.245613982833E-04]
+        me.vel.values[:, 9] = [+3.034112963576E-03,	-1.111317562971E-03, -1.261841468083E-03]
+
+        # masses relative to the sun taken from
+        # https://en.wikipedia.org/wiki/Planetary_mass#Values_from_the_DE405_ephemeris
+        me.m[0] = 1.0                   # Sun
+        me.m[1] = 0.1660100 * 1E-06     # Mercury
+        me.m[2] = 2.4478383 * 1E-06     # Venus
+        me.m[3] = 3.0404326 * 1E-06    # Earth+Moon
+        me.m[4] = 0.3227151 * 1E-06    # Mars
+        me.m[5] = 954.79194 * 1E-06    # Jupiter
+        me.m[6] = 285.88600 * 1E-06     # Saturn
+        me.m[7] = 43.662440 * 1E-06    # Uranus
+        me.m[8] = 51.513890 * 1E-06    # Neptune
+        me.m[9] = 0.0073960 * 1E-06     # Pluto
+
+        return me

pySDC/projects/Hamiltonian/solar_system.py

@@ -12,6 +12,7 @@
 from pySDC.implementations.sweeper_classes.verlet import verlet
 from pySDC.implementations.datatype_classes.particles import particles, acceleration
 from pySDC.implementations.problem_classes.OuterSolarSystem import outer_solar_system
+from pySDC.implementations.problem_classes.FullSolarSystem import full_solar_system
 from pySDC.implementations.transfer_classes.TransferParticles_NoCoarse import particles_to_particles
 
 from pySDC.helpers.stats_helper import filter_stats, sort_stats
@@ -67,6 +68,55 @@ def setup_outer_solar_system():
     return description, controller_params
 
 
+def setup_full_solar_system():
+    """
+    Helper routine for setting up everything for the full solar system problem
+
+    Returns:
+        description (dict): description of the controller
+        controller_params (dict): controller parameters
+    """
+
+    # initialize level parameters
+    level_params = dict()
+    level_params['restol'] = 1E-10
+    level_params['dt'] = 10.0
+
+    # initialize sweeper parameters
+    sweeper_params = dict()
+    sweeper_params['collocation_class'] = CollGaussLobatto
+    sweeper_params['num_nodes'] = [5, 3]
+    sweeper_params['spread'] = True
+
+    # initialize problem parameters for the Penning trap
+    problem_params = dict()
+    problem_params['sun_only'] = [False, True]
+
+    # initialize step parameters
+    step_params = dict()
+    step_params['maxiter'] = 50
+
+    # initialize controller parameters
+    controller_params = dict()
+    controller_params['hook_class'] = hamiltonian_output  # specialized hook class for more statistics and output
+    controller_params['logger_level'] = 30
+    controller_params['predict'] = False
+
+    # Fill description dictionary for easy hierarchy creation
+    description = dict()
+    description['problem_class'] = full_solar_system
+    description['problem_params'] = problem_params
+    description['dtype_u'] = particles
+    description['dtype_f'] = acceleration
+    description['sweeper_class'] = verlet
+    description['sweeper_params'] = sweeper_params
+    description['level_params'] = level_params
+    description['step_params'] = step_params
+    description['space_transfer_class'] = particles_to_particles
+
+    return description, controller_params
+
+
 def run_simulation(prob=None):
     """
     Routine to run the simulation of a second order problem
@@ -80,11 +130,24 @@ def run_simulation(prob=None):
         description, controller_params = setup_outer_solar_system()
         # set time parameters
         t0 = 0.0
-        Tend = 20000.0
+        Tend = 10000.0
         num_procs = 100
+        maxmeaniter = 4.0
+    elif prob == 'full_solar_system':
+        description, controller_params = setup_full_solar_system()
+        # set time parameters
+        t0 = 0.0
+        Tend = 1000.0
+        num_procs = 100
+        maxmeaniter = 16.0
     else:
         raise NotImplemented('Problem type not implemented, got %s' % prob)
 
+    f = open(prob + '_out.txt', 'w')
+    out = 'Running ' + prob + ' problem with %s processors...' % num_procs
+    f.write(out + '\n')
+    print(out)
+
     # instantiate the controller
     controller = allinclusive_classic_nonMPI(num_procs=num_procs, controller_params=controller_params,
                                              description=description)
@@ -103,22 +166,28 @@ def run_simulation(prob=None):
     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)
+    # for item in iter_counts:
+    #     out = 'Number of iterations for time %4.2f: %2i' % item
+    #     f.write(out)
+    #     print(out)
 
     niters = np.array([item[1] for item in iter_counts])
     out = '   Mean number of iterations: %4.2f' % np.mean(niters)
+    f.write(out + '\n')
     print(out)
     out = '   Range of values for number of iterations: %2i ' % np.ptp(niters)
+    f.write(out + '\n')
     print(out)
     out = '   Position of max/min number of iterations: %2i -- %2i' % \
           (int(np.argmax(niters)), int(np.argmin(niters)))
+    f.write(out + '\n')
     print(out)
     out = '   Std and var for number of iterations: %4.2f -- %4.2f' % (float(np.std(niters)), float(np.var(niters)))
+    f.write(out + '\n')
     print(out)
+    f.close()
 
-    assert np.mean(niters) <= 4.0, 'Mean number of iterations is too high, got %s' % np.mean(niters)
+    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')
@@ -213,6 +282,9 @@ def main():
     prob = 'outer_solar_system'
     run_simulation(prob)
     show_results(prob)
+    prob = 'full_solar_system'
+    run_simulation(prob)
+    show_results(prob)
 
 
 if __name__ == "__main__":