[fix] Heun works correctly for p2p.

the current solution is ugly, I do a lot of copying of temporary arrays,
I believe we can do better.
But at least it works

cpp/particles/interpolation/abstract_particle_set.hpp

@@ -97,6 +97,28 @@ class abstract_particle_set
             return EXIT_SUCCESS;
         }
 
+        int copy_state_tofrom(
+                particle_rnumber *dst,
+                const particle_rnumber *src)
+        {
+            return this->LOOP(
+                    [&](const partsize_t idx)
+                    {
+                        dst[idx] = src[idx];
+                    });
+        }
+
+        int copy_state_tofrom(
+                std::unique_ptr<particle_rnumber[]> &dst,
+                const std::unique_ptr<particle_rnumber[]> &src)
+        {
+            return this->LOOP(
+                    [&](const partsize_t idx)
+                    {
+                        dst[idx] = src[idx];
+                    });
+        }
+
         /** Reorganize particles within MPI domain.
          *
          * Based on current particle locations, redistribute the full state

cpp/particles/particle_solver.cpp

@@ -90,25 +90,27 @@ int particle_solver::Heun(
     // so we must put everything in the same vector.
     std::vector<std::unique_ptr<particle_rnumber[]>> tvalues;
     std::unique_ptr<particle_rnumber[]> xx, x0, k1, k2;
-    // allocate k0
-    x0.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
-    // allocate k1
-    k1.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
-    // allocate k2
-    k2.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
-    // allocate xx
-    xx.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
 
+    // allocate x0
+    x0.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
     // store initial condition
     this->pset->getParticleState(x0.get());
 
     tvalues.resize(1);
-    tvalues[0].reset(new particle_rnumber[1]);
-    tvalues[0].swap(x0);
+    tvalues[0].reset(new particle_rnumber[
+            this->pset->getLocalNumberOfParticles()*
+            this->pset->getStateSize()]);
+    this->pset->copy_state_tofrom(tvalues[0], x0);
+    // allocate k1
+    k1.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
     // compute k1
+    // possibly change local ordering, but size of particle state remains the same
     rhs(0, *(this->pset), k1.get(), tvalues);
-    x0.swap(tvalues[0]);
+    // copy back initial date to x0
+    this->pset->copy_state_tofrom(x0, tvalues[0]);
 
+    // allocate xx
+    xx.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
     // compute xn+h*k1 in particle set
     this->pset->LOOP(
             [&](const partsize_t idx){
@@ -120,18 +122,28 @@ int particle_solver::Heun(
     // the subsequent computations are meaningless (and can lead to segfaults).
     this->pset->setParticleState(xx.get());
     tvalues.resize(2);
-    tvalues[0].reset(new particle_rnumber[1]);
-    tvalues[1].reset(new particle_rnumber[1]);
-    tvalues[0].swap(x0);
-    tvalues[1].swap(k1);
+    tvalues[0].reset(new particle_rnumber[
+            this->pset->getLocalNumberOfParticles()*
+            this->pset->getStateSize()]);
+    tvalues[1].reset(new particle_rnumber[
+            this->pset->getLocalNumberOfParticles()*
+            this->pset->getStateSize()]);
+    this->pset->copy_state_tofrom(tvalues[0], x0);
+    this->pset->copy_state_tofrom(tvalues[1], k1);
     this->pset->redistribute(tvalues);
     // impose model constraints, so that rhs is computed correctly
     rhs.imposeModelConstraints(*(this->pset));
+    // allocate k2
+    k2.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
     // compute k2
     // note the temporal interpolation is at `1`, i.e. end of time-step
     rhs(1, *(this->pset), k2.get(), tvalues);
-    x0.swap(tvalues[0]);
-    k1.swap(tvalues[1]);
+    // reallocate k0
+    x0.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
+    // reallocate k1
+    k1.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
+    this->pset->copy_state_tofrom(x0, tvalues[0]);
+    this->pset->copy_state_tofrom(k1, tvalues[1]);
     // tvalues no longer required
     tvalues.resize(0);
 

cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp

@@ -56,20 +56,16 @@ class tracer_with_collision_counter_rhs: public tracer_rhs<rnumber, be, tt>
             additional_data.push_back(std::unique_ptr<abstract_particle_set::particle_rnumber[]>(
                     new abstract_particle_set::particle_rnumber[pset.getLocalNumberOfParticles()*pset.getStateSize()]));
             // copy rhs values to temporary array
-            pset.LOOP(
-                    [&](const abstract_particle_set::partsize_t idx)
-                    {
-                        additional_data[additional_data.size()-1][idx] = result[idx];
-                    });
+            pset.copy_state_tofrom(
+                    additional_data[additional_data.size()-1].get(),
+                    result);
             pset.template applyP2PKernel<3, p2p_ghost_collisions<abstract_particle_set::particle_rnumber, abstract_particle_set::partsize_t>>(
                     this->p2p_gc,
                     additional_data);
             // copy shuffled rhs values
-            pset.LOOP(
-                    [&](const abstract_particle_set::partsize_t idx)
-                    {
-                        result[idx] = additional_data[additional_data.size()-1][idx];
-                    });
+            pset.copy_state_tofrom(
+                    result,
+                    additional_data[additional_data.size()-1].get());
             // clear temporary array
             additional_data.pop_back();
             return interpolation_result;