diff --git a/TurTLE/test/particle_set/ADNS.py b/TurTLE/test/particle_set/ADNS.py
index 1c6a66a9ae2a1628994311cff1256f3096959a01..8e60c82b5cb675301a4b79ef38f1b35ec2c0ceb9 100644
--- a/TurTLE/test/particle_set/ADNS.py
+++ b/TurTLE/test/particle_set/ADNS.py
@@ -81,7 +81,7 @@ class ADNS(TurTLE.DNS):
self.parameters['niter_stat'] = int(20)
self.parameters['niter_out'] = int(200)
self.parameters['checkpoints_per_file'] = int(1)
- self.parameters['dt'] = float(0.001)
+ self.parameters['dt'] = float(0.0001)
self.parameters['histogram_bins'] = int(64)
self.parameters['max_velocity_estimate'] = float(1)
self.parameters['max_vorticity_estimate'] = float(1)
diff --git a/cpp/particles/abstract_particle_rhs.hpp b/cpp/particles/abstract_particle_rhs.hpp
index 3687997809c041c16ca2243f2d9f624c3b79400c..a948a5effd9bf22b1b9534e150dfde01d08fa2e0 100644
--- a/cpp/particles/abstract_particle_rhs.hpp
+++ b/cpp/particles/abstract_particle_rhs.hpp
@@ -38,10 +38,13 @@ class abstract_particle_rhs
virtual ~abstract_particle_rhs(){}
// important bit
+ // may resort particles, therefore it takes the
+ // parameter `additional_data` to possibly resort other arrays
virtual int operator()(
double t,
abstract_particle_set &pset,
- particle_rnumber *result) const = 0;
+ particle_rnumber *result,
+ std::vector<std::unique_ptr<abstract_particle_set::particle_rnumber[]>> &additional_data) const = 0;
// for post-timestep actions such as
// normalization of orientation vector
diff --git a/cpp/particles/particle_solver.cpp b/cpp/particles/particle_solver.cpp
index d817079160a1c369cb819f779d5d869aa21210bf..65ec05704c8be78f2b4a2de7643a9832ed831638 100644
--- a/cpp/particles/particle_solver.cpp
+++ b/cpp/particles/particle_solver.cpp
@@ -34,9 +34,11 @@ int particle_solver::Euler(
auto *rhs_val = new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())];
// temporary array for new particle state
auto *xx = new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())];
+ // empty array required by rhs call
+ std::vector<std::unique_ptr<particle_rnumber[]>> tvalues;
// compute right hand side
- rhs(0, *(this->pset), rhs_val);
+ rhs(0, *(this->pset), rhs_val, tvalues);
// compute new particle state
@@ -50,9 +52,7 @@ int particle_solver::Euler(
//std::vector<std::unique_ptr<particle_rnumber[]>> temporary;
this->pset->setParticleState(xx);
- // fake kvalues for use in redistribute
- std::vector<std::unique_ptr<particle_rnumber[]>> kvalues;
- this->pset->redistribute(kvalues);
+ this->pset->redistribute(tvalues);
// ensure new state is consistent with physical model
rhs.imposeModelConstraints(*(this->pset));
@@ -86,57 +86,65 @@ int particle_solver::Heun(
{
TIMEZONE("particle_solver::Euler");
// temporary arrays for storing the two right-hand-side values
- // I allocate 3 values, because 0 will be used for the initial condition
// We need to be able to redistribute the initial condition for computations to make sense,
// so we must put everything in the same vector.
- std::vector<std::unique_ptr<particle_rnumber[]>> kvalues;
- kvalues.resize(3);
+ std::vector<std::unique_ptr<particle_rnumber[]>> tvalues;
+ std::unique_ptr<particle_rnumber[]> xx, x0, k1, k2;
// allocate k0
- kvalues[0].reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
+ x0.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
// allocate k1
- kvalues[1].reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
+ k1.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
// allocate k2
- kvalues[2].reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
- // temporary array for particle states
- std::unique_ptr<particle_rnumber[]> xx;
+ k2.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
+ // allocate xx
xx.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
// store initial condition
- this->pset->getParticleState(kvalues[0].get());
+ this->pset->getParticleState(x0.get());
+ tvalues.resize(1);
+ tvalues[0].reset(new particle_rnumber[1]);
+ tvalues[0].swap(x0);
// compute k1
- rhs(0, *(this->pset), kvalues[1].get());
+ rhs(0, *(this->pset), k1.get(), tvalues);
+ x0.swap(tvalues[0]);
// compute xn+h*k1 in particle set
this->pset->LOOP(
[&](const partsize_t idx){
- xx[idx] = kvalues[0][idx] + dt*kvalues[1][idx];
+ xx[idx] = x0[idx] + dt*k1[idx];
});
// copy the temporary state to the particle set, and redistribute
// we MUST redistribute, because we want to interpolate.
- // we MUST shuffle kvalues according to the same redistribution, otherwise
+ // we MUST shuffle values according to the same redistribution, otherwise
// the subsequent computations are meaningless (and can lead to segfaults).
this->pset->setParticleState(xx.get());
- this->pset->redistribute(kvalues);
+ 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);
+ this->pset->redistribute(tvalues);
// impose model constraints, so that rhs is computed correctly
rhs.imposeModelConstraints(*(this->pset));
// compute k2
// note the temporal interpolation is at `1`, i.e. end of time-step
- rhs(1, *(this->pset), kvalues[2].get());
+ rhs(1, *(this->pset), k2.get(), tvalues);
+ x0.swap(tvalues[0]);
+ k1.swap(tvalues[1]);
+ // tvalues no longer required
+ tvalues.resize(0);
// local number of particles may have changed
xx.reset(new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())]);
// compute final state
this->pset->LOOP(
[&](const partsize_t idx){
- xx[idx] = kvalues[0][idx] + 0.5*dt*(kvalues[1][idx] + kvalues[2][idx]);
+ xx[idx] = x0[idx] + 0.5*dt*(k1[idx] + k2[idx]);
});
this->pset->setParticleState(xx.get());
- // kvalues no longer required, no point in passing them around with MPI
- kvalues.resize(0);
// redistribute, impose model constraints
- this->pset->redistribute(kvalues);
- DEBUG_MSG("check 03\n");
+ this->pset->redistribute(tvalues);
rhs.imposeModelConstraints(*(this->pset));
// and we are done
diff --git a/cpp/particles/rhs/tracer_rhs.hpp b/cpp/particles/rhs/tracer_rhs.hpp
index 676789ae1ad5ffdcce91d13d2e456fadcce591e0..50da9e90ac7c5f48b2a28b62efd0219de674add6 100644
--- a/cpp/particles/rhs/tracer_rhs.hpp
+++ b/cpp/particles/rhs/tracer_rhs.hpp
@@ -52,7 +52,8 @@ class tracer_rhs: public abstract_particle_rhs
int operator()(
double t,
abstract_particle_set &pset,
- particle_rnumber *result) const
+ particle_rnumber *result,
+ std::vector<std::unique_ptr<abstract_particle_set::particle_rnumber[]>> &additional_data) const
{
// interpolation adds on top of existing values, so result must be cleared.
std::fill_n(result, pset.getLocalNumberOfParticles()*3, 0);
diff --git a/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp b/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
index a21a8113b551f2f08927765977ed2c3d8337b341..200fcdeb80c5bbe39568264b83b7c2591d4076f5 100644
--- a/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
+++ b/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
@@ -47,29 +47,31 @@ class tracer_with_collision_counter_rhs: public tracer_rhs<rnumber, be, tt>
int operator()(
double t,
abstract_particle_set &pset,
- abstract_particle_set::particle_rnumber *result) const
+ abstract_particle_set::particle_rnumber *result,
+ std::vector<std::unique_ptr<abstract_particle_set::particle_rnumber[]>> &additional_data) const
{
// interpolation adds on top of existing values, so result must be cleared.
std::fill_n(result, pset.getLocalNumberOfParticles()*3, 0);
int interpolation_result = (*(this->velocity))(t, pset, result);
- std::vector<std::unique_ptr<abstract_particle_set::particle_rnumber[]>> bla;
- bla.resize(1);
- bla[0].reset(new abstract_particle_set::particle_rnumber[pset.getLocalNumberOfParticles()*pset.getStateSize()]);
+ 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)
{
- bla[0][idx] = result[idx];
+ additional_data[additional_data.size()-1][idx] = result[idx];
});
pset.template applyP2PKernel<3, p2p_ghost_collisions<abstract_particle_set::particle_rnumber, abstract_particle_set::partsize_t>>(
this->p2p_gc,
- bla);
+ additional_data);
// copy shuffled rhs values
pset.LOOP(
[&](const abstract_particle_set::partsize_t idx)
{
- result[idx] = bla[0][idx];
+ result[idx] = additional_data[additional_data.size()-1][idx];
});
+ // clear temporary array
+ additional_data.pop_back();
return interpolation_result;
}