217 if (particles.empty()) {
220 enum coupling_modes { none, particle_force, swimmer_force_on_fluid };
221 auto const halo = 0.5 * m_lb.
get_agrid();
223 auto const fully_inside_lower = m_local_box.
my_left() + 2. * halo_vec;
224 auto const fully_inside_upper = m_local_box.
my_right() - 2. * halo_vec;
225 auto const halo_lower_corner = m_local_box.
my_left() - halo_vec;
226 auto const halo_upper_corner = m_local_box.
my_right() + halo_vec;
227 std::vector<Utils::Vector3d> positions_velocity_coupling;
228 std::vector<Utils::Vector3d> positions_force_coupling;
229 std::vector<Utils::Vector3d> force_coupling_forces;
230 std::vector<uint8_t> positions_force_coupling_counter;
231 std::vector<Particle *> coupled_particles;
232 for (
auto ptr : particles) {
236 if (
in_box(folded_pos, fully_inside_lower, fully_inside_upper)) {
237 positions_force_coupling.emplace_back(folded_pos);
239 auto const old_size = positions_force_coupling.size();
241 m_box_geo, positions_force_coupling);
242 auto const new_size = positions_force_coupling.size();
243 span_size =
static_cast<uint8_t
>(new_size - old_size);
245 auto coupling_mode = none;
247 if (p.swimming().is_engine_force_on_fluid) {
248 coupling_mode = swimmer_force_on_fluid;
251 if (coupling_mode == none) {
252 for (
auto end = positions_force_coupling.end(), it = end - span_size;
254 auto const &pos = *it;
255 if (pos >= halo_lower_corner and pos < halo_upper_corner) {
256 positions_velocity_coupling.emplace_back(pos);
257 coupling_mode = particle_force;
262 if (coupling_mode == none) {
263 positions_force_coupling.erase(positions_force_coupling.end() - span_size,
264 positions_force_coupling.end());
266 coupled_particles.emplace_back(ptr);
267 positions_force_coupling_counter.emplace_back(span_size);
271 if (coupled_particles.empty()) {
274 auto interpolated_velocities =
277 auto const &domain_lower_corner = m_local_box.
my_left();
278 auto const &domain_upper_corner = m_local_box.
my_right();
279 auto it_interpolated_velocities = interpolated_velocities.begin();
280 auto it_positions_force_coupling = positions_force_coupling.begin();
281 auto it_positions_velocity_coupling = positions_velocity_coupling.begin();
282 auto it_positions_force_coupling_counter =
283 positions_force_coupling_counter.begin();
284 for (
auto ptr : coupled_particles) {
286 auto coupling_mode = particle_force;
288 if (p.swimming().is_engine_force_on_fluid) {
289 coupling_mode = swimmer_force_on_fluid;
293 if (coupling_mode == particle_force) {
294 auto &v_fluid = *it_interpolated_velocities;
300 *it_positions_velocity_coupling, p.pos(), m_box_geo);
301 v_fluid += vel_correction;
305 force_on_particle = drag_force + random_force;
306 ++it_interpolated_velocities;
307 ++it_positions_velocity_coupling;
310 auto force_on_fluid = -force_on_particle;
312 if (coupling_mode == swimmer_force_on_fluid) {
313 force_on_fluid = p.calc_director() * p.swimming().f_swim;
317 auto const span_size = *it_positions_force_coupling_counter;
318 ++it_positions_force_coupling_counter;
319 for (uint8_t i{0u}; i < span_size; ++i) {
320 auto &pos = *it_positions_force_coupling;
321 if (pos >= domain_lower_corner and pos < domain_upper_corner) {
324 p.force() += force_on_particle;
326 force_coupling_forces.emplace_back(force_on_fluid);
327 ++it_positions_force_coupling;
351 CALI_CXX_MARK_FUNCTION;
353 assert(thermostat->lb !=
nullptr);
354 if (thermostat->lb->couple_to_md) {
355 if (not lb.is_solver_set()) {
359 auto const real_particles = cell_structure->local_particles();
360 auto const ghost_particles = cell_structure->ghost_particles();
363 lb.ghost_communication_vel();
364 std::vector<Particle *> particles{};
365 for (
auto const *particle_range : {&real_particles, &ghost_particles}) {
366 for (
auto &p : *particle_range) {
367 if (not
LB::is_tracer(p) and bookkeeping.should_be_coupled(p)) {
368#if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY)
371 particles.emplace_back(&p);
375 coupling.kernel(particles);
void kernel(std::vector< Particle * > const &particles)
Utils::Vector3d get_noise_term(Particle const &p) const
bool in_local_halo(LocalBox const &local_box, Utils::Vector3d const &pos, double agrid)
Check if a position is within the local LB domain plus halo.
static auto lees_edwards_vel_shift(Utils::Vector3d const &pos_shifted_by_box_l, Utils::Vector3d const &orig_pos, BoxGeometry const &box_geo)
static bool in_box(Utils::Vector3d const &pos, Utils::Vector3d const &lower_corner, Utils::Vector3d const &upper_corner)
static bool in_local_domain(LocalBox const &local_box, Utils::Vector3d const &pos, double halo=0.)
Check if a position is within the local box + halo.
static void positions_in_halo_impl(Utils::Vector3d const &pos_folded, Utils::Vector3d const &halo_lower_corner, Utils::Vector3d const &halo_upper_corner, BoxGeometry const &box_geo, std::vector< Utils::Vector3d > &res)
std::vector< Utils::Vector3d > positions_in_halo(Utils::Vector3d const &pos, BoxGeometry const &box_geo, LocalBox const &local_box, double agrid)
Return a vector of positions shifted by +,- box length in each coordinate.
static Utils::Vector3d lb_drag_force(Particle const &p, double lb_gamma, Utils::Vector3d const &v_fluid)
static Thermostat::GammaType lb_handle_particle_anisotropy(Particle const &p, double lb_gamma)
Struct holding all information for one particle.
auto const & gamma() const
auto const & mu_E() const