ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
particle_coupling.hpp
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
20#pragma once
21
22#include "BoxGeometry.hpp"
23#include "LocalBox.hpp"
24#include "Particle.hpp"
25#include "ParticleRange.hpp"
26#include "PropagationMode.hpp"
28#include "lb/Solver.hpp"
29#include "thermostat.hpp"
30
31#include <utils/Vector.hpp>
32#include <utils/math/sqr.hpp>
33
34#include <cmath>
35#include <unordered_set>
36#include <vector>
37
38/**
39 * @brief Check if a position is within the local LB domain plus halo.
40 *
41 * @return True iff the point is inside of the domain.
42 */
43bool in_local_halo(LocalBox const &local_box, Utils::Vector3d const &pos,
44 double agrid);
45
46// internal function exposed for unit testing
47std::vector<Utils::Vector3d> positions_in_halo(Utils::Vector3d const &pos,
48 BoxGeometry const &box_geo,
49 LocalBox const &local_box,
50 double agrid);
51
52/** @brief Calculate drag force on a single particle.
53 *
54 * See section II.C. and eq. 9 in @cite ahlrichs99a.
55 *
56 * @param[in] lb The coupled fluid
57 * @param[in] lb_gamma The friction coefficient
58 * @param[in] p The coupled particle
59 * @param[in] shifted_pos The particle position in MD units with optional shift
60 *
61 * @return The viscous coupling force
62 */
63Utils::Vector3d lb_drag_force(LB::Solver const &lb, double lb_gamma,
64 Particle const &p,
65 Utils::Vector3d const &shifted_pos);
66
67namespace LB {
68
70 LBThermostat const &m_thermostat;
71 LB::Solver &m_lb;
72 BoxGeometry const &m_box_geo;
73 LocalBox const &m_local_box;
74 double m_noise_pref_wo_gamma;
75 bool m_thermalized;
76
77public:
79 BoxGeometry const &box_geo, LocalBox const &local_box)
80 : m_thermostat{thermostat}, m_lb{lb}, m_box_geo{box_geo},
81 m_local_box{local_box} {
82 /* Eq. (16) @cite ahlrichs99a, without the gamma term.
83 * The factor 12 comes from the fact that we use random numbers
84 * from -0.5 to 0.5 (equally distributed) which have variance 1/12.
85 * The time step comes from the discretization.
86 */
87 auto constexpr variance_inv = 12.;
88 auto const kT = lb.get_kT() * Utils::sqr(lb.get_lattice_speed());
89 m_thermalized = (kT != 0.);
90 m_noise_pref_wo_gamma = std::sqrt(variance_inv * 2. * kT / lb.get_tau());
91 }
92
94 void kernel(std::vector<Particle *> const &particles);
95};
96
97/**
98 * @brief Keep track of ghost particles that have already been coupled once.
99 * In certain cases, there may be more than one ghost for the same particle.
100 * To make sure these are only coupled once, ghosts' ids are recorded.
101 */
103 std::unordered_set<int> m_coupled_ghosts;
104 CellStructure const &m_cell_structure;
105
106 /** @brief Check if there is locally a real particle for the given ghost. */
107 bool is_ghost_for_local_particle(Particle const &p) const {
108 return not m_cell_structure.get_local_particle(p.id())->is_ghost();
109 }
110
111public:
112 explicit CouplingBookkeeping(CellStructure const &cell_structure)
113 : m_cell_structure{cell_structure} {}
114
115 /** @brief Determine if a given particle should be coupled. */
117 auto const propagation = p.propagation();
118 if ((propagation & PropagationMode::TRANS_LB_MOMENTUM_EXCHANGE) == 0 and
119 (propagation & PropagationMode::SYSTEM_DEFAULT) == 0 and
120 (propagation & PropagationMode::TRANS_LB_TRACER) == 0) {
121 return false;
122 }
123#ifdef VIRTUAL_SITES_RELATIVE
124 if ((propagation & PropagationMode::TRANS_LB_MOMENTUM_EXCHANGE) == 0 and
127 return false;
128 }
129#endif
130 // real particles: always couple
131 if (not p.is_ghost()) {
132 return true;
133 }
134 // ghosts: check we don't have the corresponding real particle on the same
135 // node, and that a ghost for the same particle hasn't been coupled already
136 if (m_coupled_ghosts.count(p.id()) != 0 or is_ghost_for_local_particle(p)) {
137 return false;
138 }
139 m_coupled_ghosts.insert(p.id());
140 return true;
141 }
142};
143
144inline bool is_tracer(Particle const &p) {
146}
147
148} // namespace LB
Vector implementation and trait types for boost qvm interoperability.
Keep track of ghost particles that have already been coupled once.
CouplingBookkeeping(CellStructure const &cell_structure)
bool should_be_coupled(Particle const &p)
Determine if a given particle should be coupled.
void kernel(std::vector< Particle * > const &particles)
ParticleCoupling(LBThermostat const &thermostat, LB::Solver &lb, BoxGeometry const &box_geo, LocalBox const &local_box)
Utils::Vector3d get_noise_term(Particle const &p) const
bool is_tracer(Particle const &p)
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
Utils::Vector3d lb_drag_force(LB::Solver const &lb, double lb_gamma, Particle const &p, Utils::Vector3d const &shifted_pos)
Calculate drag force on a single particle.
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.
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.
Describes a cell structure / cell system.
Particle * get_local_particle(int id)
Get a local particle by id.
Thermostat for lattice-Boltzmann particle coupling.
double get_tau() const
Get the LB time step.
auto get_lattice_speed() const
Get the lattice speed (agrid/tau).
double get_kT() const
Get the thermal energy parameter of the LB fluid.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & propagation() const
Definition Particle.hpp:421
bool is_ghost() const
Definition Particle.hpp:440
auto const & id() const
Definition Particle.hpp:414