Commit e7d42ead authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

fix particle vorticity data access

particle system needs vorticity field, then does its own interpolation.
parent 913d5d14
Pipeline #22905 passed with stage
in 11 minutes and 19 seconds
......@@ -14,8 +14,8 @@ int NSVEparticlesP2P<rnumber>::initialize(void)
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);
//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);
......@@ -64,10 +64,7 @@ int NSVEparticlesP2P<rnumber>::step(void)
if(enable_vorticity_omega){
*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());
this->ps->completeLoopWithVorticity(this->dt, *this->tmp_vec_field);
}
else{
this->ps->completeLoop(this->dt);
......
......@@ -20,6 +20,8 @@ public:
virtual void enforce_unit_orientation() = 0;
virtual void add_Lagrange_multipliers() = 0;
virtual void compute_particles_inner(const real_number particle_extra_rhs[]) = 0;
virtual void move(const real_number dt) = 0;
......
......@@ -3,6 +3,7 @@
#include <cstring>
#include <cassert>
#include <iostream>
template <class real_number, class partsize_t>
class particles_inner_computer{
......@@ -27,8 +28,12 @@ public:
}
}
// for given orientation and right-hand-side, recompute right-hand-side such
// that it is perpendicular to the current orientation.
// this is the job of the Lagrange multiplier terms, hence the
// "add_Lagrange_multipliers" name of the method.
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{
void add_Lagrange_multipliers(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");
......@@ -42,6 +47,12 @@ public:
pos_part[idx0+IDX_X]*pos_part[idx0+IDX_X] +
pos_part[idx0+IDX_Y]*pos_part[idx0+IDX_Y] +
pos_part[idx0+IDX_Z]*pos_part[idx0+IDX_Z]);
//DEBUG_MSG("particle ID %d\n"
// "pos_part[%d] = %g, pos_part[%d] = %g, pos_part[%d] = %g\n",
// idx_part,
// IDX_X, pos_part[idx0 + IDX_X],
// IDX_Y, pos_part[idx0 + IDX_Y],
// IDX_Z, pos_part[idx0 + 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
......@@ -89,8 +100,9 @@ public:
IDX_X, rhs_part[idx1 + IDX_X],
IDX_Y, rhs_part[idx1 + IDX_Y],
IDX_Z, rhs_part[idx1 + IDX_Z]);
assert(false);
}
assert(dotproduct <= 0.1);
//assert(dotproduct <= 0.1);
}
}
......@@ -115,6 +127,27 @@ public:
}
}
// meant to be called AFTER executing the time-stepping operation.
// once the particles have been moved, ensure that the orientation is a unit vector.
template <int size_particle_positions>
void enforce_unit_orientation(const partsize_t nb_particles, real_number pos_part[]) const{
static_assert(size_particle_positions == 6, "This kernel works only with 6 values for one particle's position");
#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;
// compute orientation size:
real_number orientation_size = sqrt(
pos_part[idx0+IDX_X]*pos_part[idx0+IDX_X] +
pos_part[idx0+IDX_Y]*pos_part[idx0+IDX_Y] +
pos_part[idx0+IDX_Z]*pos_part[idx0+IDX_Z]);
// now renormalize
pos_part[idx0 + IDX_X] /= orientation_size;
pos_part[idx0 + IDX_Y] /= orientation_size;
pos_part[idx0 + IDX_Z] /= orientation_size;
}
}
bool isEnable() const {
return isActive;
......
......@@ -11,8 +11,12 @@ public:
void compute_interaction(const partsize_t /*nb_particles*/, real_number /*pos_part*/[], real_number /*rhs_part*/[]) const{
}
template <int size_particle_positions>
void enforce_unit_orientation(const partsize_t /*nb_particles*/, real_number /*pos_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{
void add_Lagrange_multipliers(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>
......
......@@ -162,11 +162,19 @@ public:
}
}
void add_Lagrange_multipliers() final {
if(computer_particules_inner.isEnable() == true){
TIMEZONE("particles_system::add_Lagrange_multipliers");
computer_particules_inner.template add_Lagrange_multipliers<size_particle_positions, size_particle_rhs>(
my_nb_particles, my_particles_positions.get(), my_particles_rhs.front().get());
}
}
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());
computer_particules_inner.template enforce_unit_orientation<size_particle_positions>(
my_nb_particles, my_particles_positions.get());
}
}
......@@ -260,8 +268,9 @@ public:
compute();
compute_p2p();
compute_particles_inner();
enforce_unit_orientation();
add_Lagrange_multipliers();
move(dt);
enforce_unit_orientation();
redistribute();
inc_step_idx();
shift_rhs_vectors();
......@@ -273,8 +282,9 @@ public:
compute();
compute_p2p();
compute_particles_inner(particle_extra_rhs);
enforce_unit_orientation();
add_Lagrange_multipliers();
move(dt);
enforce_unit_orientation();
redistribute();
inc_step_idx();
shift_rhs_vectors();
......
......@@ -43,6 +43,22 @@ def main():
pf = h5py.File(
'test_particles.h5',
'r')
# initial condition:
# show a histogram of the orientations
f = plt.figure()
a = f.add_subplot(111)
for iteration in range(1):
x = cf['tracers0/state/{0}'.format(iteration)][:, 3:]
hist, bins = np.histogram(
np.sum(x**2, axis = -1).flatten()**.5,
bins = np.linspace(0, 2, 40))
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('orientation_histogram.pdf')
plt.close(f)
# show a histogram of the positions
f = plt.figure()
a = f.add_subplot(111)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment