35#include <boost/mpi/collectives/all_reduce.hpp>
36#include <boost/range/algorithm.hpp>
50 auto save_pos = [](
Particle &p) { p.pos_last_time_step() = p.pos(); };
52 boost::for_each(particles, save_pos);
53 boost::for_each(ghost_particles, save_pos);
64 auto reset_force = [](
Particle &p) { p.rattle_params().correction.fill(0); };
66 boost::for_each(particles, reset_force);
67 boost::for_each(ghost_particles, reset_force);
83 auto const r_ij2 = r_ij.norm2();
85 if (std::abs(1.0 - r_ij2 / ia_params.
d2) > ia_params.
p_tol) {
88 auto const r_ij_dot = r_ij_t * r_ij;
90 0.50 * (ia_params.
d2 - r_ij2) / r_ij_dot / (p1.
mass() + p2.
mass());
92 auto const pos_corr = G * r_ij_t;
111template <
typename Kernel>
116 bool correction =
false;
117 cs.
bond_loop([&correction, &kernel, &box_geo, &bonded_ias](
118 Particle &p1,
int bond_id, std::span<Particle *> partners) {
119 auto const &iaparams = *bonded_ias.
at(bond_id);
121 if (
auto const *bond = boost::get<RigidBond>(&iaparams)) {
122 auto const corrected = kernel(*bond, box_geo, p1, *partners[0]);
140 boost::for_each(particles, [](
Particle &p) {
160 boost::mpi::all_reduce(
comm_cart, repeat_, std::logical_or<bool>());
176 auto const resort_level =
196 auto const v_ij = p1.
v() - p2.
v();
199 auto const v_proj = v_ij * r_ij;
200 if (std::abs(v_proj) > ia_params.
v_tol) {
201 auto const K = v_proj / ia_params.
d2 / (p1.
mass() + p2.
mass());
203 auto const vel_corr = K * r_ij;
220 boost::for_each(particles,
237 boost::mpi::all_reduce(
comm_cart, repeat_, std::logical_or<bool>());
Data structures for bonded interactions.
container for bonded interactions.
mapped_type at(key_type const &key) const
Utils::Vector< T, 3 > get_mi_vector(const Utils::Vector< T, 3 > &a, const Utils::Vector< T, 3 > &b) const
Get the minimum-image vector between two coordinates.
boost::mpi::communicator comm_cart
The communicator.
#define SHAKE_MAX_ITERATIONS
Maximal number of iterations in the RATTLE algorithm before it bails out.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
@ DATA_PART_MOMENTUM
Particle::m.
@ DATA_PART_PROPERTIES
Particle::p.
@ DATA_PART_POSITION
Particle::r.
void correct_velocity_shake(CellStructure &cs, BoxGeometry const &box_geo, BondedInteractionsMap const &bonded_ias)
Correction of current velocities using RATTLE algorithm.
void save_old_position(const ParticleRange &particles, const ParticleRange &ghost_particles)
copy current position
void correct_position_shake(CellStructure &cs, BoxGeometry const &box_geo, BondedInteractionsMap const &bonded_ias)
Propagate velocity and position while using SHAKE algorithm for bond constraint.
static void init_correction_vector(const ParticleRange &particles, const ParticleRange &ghost_particles)
reset correction vectors to zero
static void apply_positional_correction(const ParticleRange &particles)
Apply positional corrections.
static bool calculate_velocity_correction(RigidBond const &ia_params, BoxGeometry const &box_geo, Particle &p1, Particle &p2)
Calculate the velocity correction for the particles.
static void apply_velocity_correction(ParticleRange const &particles)
Apply velocity corrections.
static bool compute_correction_vector(CellStructure &cs, BoxGeometry const &box_geo, BondedInteractionsMap const &bonded_ias, Kernel kernel)
Compute the correction vectors using given kernel.
static bool calculate_positional_correction(RigidBond const &ia_params, BoxGeometry const &box_geo, Particle &p1, Particle &p2)
Calculate the positional correction for the particles.
Definition of the rigid bond data type for the Rattle algorithm.
Describes a cell structure / cell system.
ParticleRange ghost_particles() const
void update_ghosts_and_resort_particle(unsigned data_parts)
Update ghost particles, with particle resort if needed.
void ghosts_update(unsigned data_parts)
Update ghost particles.
bool check_resort_required(Utils::Vector3d const &additional_offset={}) const
Check whether a particle has moved further than half the skin since the last Verlet list update,...
void bond_loop(BondKernel const &bond_kernel)
Bonded pair loop.
void set_resort_particles(Cells::Resort level)
Increase the local resort level at least to level.
ParticleRange local_particles() const
void ghosts_reduce_rattle_correction()
Add rattle corrections from ghost particles to real particles.
Struct holding all information for one particle.
auto const & rattle_params() const
auto const & mass() const
auto const & pos_last_time_step() const
Parameters for the rigid_bond/SHAKE/RATTLE ALGORITHM.
double d2
Square of the length of Constrained Bond.
double v_tol
Velocity Tolerance/Accuracy for termination of RATTLE/SHAKE iterations during velocity corrections.
double p_tol
Positional Tolerance/Accuracy value for termination of RATTLE/SHAKE iterations during position correc...