Commit 75c2c631 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

Merge branch 'feature/particles-with-own-drag' into develop

parents 815f1cc8 9a4431f7
Pipeline #70581 passed with stage
in 6 minutes and 34 seconds
......@@ -29,6 +29,16 @@
#include <cstring>
#include <cassert>
/** \brief Computes action of Stokes drag on particles.
*
* Either 6 state variables (positions + momenta) or 7 state variables (positions + momenta + value of drag coefficient).
* The case of particles with individual values of the drag coefficient is relevant for the case of plastic collisions,
* when we do not have a priori knowledge of what the mass/size/drag coefficient of the particle is.
*
* `compute_interaction_with_extra` is currently only defined for 6 or 7 state variables.
*
*/
template <class real_number, class partsize_t>
class particles_inner_computer_2nd_order_Stokes{
double drag_coefficient;
......@@ -51,17 +61,31 @@ public:
real_number particle_state[],
real_number particle_rhs[],
const real_number sampled_velocity[]) const{
assert(size_particle_state == 6);
assert(size_particle_rhs == 6);
assert(size_particle_rhs_extra == 3);
#pragma omp parallel for
for(partsize_t idx_part = 0 ; idx_part < number_of_particles ; ++idx_part){
particle_rhs[idx_part*size_particle_rhs + IDXC_X] = particle_state[idx_part*size_particle_state + 3 + IDXC_X];
particle_rhs[idx_part*size_particle_rhs + IDXC_Y] = particle_state[idx_part*size_particle_state + 3 + IDXC_Y];
particle_rhs[idx_part*size_particle_rhs + IDXC_Z] = particle_state[idx_part*size_particle_state + 3 + IDXC_Z];
particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_X] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_X] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_X]);
particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_Y] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_Y] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_Y]);
particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_Z] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_Z] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_Z]);
assert(size_particle_state == size_particle_rhs);
assert((size_particle_state == 6) || (size_particle_state == 7));
if (size_particle_state == 6){
#pragma omp parallel for
for(partsize_t idx_part = 0 ; idx_part < number_of_particles ; ++idx_part){
particle_rhs[idx_part*size_particle_rhs + IDXC_X] = particle_state[idx_part*size_particle_state + 3 + IDXC_X];
particle_rhs[idx_part*size_particle_rhs + IDXC_Y] = particle_state[idx_part*size_particle_state + 3 + IDXC_Y];
particle_rhs[idx_part*size_particle_rhs + IDXC_Z] = particle_state[idx_part*size_particle_state + 3 + IDXC_Z];
particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_X] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_X] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_X]);
particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_Y] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_Y] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_Y]);
particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_Z] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_Z] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_Z]);
}
}
else if (size_particle_state == 7){
// particle stores its own drag coefficient in the 7th state variable
#pragma omp parallel for
for(partsize_t idx_part = 0 ; idx_part < number_of_particles ; ++idx_part){
particle_rhs[idx_part*size_particle_rhs + IDXC_X] = particle_state[idx_part*size_particle_state + 3 + IDXC_X];
particle_rhs[idx_part*size_particle_rhs + IDXC_Y] = particle_state[idx_part*size_particle_state + 3 + IDXC_Y];
particle_rhs[idx_part*size_particle_rhs + IDXC_Z] = particle_state[idx_part*size_particle_state + 3 + IDXC_Z];
particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_X] = - particle_state[idx_part*size_particle_state + 6] * (particle_state[idx_part*size_particle_state + 3 + IDXC_X] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_X]);
particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_Y] = - particle_state[idx_part*size_particle_state + 6] * (particle_state[idx_part*size_particle_state + 3 + IDXC_Y] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_Y]);
particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_Z] = - particle_state[idx_part*size_particle_state + 6] * (particle_state[idx_part*size_particle_state + 3 + IDXC_Z] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_Z]);
}
}
}
......
......@@ -314,8 +314,8 @@ public:
}
void complete2ndOrderLoop(const real_number dt) final {
assert(size_particle_positions == 6);
assert(size_particle_rhs == 6);
assert((size_particle_positions == 6) || (size_particle_positions == 7));
assert(size_particle_rhs == size_particle_positions);
std::unique_ptr<real_number[]> sampled_velocity(new real_number[getLocalNbParticles()*3]());
std::fill_n(sampled_velocity.get(), 3*getLocalNbParticles(), 0);
sample_compute_field(default_field, sampled_velocity.get());
......
Markdown is supported
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