ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
particle_coupling.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
21#include "BoxGeometry.hpp"
22#include "LocalBox.hpp"
23#include "Particle.hpp"
25#include "communication.hpp"
26#include "config/config.hpp"
27#include "errorhandling.hpp"
28#include "random.hpp"
29#include "system/System.hpp"
30#include "thermostat.hpp"
31
32#include <utils/Counter.hpp>
33#include <utils/Vector.hpp>
34
35#include <boost/mpi.hpp>
36
37#ifdef CALIPER
38#include <caliper/cali.h>
39#endif
40
41#include <cmath>
42#include <cstdint>
43#include <initializer_list>
44#include <limits>
45#include <stdexcept>
46#include <vector>
47
49 double lb_gamma) {
50#ifdef THERMOSTAT_PER_PARTICLE
51 auto const &partcl_gamma = p.gamma();
52#ifdef PARTICLE_ANISOTROPY
53 auto const default_gamma = Thermostat::GammaType::broadcast(lb_gamma);
54#else
55 auto const default_gamma = lb_gamma;
56#endif // PARTICLE_ANISOTROPY
57 return Thermostat::handle_particle_gamma(partcl_gamma, default_gamma);
58#else
59 return lb_gamma;
60#endif // THERMOSTAT_PER_PARTICLE
61}
62
63static Utils::Vector3d lb_drag_force(Particle const &p, double lb_gamma,
64 Utils::Vector3d const &v_fluid) {
65#ifdef LB_ELECTROHYDRODYNAMICS
66 auto const v_drift = v_fluid + p.mu_E();
67#else
68 auto const &v_drift = v_fluid;
69#endif
70 auto const gamma = lb_handle_particle_anisotropy(p, lb_gamma);
71
72 /* calculate viscous force (eq. (9) @cite ahlrichs99a) */
73 return Utils::hadamard_product(gamma, v_drift - p.v());
74}
75
76Utils::Vector3d lb_drag_force(LB::Solver const &lb, double lb_gamma,
77 Particle const &p,
78 Utils::Vector3d const &shifted_pos) {
79 auto const v_fluid = lb.get_coupling_interpolated_velocity(shifted_pos);
80 return lb_drag_force(p, lb_gamma, v_fluid);
81}
82
83/**
84 * @brief Check if a position is within the local box + halo.
85 *
86 * @param local_box Local geometry
87 * @param pos Position to check
88 * @param halo Halo
89 *
90 * @return True iff the point is inside of the box up to halo.
91 */
92static bool in_local_domain(LocalBox const &local_box,
93 Utils::Vector3d const &pos, double halo = 0.) {
94 auto const halo_vec = Utils::Vector3d::broadcast(halo);
95 auto const lower_corner = local_box.my_left() - halo_vec;
96 auto const upper_corner = local_box.my_right() + halo_vec;
97
98 return pos >= lower_corner and pos < upper_corner;
99}
100
101static bool in_box(Utils::Vector3d const &pos,
102 Utils::Vector3d const &lower_corner,
103 Utils::Vector3d const &upper_corner) {
104 return pos >= lower_corner and pos < upper_corner;
105}
106
107bool in_local_halo(LocalBox const &local_box, Utils::Vector3d const &pos,
108 double agrid) {
109 auto const halo = 0.5 * agrid;
110 return in_local_domain(local_box, pos, halo);
111}
112
113static void positions_in_halo_impl(Utils::Vector3d const &pos_folded,
114 Utils::Vector3d const &halo_lower_corner,
115 Utils::Vector3d const &halo_upper_corner,
116 BoxGeometry const &box_geo,
117 std::vector<Utils::Vector3d> &res) {
118
119 // Lees-Edwards: pre-calc positional offset folded into the simulation box
120 double folded_le_offset = 0.;
121 if (box_geo.type() == BoxType::LEES_EDWARDS) {
122 auto const &le = box_geo.lees_edwards_bc();
123 folded_le_offset = Algorithm::periodic_fold(
124 le.pos_offset, box_geo.length()[le.shear_direction]);
125 }
126
127 for (int i : {-1, 0, 1}) {
128 for (int j : {-1, 0, 1}) {
129 for (int k : {-1, 0, 1}) {
130 Utils::Vector3d shift{{double(i), double(j), double(k)}};
131
132 // Lees Edwards: folded position incl. LE pos offset
133 // This is needed to ensure that the position from which `pos_shifted`
134 // is calculated below, is always in the primary simulation box.
135 auto with_le_offset = [&](auto pos) {
136 auto const &le = box_geo.lees_edwards_bc();
137 pos[le.shear_direction] = Algorithm::periodic_fold(
138 pos[le.shear_direction] +
139 shift[le.shear_plane_normal] * folded_le_offset,
140 box_geo.length()[le.shear_direction]);
141 return pos;
142 };
143
144 Utils::Vector3d pos_shifted =
145 (box_geo.type() != BoxType::LEES_EDWARDS) ? // no Lees Edwards
146 pos_folded + Utils::hadamard_product(box_geo.length(), shift)
147 : // Lees Edwards
148 with_le_offset(pos_folded) +
149 Utils::hadamard_product(box_geo.length(), shift);
150
151 if (in_box(pos_shifted, halo_lower_corner, halo_upper_corner)) {
152 res.emplace_back(pos_shifted);
153 }
154 }
155 }
156 }
157}
158
159/**
160 * @brief Return a vector of positions shifted by +,- box length in each
161 * coordinate
162 */
163std::vector<Utils::Vector3d> positions_in_halo(Utils::Vector3d const &pos,
164 BoxGeometry const &box_geo,
165 LocalBox const &local_box,
166 double agrid) {
167 auto const halo = 0.5 * agrid;
168 auto const halo_vec = Utils::Vector3d::broadcast(halo);
169 auto const fully_inside_lower = local_box.my_left() + 2. * halo_vec;
170 auto const fully_inside_upper = local_box.my_right() - 2. * halo_vec;
171 auto const pos_folded = box_geo.folded_position(pos);
172 if (in_box(pos_folded, fully_inside_lower, fully_inside_upper)) {
173 return {pos_folded};
174 }
175 auto const halo_lower_corner = local_box.my_left() - halo_vec;
176 auto const halo_upper_corner = local_box.my_right() + halo_vec;
177 std::vector<Utils::Vector3d> res;
178 positions_in_halo_impl(pos_folded, halo_lower_corner, halo_upper_corner,
179 box_geo, res);
180 return res;
181}
182
183static auto lees_edwards_vel_shift(Utils::Vector3d const &pos_shifted_by_box_l,
184 Utils::Vector3d const &orig_pos,
185 BoxGeometry const &box_geo) {
186 Utils::Vector3d vel_shift{{0., 0., 0.}};
187 if (box_geo.type() == BoxType::LEES_EDWARDS) {
188 auto le = box_geo.lees_edwards_bc();
189 auto normal_shift =
190 (pos_shifted_by_box_l - orig_pos)[le.shear_plane_normal];
191 // normal_shift is +,- box_l or 0 up to floating point errors
192 if (normal_shift > std::numeric_limits<double>::epsilon()) {
193 vel_shift[le.shear_direction] -= le.shear_velocity;
194 }
195 if (normal_shift < -std::numeric_limits<double>::epsilon()) {
196 vel_shift[le.shear_direction] += le.shear_velocity;
197 }
198 }
199 return vel_shift;
200}
201
202namespace LB {
203
205 if (not m_thermalized) {
206 return Utils::Vector3d{};
207 }
208 using std::sqrt;
209 using Utils::sqrt;
210 auto const gamma = lb_handle_particle_anisotropy(p, m_thermostat.gamma);
211 auto const noise = Random::noise_uniform<RNGSalt::PARTICLES>(
212 m_thermostat.rng_counter(), m_thermostat.rng_seed(), p.id());
213 return m_noise_pref_wo_gamma * Utils::hadamard_product(sqrt(gamma), noise);
214}
215
216void ParticleCoupling::kernel(std::vector<Particle *> const &particles) {
217 if (particles.empty()) {
218 return;
219 }
220 enum coupling_modes { none, particle_force, swimmer_force_on_fluid };
221 auto const halo = 0.5 * m_lb.get_agrid();
222 auto const halo_vec = Utils::Vector3d::broadcast(halo);
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) {
233 auto &p = *ptr;
234 auto span_size = 1u;
235 auto const folded_pos = m_box_geo.folded_position(p.pos());
236 if (in_box(folded_pos, fully_inside_lower, fully_inside_upper)) {
237 positions_force_coupling.emplace_back(folded_pos);
238 } else {
239 auto const old_size = positions_force_coupling.size();
240 positions_in_halo_impl(folded_pos, halo_lower_corner, halo_upper_corner,
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);
244 }
245 auto coupling_mode = none;
246#ifdef ENGINE
247 if (p.swimming().is_engine_force_on_fluid) {
248 coupling_mode = swimmer_force_on_fluid;
249 }
250#endif
251 if (coupling_mode == none) {
252 for (auto end = positions_force_coupling.end(), it = end - span_size;
253 it != end; ++it) {
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;
258 break;
259 }
260 }
261 }
262 if (coupling_mode == none) {
263 positions_force_coupling.erase(positions_force_coupling.end() - span_size,
264 positions_force_coupling.end());
265 } else {
266 coupled_particles.emplace_back(ptr);
267 positions_force_coupling_counter.emplace_back(span_size);
268 }
269 }
270
271 if (coupled_particles.empty()) {
272 return;
273 }
274 auto interpolated_velocities =
275 m_lb.get_coupling_interpolated_velocities(positions_velocity_coupling);
276
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) {
285 auto &p = *ptr;
286 auto coupling_mode = particle_force;
287#ifdef ENGINE
288 if (p.swimming().is_engine_force_on_fluid) {
289 coupling_mode = swimmer_force_on_fluid;
290 }
291#endif
292 Utils::Vector3d force_on_particle = {};
293 if (coupling_mode == particle_force) {
294 auto &v_fluid = *it_interpolated_velocities;
295 if (m_box_geo.type() == BoxType::LEES_EDWARDS) {
296 // Account for the case where the interpolated velocity has been read
297 // from a ghost of the particle across the LE boundary (or vice verssa)
298 // Then the particle velocity is shifted by +,- the LE shear velocity
299 auto const vel_correction = lees_edwards_vel_shift(
300 *it_positions_velocity_coupling, p.pos(), m_box_geo);
301 v_fluid += vel_correction;
302 }
303 auto const drag_force = lb_drag_force(p, m_thermostat.gamma, v_fluid);
304 auto const random_force = get_noise_term(p);
305 force_on_particle = drag_force + random_force;
306 ++it_interpolated_velocities;
307 ++it_positions_velocity_coupling;
308 }
309
310 auto force_on_fluid = -force_on_particle;
311#ifdef ENGINE
312 if (coupling_mode == swimmer_force_on_fluid) {
313 force_on_fluid = p.calc_director() * p.swimming().f_swim;
314 }
315#endif
316
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) {
322 /* Particle is in our LB volume, so this node
323 * is responsible to adding its force */
324 p.force() += force_on_particle;
325 }
326 force_coupling_forces.emplace_back(force_on_fluid);
327 ++it_positions_force_coupling;
328 }
329 }
330 m_lb.add_forces_at_pos(positions_force_coupling, force_coupling_forces);
331}
332
333#if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY)
334static void lb_coupling_sanity_checks(Particle const &p) {
335 /*
336 lb does (at the moment) not support rotational particle coupling.
337 Consequently, anisotropic particles are also not supported.
338 */
339 auto const &p_gamma = p.gamma();
340 if (p_gamma[0] != p_gamma[1] or p_gamma[1] != p_gamma[2]) {
341 runtimeErrorMsg() << "anisotropic particle (id " << p.id()
342 << ") coupled to LB.";
343 }
344}
345#endif
346
347} // namespace LB
348
350#ifdef CALIPER
351 CALI_CXX_MARK_FUNCTION;
352#endif
353 assert(thermostat->lb != nullptr);
354 if (thermostat->lb->couple_to_md) {
355 if (not lb.is_solver_set()) {
356 runtimeErrorMsg() << "The LB thermostat requires a LB fluid";
357 return;
358 }
359 auto const real_particles = cell_structure->local_particles();
360 auto const ghost_particles = cell_structure->ghost_particles();
361 LB::ParticleCoupling coupling{*thermostat->lb, lb, *box_geo, *local_geo};
362 LB::CouplingBookkeeping bookkeeping{*cell_structure};
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)
370#endif
371 particles.emplace_back(&p);
372 }
373 }
374 }
375 coupling.kernel(particles);
376 }
377}
@ LEES_EDWARDS
Vector implementation and trait types for boost qvm interoperability.
auto folded_position(Utils::Vector3d const &pos) const
Calculate coordinates folded to primary simulation box.
Utils::Vector3d const & length() const
Box length.
LeesEdwardsBC const & lees_edwards_bc() const
BoxType type() const
Keep track of ghost particles that have already been coupled once.
void kernel(std::vector< Particle * > const &particles)
Utils::Vector3d get_noise_term(Particle const &p) const
auto const & my_right() const
Right (top, back) corner of this nodes local box.
Definition LocalBox.hpp:47
auto const & my_left() const
Left (bottom, front) corner of this nodes local box.
Definition LocalBox.hpp:45
void lb_couple_particles()
Calculate particle-lattice interactions.
static DEVICE_QUALIFIER constexpr Vector< T, N > broadcast(typename Base::value_type const &value)
Create a vector that has all entries set to the same value.
Definition Vector.hpp:110
This file contains the defaults for ESPResSo.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
std::pair< T, I > periodic_fold(T x, I i, T const l)
Fold value into primary interval.
static void lb_coupling_sanity_checks(Particle const &p)
bool is_tracer(Particle const &p)
auto const & handle_particle_gamma(GammaType const &particle_gamma, GammaType const &default_gamma)
Vector< T, N > sqrt(Vector< T, N > const &a)
Definition Vector.hpp:357
auto hadamard_product(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:379
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)
Random number generation using Philox.
uint64_t rng_counter() const
Get current value of the RNG.
uint32_t rng_seed() const
double gamma
Friction coefficient.
Utils::Vector3d get_coupling_interpolated_velocity(Utils::Vector3d const &pos) const
Calculate the interpolated fluid velocity in MD units.
void add_forces_at_pos(std::vector< Utils::Vector3d > const &pos, std::vector< Utils::Vector3d > const &forces)
double get_agrid() const
Get the LB grid spacing.
std::vector< Utils::Vector3d > get_coupling_interpolated_velocities(std::vector< Utils::Vector3d > const &pos) const
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & gamma() const
Definition Particle.hpp:531
auto const & v() const
Definition Particle.hpp:433
auto const & mu_E() const
Definition Particle.hpp:514
auto const & id() const
Definition Particle.hpp:414