diff --git a/cpp/full_code/test_tracer_set.cpp b/cpp/full_code/test_tracer_set.cpp
index 836a3311917b00ad5c8cc29d9c0dc1e656a91335..74112e870fcfdaaac173f23aafae28b17fe8e8bb 100644
--- a/cpp/full_code/test_tracer_set.cpp
+++ b/cpp/full_code/test_tracer_set.cpp
@@ -80,15 +80,19 @@ int test_tracer_set<rnumber>::do_work(void)
this->dkz);
particles_input_random<long long int, double, 3> bla(
- this->comm, 1000, 1, pset.getSpatialLowLimitZ(), pset.getSpatialUpLimitZ());
+ this->comm,
+ 1000,
+ 1,
+ pset.getSpatialLowLimitZ(),
+ pset.getSpatialUpLimitZ());
fti->set_field(vec_field);
trhs->setVelocity(fti);
pset.init(bla);
- //particle_solver psolver(pset);
+ particle_solver psolver(pset);
- //psolver.Euler(0.1, *trhs);
+ psolver.Euler(0.1, *trhs);
// deallocate
delete trhs;
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index 6f7c82c1ca3a5f66aff90c5f4dd241d4a53b45e0..7330cfe80a8f50490879eb38c8dbb8ac6f62c949 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -135,6 +135,10 @@ class particle_set: public abstract_particle_set
assert(IDXC_Y == 1);
assert(IDXC_Z == 2);
+ DEBUG_MSG("current partition interval %d %d\n",
+ current_partition_interval.first,
+ current_partition_interval.second);
+
this->number_particles_per_partition.reset(new partsize_t[partition_interval_size]);
this->offset_particles_for_partition.reset(new partsize_t[partition_interval_size+1]);
}
@@ -279,22 +283,30 @@ class particle_set: public abstract_particle_set
this->local_index = particles_input.getMyParticlesIndexes();
this->local_number_of_particles = particles_input.getLocalNbParticles();
- //particles_utils::partition_extra_z<partsize_t, state_size>(&this->local_state[0], this->local_number_of_particles, partition_interval_size,
- // this->local_number_of_particles_per_partition.get(), current_offset_particles_for_partition.get(),
- //[&](const particle_rnumber& z_pos){
- // const int partition_level = this->pinterpolator.pbc_field_layer(z_pos, IDXC_Z);
- // assert(current_partition_interval.first <= partition_level && partition_level < current_partition_interval.second);
- // return partition_level - current_partition_interval.first;
- //},
- //[&](const partsize_t idx1, const partsize_t idx2){
- // std::swap(this->local_index[idx1], this->local_index[idx2]);
- // for(int idx_rhs = 0 ; idx_rhs < int(my_particles_rhs.size()) ; ++idx_rhs){
- // for(int idx_val = 0 ; idx_val < size_particle_rhs ; ++idx_val){
- // std::swap(my_particles_rhs[idx_rhs][idx1*size_particle_rhs + idx_val],
- // my_particles_rhs[idx_rhs][idx2*size_particle_rhs + idx_val]);
- // }
- // }
- //});
+ for (partsize_t idx=0; idx < this->local_number_of_particles; idx++)
+ {
+ DEBUG_MSG_WAIT(this->mpi_comm, "idx = %d, z = %g\n",
+ this->local_index[idx],
+ this->local_state[idx*state_size+2]);
+ }
+
+ particles_utils::partition_extra_z<partsize_t, state_size>(
+ &this->local_state[0],
+ this->local_number_of_particles,
+ partition_interval_size,
+ this->number_particles_per_partition.get(),
+ this->offset_particles_for_partition.get(),
+ [&](const particle_rnumber& z_pos){
+ const int partition_level = this->pinterpolator.pbc_field_layer(z_pos, IDXC_Z);
+ DEBUG_MSG("z_pos = %g, partition_level is %d\n",
+ z_pos,
+ partition_level);
+ assert(current_partition_interval.first <= partition_level && partition_level < current_partition_interval.second);
+ return partition_level - current_partition_interval.first;
+ },
+ [&](const partsize_t idx1, const partsize_t idx2){
+ std::swap(this->local_index[idx1], this->local_index[idx2]);
+ });
return EXIT_SUCCESS;
}
diff --git a/cpp/particles/particles_distr_mpi.hpp b/cpp/particles/particles_distr_mpi.hpp
index 2cc73cd5d732ec3fd5877ae5d5a5fbc694f2f913..17dc121646ae54883a709a8a5009b60a64d209b7 100644
--- a/cpp/particles/particles_distr_mpi.hpp
+++ b/cpp/particles/particles_distr_mpi.hpp
@@ -182,6 +182,10 @@ public:
partsize_t myTotalNbParticles = 0;
for(int idxPartition = 0 ; idxPartition < current_partition_size ; ++idxPartition){
myTotalNbParticles += current_my_nb_particles_per_partition[idxPartition];
+ DEBUG_MSG("idxPartition = %d, current_offset_particles_for_partition = %d, current_my_nb_particles_per_partition = %d\n",
+ idxPartition,
+ current_offset_particles_for_partition[idxPartition],
+ current_my_nb_particles_per_partition[idxPartition]);
current_offset_particles_for_partition[idxPartition+1] = current_offset_particles_for_partition[idxPartition] + current_my_nb_particles_per_partition[idxPartition];
}
@@ -211,7 +215,10 @@ public:
const int nbPartitionsToSend = std::min(current_partition_size, interpolation_size-(idxLower-1));
assert(nbPartitionsToSend >= 0);
+ DEBUG_MSG("idxLower = %d, nbPartitionsToSend = %d, current_offset_particles_for_partition = %d, current_offset_particles_for_partition[0] = %d\n",
+ idxLower, nbPartitionsToSend, current_offset_particles_for_partition[nbPartitionsToSend], current_offset_particles_for_partition[0]);
const partsize_t nbParticlesToSend = current_offset_particles_for_partition[nbPartitionsToSend] - current_offset_particles_for_partition[0];
+ assert(nbParticlesToSend >= 0);
const int nbPartitionsToRecv = std::min(partition_interval_size_per_proc[destProc], (interpolation_size+1)-(idxLower-1));
assert(nbPartitionsToRecv > 0);
@@ -558,6 +565,11 @@ public:
for(partsize_t idx = 0 ; idx < current_my_nb_particles_per_partition[idxPartition] ; ++idx){
const int partition_level = in_computer.pbc_field_layer((*inout_positions_particles)[(idx+partOffset)*size_particle_positions+IDXC_Z], IDXC_Z);
variable_used_only_in_assert(partition_level);
+ DEBUG_MSG("partition_level = %d, current_parition_interval.first = %d, idxPartition = %d\n",
+ partition_level,
+ current_partition_interval.first,
+ idxPartition
+ );
assert(partition_level == current_partition_interval.first + idxPartition
|| partition_level == (current_partition_interval.first + idxPartition-1+int(field_grid_dim[IDXC_Z]))%int(field_grid_dim[IDXC_Z])
|| partition_level == (current_partition_interval.first + idxPartition+1)%int(field_grid_dim[IDXC_Z]));
diff --git a/cpp/particles/particles_input_random.hpp b/cpp/particles/particles_input_random.hpp
index 7d5bacefddb3c9a81fe93ecd007c7cdcad657549..d862acc30ada993b1f8cefc1cd68a2bf0f274218 100644
--- a/cpp/particles/particles_input_random.hpp
+++ b/cpp/particles/particles_input_random.hpp
@@ -64,6 +64,7 @@ class particles_input_random: public abstract_particles_input<partsize_t, partic
// initialize a separate random number generator for each MPI process
std::mt19937_64 rgen;
std::uniform_real_distribution<particle_rnumber> udist(0.0, 1.0);
+ std::uniform_real_distribution<particle_rnumber> uzdist(my_spatial_low_limit, my_spatial_up_limit);
std::normal_distribution<particle_rnumber> gdist;
// seed random number generators such that no seed is ever repeated if we change the value of rseed.
// basically use a multi-dimensional array indexing technique to come up with actual seed.
@@ -96,10 +97,11 @@ class particles_input_random: public abstract_particles_input<partsize_t, partic
this->nprocs*(this->total_number_of_particles/this->nprocs) +
std::min(partsize_t(this->myrank), partsize_t(this->total_number_of_particles%this->nprocs)));
for (int cc=0; cc < 2; cc++)
- this->local_particle_state[cc] = twopi*udist(rgen);
- this->local_particle_state[2] = udist(rgen)*(my_spatial_up_limit - my_spatial_low_limit) + my_spatial_low_limit;
+ this->local_particle_state[idx*state_size + cc] = twopi*udist(rgen);
+ this->local_particle_state[idx*state_size + 2] = uzdist(rgen);
+ DEBUG_MSG("zpos = %g\n", this->local_particle_state[idx*state_size + 2]);
for (int cc = 3; cc < state_size; cc++)
- this->local_particle_state[cc] = gdist(rgen);
+ this->local_particle_state[idx*state_size + cc] = gdist(rgen);
}
}