Skip to content
GitLab
Commit 94f89084 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

add dumb version of ellipsoid inner computer

parent 13164286
Hide whitespace changes
Inline Side-by-side
bfps/cpp/particles/particles_inner_computer.cpp
View file @ 94f89084
......@@ -16,9 +16,9 @@ void particles_inner_computer<real_number, partsize_t>::compute_interaction(
#pragma omp parallel for
for(partsize_t idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
// Add attr × V0 to the field interpolation
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;
rhs_part[idx_part*size_particle_rhs + IDXC_X] += pos_part[idx_part*size_particle_positions + 3+IDXC_X]*v0;
rhs_part[idx_part*size_particle_rhs + IDXC_Y] += pos_part[idx_part*size_particle_positions + 3+IDXC_Y]*v0;
rhs_part[idx_part*size_particle_rhs + IDXC_Z] += pos_part[idx_part*size_particle_positions + 3+IDXC_Z]*v0;
}
}
......@@ -39,31 +39,31 @@ void particles_inner_computer<double, long long>::add_Lagrange_multipliers<6,6>(
const long long idx1 = idx_part*6 + 3;
// check that orientation is unit vector:
double 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]);
pos_part[idx0+IDXC_X]*pos_part[idx0+IDXC_X] +
pos_part[idx0+IDXC_Y]*pos_part[idx0+IDXC_Y] +
pos_part[idx0+IDXC_Z]*pos_part[idx0+IDXC_Z]);
variable_used_only_in_assert(orientation_size);
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:
double projection = (
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]);
pos_part[idx0+IDXC_X]*rhs_part[idx1+IDXC_X] +
pos_part[idx0+IDXC_Y]*rhs_part[idx1+IDXC_Y] +
pos_part[idx0+IDXC_Z]*rhs_part[idx1+IDXC_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;
rhs_part[idx1+IDXC_X] -= pos_part[idx0+IDXC_X]*projection;
rhs_part[idx1+IDXC_Y] -= pos_part[idx0+IDXC_Y]*projection;
rhs_part[idx1+IDXC_Z] -= pos_part[idx0+IDXC_Z]*projection;
// DEBUG
// sanity check, for debugging purposes
// compute dot product between orientation and orientation change
//double 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]);
// rhs_part[idx1 + IDXC_X]*pos_part[idx0 + IDXC_X] +
// rhs_part[idx1 + IDXC_Y]*pos_part[idx0 + IDXC_Y] +
// rhs_part[idx1 + IDXC_Z]*pos_part[idx0 + IDXC_Z]);
//if (dotproduct > 0.1)
//{
// DEBUG_MSG("dotproduct = %g, projection = %g\n"
......@@ -71,12 +71,12 @@ void particles_inner_computer<double, long long>::add_Lagrange_multipliers<6,6>(
// "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]);
// IDXC_X, pos_part[idx0 + IDXC_X],
// IDXC_Y, pos_part[idx0 + IDXC_Y],
// IDXC_Z, pos_part[idx0 + IDXC_Z],
// IDXC_X, rhs_part[idx1 + IDXC_X],
// IDXC_Y, rhs_part[idx1 + IDXC_Y],
// IDXC_Z, rhs_part[idx1 + IDXC_Z]);
// assert(false);
//}
//assert(dotproduct <= 0.1);
......@@ -97,12 +97,35 @@ void particles_inner_computer<double, long long>::compute_interaction_with_extra
#pragma omp parallel for
for(long long idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
// Cross product vorticity/orientation
rhs_part[idx_part*6 + 3+IDX_X] += 0.5*(rhs_part_extra[idx_part*3 + IDX_Y]*pos_part[idx_part*6 + 3+IDX_Z] -
rhs_part_extra[idx_part*3 + IDX_Z]*pos_part[idx_part*6 + 3+IDX_Y]);
rhs_part[idx_part*6 + 3+IDX_Y] += 0.5*(rhs_part_extra[idx_part*3 + IDX_Z]*pos_part[idx_part*6 + 3+IDX_X] -
rhs_part_extra[idx_part*3 + IDX_X]*pos_part[idx_part*6 + 3+IDX_Z]);
rhs_part[idx_part*6 + 3+IDX_Z] += 0.5*(rhs_part_extra[idx_part*3 + IDX_X]*pos_part[idx_part*6 + 3+IDX_Y] -
rhs_part_extra[idx_part*3 + IDX_Y]*pos_part[idx_part*6 + 3+IDX_X]);
rhs_part[idx_part*6 + 3+IDXC_X] += 0.5*(rhs_part_extra[idx_part*3 + IDXC_Y]*pos_part[idx_part*6 + 3+IDXC_Z] -
rhs_part_extra[idx_part*3 + IDXC_Z]*pos_part[idx_part*6 + 3+IDXC_Y]);
rhs_part[idx_part*6 + 3+IDXC_Y] += 0.5*(rhs_part_extra[idx_part*3 + IDXC_Z]*pos_part[idx_part*6 + 3+IDXC_X] -
rhs_part_extra[idx_part*3 + IDXC_X]*pos_part[idx_part*6 + 3+IDXC_Z]);
rhs_part[idx_part*6 + 3+IDXC_Z] += 0.5*(rhs_part_extra[idx_part*3 + IDXC_X]*pos_part[idx_part*6 + 3+IDXC_Y] -
rhs_part_extra[idx_part*3 + IDXC_Y]*pos_part[idx_part*6 + 3+IDXC_X]);
}
}
template <>
template <>
void particles_inner_computer<double, long long>::compute_interaction_with_extra<6,6,9>(
const long long nb_particles,
const double pos_part[],
double rhs_part[],
const double rhs_part_extra[]) const{
// call plain compute_interaction first
compute_interaction<6, 6>(nb_particles, pos_part, rhs_part);
// now add vorticity term
#pragma omp parallel for
for(long long idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
// Cross product vorticity/orientation
rhs_part[idx_part*6 + 3+IDXC_X] += 0.5*(rhs_part_extra[idx_part*9 + IDXC_Y]*pos_part[idx_part*6 + 3+IDXC_Z] -
rhs_part_extra[idx_part*9 + IDXC_Z]*pos_part[idx_part*6 + 3+IDXC_Y]);
rhs_part[idx_part*6 + 3+IDXC_Y] += 0.5*(rhs_part_extra[idx_part*9 + IDXC_Z]*pos_part[idx_part*6 + 3+IDXC_X] -
rhs_part_extra[idx_part*9 + IDXC_X]*pos_part[idx_part*6 + 3+IDXC_Z]);
rhs_part[idx_part*6 + 3+IDXC_Z] += 0.5*(rhs_part_extra[idx_part*9 + IDXC_X]*pos_part[idx_part*6 + 3+IDXC_Y] -
rhs_part_extra[idx_part*9 + IDXC_Y]*pos_part[idx_part*6 + 3+IDXC_X]);
}
}
......@@ -119,13 +142,13 @@ void particles_inner_computer<double, long long>::enforce_unit_orientation<6>(
const long long idx0 = idx_part*6 + 3;
// compute orientation size:
double 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]);
pos_part[idx0+IDXC_X]*pos_part[idx0+IDXC_X] +
pos_part[idx0+IDXC_Y]*pos_part[idx0+IDXC_Y] +
pos_part[idx0+IDXC_Z]*pos_part[idx0+IDXC_Z]);
// now renormalize
pos_part[idx0 + IDX_X] /= orientation_size;
pos_part[idx0 + IDX_Y] /= orientation_size;
pos_part[idx0 + IDX_Z] /= orientation_size;
pos_part[idx0 + IDXC_X] /= orientation_size;
pos_part[idx0 + IDXC_Y] /= orientation_size;
pos_part[idx0 + IDXC_Z] /= orientation_size;
}
}
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