diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index 79ee01e27c315a7af66c0ab120cc08db3a7b7971..343991283bbdcfa1962b012adc47915de14c0d44 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/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
diff --git a/cpp/particles/particle_solver.cpp b/cpp/particles/particle_solver.cpp
index 65ec05704c8be78f2b4a2de7643a9832ed831638..3783bf1ae9382aa762f7d7892e485e0f810ee81d 100644
--- a/cpp/particles/particle_solver.cpp
+++ b/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);
diff --git a/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp b/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
index 200fcdeb80c5bbe39568264b83b7c2591d4076f5..639d3aa85393d5735c30551764bebb2760957e16 100644
--- a/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
+++ b/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;