fixes for solar system problem, WIP

Author: pancetta

pySDC/implementations/sweeper_classes/verlet.py

@@ -145,10 +145,8 @@ def integrate(self):
         # get current level and problem description
         L = self.level
         P = L.prob
-
         # create new instance of dtype_u, initialize values with 0
         p = []
-
         for m in range(1, self.coll.num_nodes + 1):
             p.append(P.dtype_u(P.init, val=0))
 
@@ -184,6 +182,8 @@ 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]
+                L.uend.m = L.u[0].m
+                L.uend.q = L.u[0].q
             # add up tau correction of the full interval (last entry)
             if L.tau is not None:
                 L.uend += L.tau[-1]

pySDC/implementations/transfer_classes/TransferParticles_NoCoarse.py

@@ -1,5 +1,5 @@
 from pySDC.core.SpaceTransfer import space_transfer
-from pySDC.implementations.datatype_classes.particles import particles, fields
+from pySDC.implementations.datatype_classes.particles import particles, fields, acceleration
 from pySDC.core.Errors import TransferError
 
 
@@ -35,6 +35,8 @@ def restrict(self, F):
             G = particles(F)
         elif isinstance(F, fields):
             G = fields(F)
+        elif isinstance(F, acceleration):
+            G = acceleration(F)
         else:
             raise TransferError("Unknown type of fine data, got %s" % type(F))
         return G
@@ -51,6 +53,8 @@ def prolong(self, G):
             F = particles(G)
         elif isinstance(G, fields):
             F = fields(G)
+        elif isinstance(G, acceleration):
+            F = acceleration(G)
         else:
             raise TransferError("Unknown type of coarse data, got %s" % type(G))
         return F

pySDC/projects/Hamiltonian/solar_system.py

@@ -33,8 +33,8 @@ def setup_outer_solar_system():
     # initialize sweeper parameters
     sweeper_params = dict()
     sweeper_params['collocation_class'] = CollGaussLobatto
-    sweeper_params['num_nodes'] = [5]
-    sweeper_params['spread'] = False
+    sweeper_params['num_nodes'] = [5, 3]
+    sweeper_params['spread'] = True
 
     # initialize problem parameters for the Penning trap
     problem_params = dict()
@@ -59,7 +59,7 @@ def setup_outer_solar_system():
     description['sweeper_params'] = sweeper_params
     description['level_params'] = level_params
     description['step_params'] = step_params
-    # description['space_transfer_class'] = particles_to_particles
+    description['space_transfer_class'] = particles_to_particles
 
     return description, controller_params
 
@@ -77,8 +77,8 @@ def run_simulation(prob=None):
         description, controller_params = setup_outer_solar_system()
         # set time parameters
         t0 = 0.0
-        Tend = 100000.0
-        num_procs = 1
+        Tend = 10000.0
+        num_procs = 100
     else:
         raise NotImplemented('Problem type not implemented, got %s' % prob)
 
@@ -118,7 +118,7 @@ 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 <= 6, 'Number of iterations is too high for %s, got %s' % (prob, k)
+        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')
@@ -131,8 +131,7 @@ def show_results(prob=None, cwd=''):
         ham = [item[1] for item in v]
         err_ham = ham[-1]
         plt_helper.plt.semilogy(time, ham, '-', lw=1, label='Iter ' + str(k))
-
-    # assert err_ham < 2.3E-08, 'Error in the Hamiltonian is too large for %s, got %s' % (prob, err_ham)
+    assert err_ham < 3.0E-15, 'Error in the Hamiltonian is too large for %s, got %s' % (prob, err_ham)
 
     plt_helper.plt.xlabel('Time')
     plt_helper.plt.ylabel('Error in Hamiltonian')

pySDC/tests/test_projects/test_hamiltonian/test_solar.py

@@ -0,0 +1,4 @@
+from pySDC.projects.Hamiltonian.solar_system import main
+
+def test_main():
+    main()