ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
ImmersedBoundaries.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
22#include "BoxGeometry.hpp"
23#include "Particle.hpp"
25#include "communication.hpp"
26#include "ibm_volcons.hpp"
27#include "system/System.hpp"
28
30
31#include <boost/mpi/collectives/all_reduce.hpp>
32#include <boost/range/algorithm/find_if.hpp>
33
34#include <algorithm>
35#include <functional>
36#include <span>
37#include <utility>
38#include <vector>
39
40/** Calculate volumes, volume force and add it to each virtual particle. */
42 if (VolumeInitDone && !BoundariesFound) {
43 return;
44 }
45 calc_volumes(cs);
46 calc_volume_force(cs);
47}
48
49/** Initialize volume conservation */
51 auto const &bonded_ias = *get_system().bonded_ias;
52 // Check since this function is called at the start of every integrate loop
53 // Also check if volume has been set due to reading of a checkpoint
54 if (not BoundariesFound) {
55 BoundariesFound = std::ranges::any_of(bonded_ias, [](auto const &kv) {
56 return (boost::get<IBMVolCons>(&(*kv.second)) != nullptr);
57 });
58 }
59
60 if (!VolumeInitDone && BoundariesFound) {
61 // Calculate volumes
62 calc_volumes(cs);
63
64 // Loop through all bonded interactions and check if we need to set the
65 // reference volume
66 for (auto &kv : bonded_ias) {
67 if (auto *v = boost::get<IBMVolCons>(kv.second.get())) {
68 // This check is important because InitVolumeConservation may be called
69 // accidentally during the integration. Then we must not reset the
70 // reference
71 BoundariesFound = true;
72 if (v->volRef == 0.) {
73 v->volRef = VolumesCurrent[v->softID];
74 }
75 }
76 }
77
78 VolumeInitDone = true;
79 }
80}
81
82static IBMVolCons const *
84 Particle const &p1) {
85 auto const it = boost::find_if(p1.bonds(), [&](auto const &bond) -> bool {
86 return boost::get<IBMVolCons>(bonded_ias.at(bond.bond_id()).get());
87 });
88
89 return (it != p1.bonds().end())
90 ? boost::get<IBMVolCons>(bonded_ias.at(it->bond_id()).get())
91 : nullptr;
92}
93
94/** Calculate partial volumes on all compute nodes and call MPI to sum up.
95 * See @cite zhang01b, @cite dupin08a, @cite kruger12a.
96 */
97void ImmersedBoundaries::calc_volumes(CellStructure &cs) {
98
99 if (!BoundariesFound)
100 return;
101
102 auto const &box_geo = *get_system().box_geo;
103 auto const &bonded_ias = *get_system().bonded_ias;
104
105 // Partial volumes for each soft particle, to be summed up
106 std::vector<double> tempVol(VolumesCurrent.size());
107
108 // Loop over all particles on local node
109 cs.bond_loop([&tempVol, &box_geo, &bonded_ias](
110 Particle &p1, int bond_id, std::span<Particle *> partners) {
111 auto const vol_cons_params = vol_cons_parameters(bonded_ias, p1);
112
113 if (vol_cons_params &&
114 boost::get<IBMTriel>(bonded_ias.at(bond_id).get()) != nullptr) {
115 // Our particle is the leading particle of a triel
116 // Get second and third particle of the triangle
117 Particle &p2 = *partners[0];
118 Particle &p3 = *partners[1];
119
120 // Unfold position of first node.
121 // This is to get a continuous trajectory with no jumps when box
122 // boundaries are crossed.
123 auto const x1 = box_geo.unfolded_position(p1.pos(), p1.image_box());
124 auto const x2 = x1 + box_geo.get_mi_vector(p2.pos(), x1);
125 auto const x3 = x1 + box_geo.get_mi_vector(p3.pos(), x1);
126
127 // Volume of this tetrahedron
128 // See @cite zhang01b
129 // The volume can be negative, but it is not necessarily the
130 // "signed volume" in the above paper (the sign of the real
131 // "signed volume" must be calculated using the normal vector; the
132 // result of the calculation here is simply a term in the sum
133 // required to calculate the volume of a particle). Again, see the
134 // paper. This should be equivalent to the formulation using
135 // vector identities in @cite kruger12a
136
137 const double v321 = x3[0] * x2[1] * x1[2];
138 const double v231 = x2[0] * x3[1] * x1[2];
139 const double v312 = x3[0] * x1[1] * x2[2];
140 const double v132 = x1[0] * x3[1] * x2[2];
141 const double v213 = x2[0] * x1[1] * x3[2];
142 const double v123 = x1[0] * x2[1] * x3[2];
143
144 tempVol[vol_cons_params->softID] +=
145 1.0 / 6.0 * (-v321 + v231 + v312 - v132 - v213 + v123);
146 }
147 return false;
148 });
149
150 // Sum up and communicate
151 boost::mpi::all_reduce(comm_cart, tempVol.data(),
152 static_cast<int>(tempVol.size()),
153 VolumesCurrent.data(), std::plus<double>());
154}
155
156/** Calculate and add the volume force to each node */
157void ImmersedBoundaries::calc_volume_force(CellStructure &cs) {
158 if (!BoundariesFound)
159 return;
160
161 auto const &box_geo = *get_system().box_geo;
162 auto const &bonded_ias = *get_system().bonded_ias;
163
164 cs.bond_loop([this, &box_geo, &bonded_ias](Particle &p1, int bond_id,
165 std::span<Particle *> partners) {
166 if (boost::get<IBMTriel>(bonded_ias.at(bond_id).get()) != nullptr) {
167 // Check if particle has an IBM Triel bonded interaction and an
168 // IBM VolCons bonded interaction. Basically this loops over all
169 // triangles, not all particles. First round to check for volume
170 // conservation.
171 auto const vol_cons_params = vol_cons_parameters(bonded_ias, p1);
172 if (not vol_cons_params)
173 return false;
174
175 auto const current_volume =
176 VolumesCurrent[static_cast<unsigned int>(vol_cons_params->softID)];
177
178 // Our particle is the leading particle of a triel
179 // Get second and third particle of the triangle
180 Particle &p2 = *partners[0];
181 Particle &p3 = *partners[1];
182
183 // Unfold position of first node.
184 // This is to get a continuous trajectory with no jumps when box
185 // boundaries are crossed.
186 auto const x1 = box_geo.unfolded_position(p1.pos(), p1.image_box());
187
188 // Unfolding seems to work only for the first particle of a triel
189 // so get the others from relative vectors considering PBC
190 auto const a12 = box_geo.get_mi_vector(p2.pos(), x1);
191 auto const a13 = box_geo.get_mi_vector(p3.pos(), x1);
192
193 // Now we have the true and good coordinates
194 // This is eq. (9) in @cite dupin08a.
195 auto const n = vector_product(a12, a13);
196 const double ln = n.norm();
197 const double A = 0.5 * ln;
198 const double fact = vol_cons_params->kappaV *
199 (current_volume - vol_cons_params->volRef) /
200 current_volume;
201
202 auto const nHat = n / ln;
203 auto const force = -fact * A * nHat;
204
205 p1.force() += force;
206 p2.force() += force;
207 p3.force() += force;
208 }
209 return false;
210 });
211}
212
214 auto const new_size = bond.softID + 1u;
215 if (new_size > VolumesCurrent.size()) {
216 VolumesCurrent.resize(new_size);
217 }
218 bond.set_volumes_view(VolumesCurrent);
219}
static IBMVolCons const * vol_cons_parameters(BondedInteractionsMap const &bonded_ias, Particle const &p1)
Data structures for bonded interactions.
container for bonded interactions.
mapped_type at(key_type const &key) const
void volume_conservation(CellStructure &cs)
Calculate volumes, volume force and add it to each virtual particle.
void init_volume_conservation(CellStructure &cs)
Initialize volume conservation.
void register_softID(IBMVolCons &bond)
boost::mpi::communicator comm_cart
The communicator.
Describes a cell structure / cell system.
void bond_loop(BondKernel const &bond_kernel)
Bonded pair loop.
Parameters for IBM volume conservation bond.
void set_volumes_view(std::vector< double > const &volumes)
unsigned int softID
ID of the large soft particle to which this node belongs.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & image_box() const
Definition Particle.hpp:444
auto const & bonds() const
Definition Particle.hpp:428
auto const & pos() const
Definition Particle.hpp:431