118 std::vector<Utils::Vector3d> &
res) {
119 for (
int i : {-1, 0, 1}) {
120 for (
int j : {-1, 0, 1}) {
121 for (
int k : {-1, 0, 1}) {
128 auto normal_shift = (pos_shifted - pos_folded)[le.shear_plane_normal];
129 if (normal_shift > std::numeric_limits<double>::epsilon())
130 pos_shifted[le.shear_direction] += le.pos_offset;
131 if (normal_shift < -std::numeric_limits<double>::epsilon())
132 pos_shifted[le.shear_direction] -= le.pos_offset;
135 if (
in_box(pos_shifted, halo_lower_corner, halo_upper_corner)) {
136 res.emplace_back(pos_shifted);
182 if (particles.empty()) {
185 enum coupling_modes { none, particle_force, swimmer_force_on_fluid };
186 auto const halo = 0.5 * m_lb.
get_agrid();
188 auto const fully_inside_lower = m_local_box.
my_left() + 2. * halo_vec;
189 auto const fully_inside_upper = m_local_box.
my_right() - 2. * halo_vec;
190 auto const halo_lower_corner = m_local_box.
my_left() - halo_vec;
191 auto const halo_upper_corner = m_local_box.
my_right() + halo_vec;
192 std::vector<Utils::Vector3d> positions_velocity_coupling;
193 std::vector<Utils::Vector3d> positions_force_coupling;
194 std::vector<Utils::Vector3d> force_coupling_forces;
195 std::vector<uint8_t> positions_force_coupling_counter;
196 std::vector<Particle *> coupled_particles;
197 for (
auto ptr : particles) {
201 if (
in_box(folded_pos, fully_inside_lower, fully_inside_upper)) {
202 positions_force_coupling.emplace_back(folded_pos);
204 auto const old_size = positions_force_coupling.size();
206 m_box_geo, positions_force_coupling);
207 auto const new_size = positions_force_coupling.size();
208 span_size =
static_cast<uint8_t
>(new_size - old_size);
210 auto coupling_mode = none;
212 if (p.swimming().is_engine_force_on_fluid) {
213 coupling_mode = swimmer_force_on_fluid;
216 if (coupling_mode == none) {
217 for (
auto end = positions_force_coupling.end(), it = end - span_size;
219 auto const &
pos = *it;
220 if (
pos >= halo_lower_corner and
pos < halo_upper_corner) {
221 positions_velocity_coupling.emplace_back(
pos);
222 coupling_mode = particle_force;
227 if (coupling_mode == none) {
228 positions_force_coupling.erase(positions_force_coupling.end() - span_size,
229 positions_force_coupling.end());
231 coupled_particles.emplace_back(ptr);
232 positions_force_coupling_counter.emplace_back(span_size);
236 if (coupled_particles.empty()) {
239 auto interpolated_velocities =
242 auto const &domain_lower_corner = m_local_box.
my_left();
243 auto const &domain_upper_corner = m_local_box.
my_right();
244 auto it_interpolated_velocities = interpolated_velocities.begin();
245 auto it_positions_force_coupling = positions_force_coupling.begin();
246 auto it_positions_force_coupling_counter =
247 positions_force_coupling_counter.begin();
248 for (
auto ptr : coupled_particles) {
250 auto coupling_mode = particle_force;
252 if (p.swimming().is_engine_force_on_fluid) {
253 coupling_mode = swimmer_force_on_fluid;
257 if (coupling_mode == particle_force) {
258 auto const &v_fluid = *it_interpolated_velocities;
261 force_on_particle = drag_force + random_force;
262 ++it_interpolated_velocities;
265 auto force_on_fluid = -force_on_particle;
267 if (coupling_mode == swimmer_force_on_fluid) {
268 force_on_fluid = p.calc_director() * p.swimming().f_swim;
272 auto const span_size = *it_positions_force_coupling_counter;
273 ++it_positions_force_coupling_counter;
274 for (uint8_t i{0
u}; i < span_size; ++i) {
275 auto &
pos = *it_positions_force_coupling;
276 if (
pos >= domain_lower_corner and
pos < domain_upper_corner) {
279 p.force() += force_on_particle;
281 force_coupling_forces.emplace_back(force_on_fluid);
282 ++it_positions_force_coupling;
306 CALI_CXX_MARK_FUNCTION;
308 assert(thermostat->lb !=
nullptr);
309 if (thermostat->lb->couple_to_md) {
310 if (not lb.is_solver_set()) {
314 auto const real_particles = cell_structure->local_particles();
315 auto const ghost_particles = cell_structure->ghost_particles();
318 std::vector<Particle *> particles{};
319 for (
auto const *particle_range : {&real_particles, &ghost_particles}) {
320 for (
auto &p : *particle_range) {
321 if (not
LB::is_tracer(p) and bookkeeping.should_be_coupled(p)) {
322#if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY)
325 particles.emplace_back(&p);
329 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 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