diff --git a/bfps/cpp/full_code/NSVEparticlesP2P.cpp b/bfps/cpp/full_code/NSVEparticlesP2P.cpp
index 12d77149a1c36b0393b89b749c8fe515e15d697b..3a8312f189c4d4e84b38c89e694803eef13efbd5 100644
--- a/bfps/cpp/full_code/NSVEparticlesP2P.cpp
+++ b/bfps/cpp/full_code/NSVEparticlesP2P.cpp
@@ -12,7 +12,10 @@ int NSVEparticlesP2P<rnumber>::initialize(void)
this->NSVE<rnumber>::initialize();
p2p_computer<double, long long int> current_p2p_computer;
+ // TODO: particle interactions are switched off manually for testing purposes.
+ // this needs to be fixed once particle interactions can be properly resolved.
current_p2p_computer.setEnable(enable_p2p);
+ //current_p2p_computer.setEnable(false);
particles_inner_computer<double, long long int> current_particles_inner_computer(inner_v0);
current_particles_inner_computer.setEnable(enable_inner);
@@ -59,7 +62,12 @@ int NSVEparticlesP2P<rnumber>::step(void)
this->fs->compute_velocity(this->fs->cvorticity);
this->fs->cvelocity->ift();
if(enable_vorticity_omega){
- this->ps->completeLoopWithVorticity(this->dt, *this->fs->cvorticity);
+ *this->tmp_vec_field = this->fs->cvorticity->get_cdata();
+ this->tmp_vec_field->ift();
+ std::unique_ptr<double[]> pdata;
+ pdata.reset(new double[ps->getLocalNbParticles()*3]);
+ this->ps->sample_compute_field(*this->tmp_vec_field, pdata.get());
+ this->ps->completeLoopWithVorticity(this->dt, pdata.get());
}
else{
this->ps->completeLoop(this->dt);
diff --git a/bfps/cpp/particles/abstract_particles_system.hpp b/bfps/cpp/particles/abstract_particles_system.hpp
index 456a756381f341a1c59dc24b468f5afbc47562ab..859432c02a090390a2a304984fd64d39842fe532 100644
--- a/bfps/cpp/particles/abstract_particles_system.hpp
+++ b/bfps/cpp/particles/abstract_particles_system.hpp
@@ -18,6 +18,8 @@ public:
virtual void compute_particles_inner() = 0;
+ virtual void enforce_unit_orientation() = 0;
+
virtual void compute_particles_inner(const real_number particle_extra_rhs[]) = 0;
virtual void move(const real_number dt) = 0;
diff --git a/bfps/cpp/particles/particles_inner_computer.hpp b/bfps/cpp/particles/particles_inner_computer.hpp
index b20abfdab437e008454de1102a1444cc8cbacc71..d58981e6204b64e5cff8b7f5e626ee08b590049d 100644
--- a/bfps/cpp/particles/particles_inner_computer.hpp
+++ b/bfps/cpp/particles/particles_inner_computer.hpp
@@ -14,7 +14,7 @@ public:
}
template <int size_particle_positions, int size_particle_rhs>
- void compute_interaction(const partsize_t nb_particles, real_number pos_part[], real_number rhs_part[]) const{
+ void compute_interaction(const partsize_t nb_particles, const real_number pos_part[], real_number rhs_part[]) const{
static_assert(size_particle_positions == 6, "This kernel works only with 6 values for one particle's position");
static_assert(size_particle_rhs == 6, "This kernel works only with 6 values per particle's rhs");
@@ -24,34 +24,85 @@ public:
rhs_part[idx_part*size_particle_rhs + IDX_X] += pos_part[idx_part*size_particle_positions + 3+IDX_X]*v0;
rhs_part[idx_part*size_particle_rhs + IDX_Y] += pos_part[idx_part*size_particle_positions + 3+IDX_Y]*v0;
rhs_part[idx_part*size_particle_rhs + IDX_Z] += pos_part[idx_part*size_particle_positions + 3+IDX_Z]*v0;
+ }
+ }
+
+ template <int size_particle_positions, int size_particle_rhs>
+ void enforce_unit_orientation(const partsize_t nb_particles, const real_number pos_part[], real_number rhs_part[]) const{
+ static_assert(size_particle_positions == 6, "This kernel works only with 6 values for one particle's position");
+ static_assert(size_particle_rhs == 6, "This kernel works only with 6 values per particle's rhs");
- real_number alpha[3] = {0, 0, 0};
+ #pragma omp parallel for
+ for(partsize_t idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
+ const partsize_t idx0 = idx_part*size_particle_positions + 3;
+ const partsize_t idx1 = idx_part*size_particle_rhs + 3;
+ //real_number alpha[3] = {0, 0, 0};
+ // check that orientation is unit vector:
+ real_number orientation_size = sqrt(
+ pos_part[idx0+IDX_X]*pos_part[idx1+IDX_X] +
+ pos_part[idx0+IDX_Y]*pos_part[idx1+IDX_Y] +
+ pos_part[idx0+IDX_Z]*pos_part[idx1+IDX_Z]);
+ assert(orientation_size > 0.99);
+ assert(orientation_size < 1.01);
// I call "rotation" to be the right hand side of the orientation part of the ODE
// project rotation on orientation:
real_number projection = (
- pos_part[idx_part*size_particle_positions + 3+IDX_X]*rhs_part[idx_part*size_particle_rhs + 3+IDX_X] +
- pos_part[idx_part*size_particle_positions + 3+IDX_Y]*rhs_part[idx_part*size_particle_rhs + 3+IDX_Y] +
- pos_part[idx_part*size_particle_positions + 3+IDX_Z]*rhs_part[idx_part*size_particle_rhs + 3+IDX_Z]);
- // alpha is the vector that makes rotation perpendicular to orientation.
- // note that the following three lines assume the current orientation is a unit vector.
- alpha[IDX_X] = -pos_part[idx_part*size_particle_positions + 3+IDX_X]*projection;
- alpha[IDX_Y] = -pos_part[idx_part*size_particle_positions + 3+IDX_Z]*projection;
- alpha[IDX_Z] = -pos_part[idx_part*size_particle_positions + 3+IDX_Y]*projection;
-
- // now add alpha term to orientation ODE right-hand side.
- rhs_part[idx_part*size_particle_rhs + 3+IDX_X] += alpha[IDX_X];
- rhs_part[idx_part*size_particle_rhs + 3+IDX_Y] += alpha[IDX_Y];
- rhs_part[idx_part*size_particle_rhs + 3+IDX_Z] += alpha[IDX_Z];
+ pos_part[idx0+IDX_X]*rhs_part[idx1+IDX_X] +
+ pos_part[idx0+IDX_Y]*rhs_part[idx1+IDX_Y] +
+ pos_part[idx0+IDX_Z]*rhs_part[idx1+IDX_Z]);
+ //// alpha is the vector that makes rotation perpendicular to orientation.
+ //// note that the following three lines assume the current orientation is a unit vector.
+ //alpha[IDX_X] = -;
+ //alpha[IDX_Y] = -;
+ //alpha[IDX_Z] = -;
+ // DEBUG_MSG("projection = %g\n"
+ // "pos_part[%d] = %g, pos_part[%d] = %g, pos_part[%d] = %g\n"
+ // "rhs_part[%d] = %g, rhs_part[%d] = %g, rhs_part[%d] = %g\n",
+ // projection,
+ // IDX_X, pos_part[idx0 + IDX_X],
+ // IDX_Y, pos_part[idx0 + IDX_Y],
+ // IDX_Z, pos_part[idx0 + IDX_Z],
+ // IDX_X, rhs_part[idx1 + IDX_X],
+ // IDX_Y, rhs_part[idx1 + IDX_Y],
+ // IDX_Z, rhs_part[idx1 + IDX_Z]);
+
+ // now remove parallel bit.
+ rhs_part[idx1+IDX_X] -= pos_part[idx0+IDX_X]*projection;
+ rhs_part[idx1+IDX_Y] -= pos_part[idx0+IDX_Y]*projection;
+ rhs_part[idx1+IDX_Z] -= pos_part[idx0+IDX_Z]*projection;
+
+ // compute dot product between orientation and orientation change
+ real_number dotproduct = (
+ rhs_part[idx1 + IDX_X]*pos_part[idx0 + IDX_X] +
+ rhs_part[idx1 + IDX_Y]*pos_part[idx0 + IDX_Y] +
+ rhs_part[idx1 + IDX_Z]*pos_part[idx0 + IDX_Z]);
+ if (dotproduct > 0.1)
+ {
+ DEBUG_MSG("dotproduct = %g, projection = %g\n"
+ "pos_part[%d] = %g, pos_part[%d] = %g, pos_part[%d] = %g\n"
+ "rhs_part[%d] = %g, rhs_part[%d] = %g, rhs_part[%d] = %g\n",
+ dotproduct,
+ projection,
+ IDX_X, pos_part[idx0 + IDX_X],
+ IDX_Y, pos_part[idx0 + IDX_Y],
+ IDX_Z, pos_part[idx0 + IDX_Z],
+ IDX_X, rhs_part[idx1 + IDX_X],
+ IDX_Y, rhs_part[idx1 + IDX_Y],
+ IDX_Z, rhs_part[idx1 + IDX_Z]);
+ }
+ assert(dotproduct <= 0.1);
}
}
template <int size_particle_positions, int size_particle_rhs, int size_particle_rhs_extra>
- void compute_interaction_with_extra(const partsize_t nb_particles, real_number pos_part[], real_number rhs_part[],
+ void compute_interaction_with_extra(const partsize_t nb_particles, const real_number pos_part[], real_number rhs_part[],
const real_number rhs_part_extra[]) const{
static_assert(size_particle_rhs_extra == 3, "This kernel works only with 3 values for one particle's rhs extra");
+ // call plain compute_interaction first
compute_interaction<size_particle_positions, size_particle_rhs>(nb_particles, pos_part, rhs_part);
+ // now add vorticity term
#pragma omp parallel for
for(partsize_t idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
// Cross product vorticity/orientation
diff --git a/bfps/cpp/particles/particles_inner_computer_empty.hpp b/bfps/cpp/particles/particles_inner_computer_empty.hpp
index 1bd3b1ecc1edbff73fa003de5bfd0932143c5b82..263d8b178ef64d034945442021700cfa702671d5 100644
--- a/bfps/cpp/particles/particles_inner_computer_empty.hpp
+++ b/bfps/cpp/particles/particles_inner_computer_empty.hpp
@@ -11,6 +11,10 @@ public:
void compute_interaction(const partsize_t /*nb_particles*/, real_number /*pos_part*/[], real_number /*rhs_part*/[]) const{
}
+ template <int size_particle_positions, int size_particle_rhs>
+ void enforce_unit_orientation(const partsize_t /*nb_particles*/, real_number /*pos_part*/[], real_number /*rhs_part*/[]) const{
+ }
+
template <int size_particle_positions, int size_particle_rhs, int size_particle_rhs_extra>
void compute_interaction_with_extra(const partsize_t /*nb_particles*/, real_number /*pos_part*/[], real_number /*rhs_part*/[],
const real_number /*rhs_part_extra*/[]) const{
diff --git a/bfps/cpp/particles/particles_system.hpp b/bfps/cpp/particles/particles_system.hpp
index 24bd228dc2a4b49fad492ec62335267b60cfbb02..a5d7878f49447f176d3766aed8c3d96f86dc38ec 100644
--- a/bfps/cpp/particles/particles_system.hpp
+++ b/bfps/cpp/particles/particles_system.hpp
@@ -162,6 +162,14 @@ public:
}
}
+ void enforce_unit_orientation() final {
+ if(computer_particules_inner.isEnable() == true){
+ TIMEZONE("particles_system::enforce_unit_orientation");
+ computer_particules_inner.template enforce_unit_orientation<size_particle_positions, size_particle_rhs>(
+ my_nb_particles, my_particles_positions.get(), my_particles_rhs.front().get());
+ }
+ }
+
void compute_particles_inner(const real_number particle_extra_rhs[]) final {
if(computer_particules_inner.isEnable() == true){
TIMEZONE("particles_system::compute_particles_inner");
@@ -252,6 +260,7 @@ public:
compute();
compute_p2p();
compute_particles_inner();
+ enforce_unit_orientation();
move(dt);
redistribute();
inc_step_idx();
@@ -264,6 +273,7 @@ public:
compute();
compute_p2p();
compute_particles_inner(particle_extra_rhs);
+ enforce_unit_orientation();
move(dt);
redistribute();
inc_step_idx();
diff --git a/bfps/test/test_particles.py b/bfps/test/test_particles.py
index 79595a23a89dd1ad22e0621f28ca69f0f87a242c..c3d4c41507da778005c52caca768f871581d7493 100644
--- a/bfps/test/test_particles.py
+++ b/bfps/test/test_particles.py
@@ -88,30 +88,23 @@ def main():
np.maximum(np.abs(x),
np.abs(s[..., 3:])))))
assert(distance < 1e-14)
- # print movement
- x0 = pf['tracers0/position/0'].value
- x1 = pf['tracers0/position/32'].value
- # compute distance travelled by first particle
- deltax = x1[0] - x0[0]
- print('distance travelled by first particle is ', np.sum(deltax**2)**.5)
- ## code relevant when velocity field is 0 everywhere.
- ## we check to see what happens to the orientation of the particles
- ## show a histogram of the orientations
- #f = plt.figure()
- #a = f.add_subplot(111)
- #for iteration in range(0, niterations*njobs+1, niterations//4):
- # x = pf['tracers0/orientation/{0}'.format(iteration)].value
- # print(x)
- # hist, bins = np.histogram(
- # x.flatten(),
- # bins = 100)
- # bb = (bins[:-1] + bins[1:])/2
- # pp = hist.astype(np.float) / (np.sum(hist) * (bb[1] - bb[0]))
- # a.plot(bb, pp, label = '{0}'.format(iteration))
- #a.legend(loc = 'best')
- #f.tight_layout()
- #f.savefig('full_orientation_histogram.pdf')
- #plt.close(f)
+ # code relevant when velocity field is 0 everywhere.
+ # we check to see what happens to the orientation of the particles
+ # show a histogram of the orientations
+ f = plt.figure()
+ a = f.add_subplot(111)
+ for iteration in range(0, niterations*njobs+1, niterations//4):
+ x = pf['tracers0/orientation/{0}'.format(iteration)].value
+ hist, bins = np.histogram(
+ x.flatten(),
+ bins = 100)
+ bb = (bins[:-1] + bins[1:])/2
+ pp = hist.astype(np.float) / (np.sum(hist) * (bb[1] - bb[0]))
+ a.plot(bb, pp, label = '{0}'.format(iteration))
+ a.legend(loc = 'best')
+ f.tight_layout()
+ f.savefig('full_orientation_histogram.pdf')
+ plt.close(f)
return None
if __name__ == '__main__':
diff --git a/setup.py b/setup.py
index 27a2993f65a36bfc3b9213fe57895ac27edb7864..34b9e16a1c7f5945513480d2b7d884c2d488ac91 100644
--- a/setup.py
+++ b/setup.py
@@ -88,7 +88,8 @@ print('This is bfps version ' + VERSION)
### lists of files and MANIFEST.in
-src_file_list = ['full_code/joint_acc_vel_stats',
+src_file_list = ['full_code/NSVEparticlesP2P',
+ 'full_code/joint_acc_vel_stats',
'full_code/test',
'full_code/filter_test',
'full_code/field_test',
@@ -128,8 +129,7 @@ src_file_list = ['full_code/joint_acc_vel_stats',
'spline_n10',
'Lagrange_polys',
'scope_timer',
- 'full_code/NSVEparticles',
- 'full_code/NSVEparticlesP2P']
+ 'full_code/NSVEparticles']
particle_headers = [
'cpp/particles/abstract_particles_input.hpp',