32#include "communication.hpp"
34#include "system/System.hpp"
41#include <boost/mpi/collectives/broadcast.hpp>
42#include <boost/mpi/collectives/reduce.hpp>
57template <
typename Archive,
typename... T>
58void serialize(Archive &ar, std::tuple<T...> &pack,
unsigned int const) {
59 std::apply([&](
auto &...item) { ((ar & item), ...); }, pack);
65 using type = std::tuple<std::invoke_result_t<F, Particle const &>...>;
69 using type = std::invoke_result_t<F, Particle const &>;
77template <
class... Trait>
79 std::vector<int>
const &p_types,
85 buffer.emplace_back(trait(p)...);
97 std::vector<int>
const &set2) {
100 std::vector<int> buf_type{};
101 std::vector<Utils::Vector3d> buf_pos{};
103 auto const &box_geo = *system.
box_geo;
104 auto const accept_all = set1.empty() or set2.empty();
106 if (accept_all or contains(set1, p.type()) or contains(set2, p.type())) {
107 buf_type.emplace_back(p.type());
108 buf_pos.emplace_back(box_geo.unfolded_position(p.pos(), p.image_box()));
119 auto mindist_sq = std::numeric_limits<double>::infinity();
120 for (std::size_t j = 0ul; j < buf_type.size(); ++j) {
123 if (set1.empty() || contains(set1, buf_type[j]))
125 if (set2.empty() || contains(set2, buf_type[j]))
129 for (
auto i = j + 1ul; i != buf_type.size(); ++i) {
132 if (((in_set & 1) and (set2.empty() or contains(set2, buf_type[i]))) or
133 ((in_set & 2) and (set1.empty() or contains(set1, buf_type[i])))) {
134 mindist_sq = std::min(
135 mindist_sq, box_geo.get_mi_vector(buf_pos[j], buf_pos[i]).norm2());
140 return std::sqrt(mindist_sq);
144 bool include_particles,
145 bool include_lbfluid) {
147 if (include_particles) {
152 return m + p.mass() * p.v();
162 auto const &box_geo = *system.
box_geo;
165 double local_mass = 0.;
167 for (
auto const &p : cell_structure.local_particles()) {
168 if ((p.type() == p_type or p_type == -1) and not p.is_virtual()) {
169 local_com += box_geo.unfolded_position(p.pos(), p.image_box()) * p.mass();
170 local_mass += p.mass();
175 boost::mpi::reduce(
::comm_cart, local_com, com, std::plus<>(), 0);
176 boost::mpi::reduce(
::comm_cart, local_mass, mass, std::plus<>(), 0);
181 auto const &box_geo = *system.
box_geo;
185 for (
auto const &p : cell_structure.local_particles()) {
186 if ((p.type() == p_type or p_type == -1) and not p.is_virtual()) {
187 auto const pos = box_geo.unfolded_position(p.pos(), p.image_box());
195 std::vector<int>
const &p_types) {
196 auto const &box_geo = *system.
box_geo;
197 auto const trait_pos = [&box_geo](
Particle const &p) {
198 return box_geo.unfolded_position(p.pos(), p.image_box());
206 static_cast<double>(buf_pos.size());
208 for (
unsigned int i = 0u; i < 3u; ++i) {
209 for (
unsigned int j = 0u; j < 3u; ++j) {
211 mat[i * 3u + j] = mat[j * 3u + i];
213 mat[i * 3u + j] = std::accumulate(
214 buf_pos.begin(), buf_pos.end(), 0.,
216 return acc + (pos[i] - center[i]) * (pos[j] - center[j]);
221 mat /=
static_cast<double>(buf_pos.size());
228 auto const &box_geo = *system.
box_geo;
234 for (
auto const &p : cell_structure.local_particles()) {
235 if (p.type() == p_type and not p.is_virtual()) {
236 auto const pos = box_geo.unfolded_position(p.pos(), p.image_box()) - com;
237 auto const mass = p.mass();
238 mat[0] += mass * (pos[1] * pos[1] + pos[2] * pos[2]);
239 mat[4] += mass * (pos[0] * pos[0] + pos[2] * pos[2]);
240 mat[8] += mass * (pos[0] * pos[0] + pos[1] * pos[1]);
241 mat[1] -= mass * (pos[0] * pos[1]);
242 mat[2] -= mass * (pos[0] * pos[2]);
243 mat[5] -= mass * (pos[1] * pos[2]);
255 std::vector<int> buf_pid{};
256 auto const dist_sq = dist * dist;
257 auto const &box_geo = *system.
box_geo;
260 auto const r_sq = box_geo.get_mi_vector(pos, p.pos()).norm2();
261 if (r_sq < dist_sq) {
262 buf_pid.push_back(p.id());
274std::vector<std::vector<double>>
276 std::vector<int>
const &p1_types,
277 std::vector<int>
const &p2_types,
double r_min,
278 double r_max,
int r_bins,
bool log_flag,
bool int_flag) {
280 auto const &box_geo = *system.
box_geo;
281 auto const trait_id = [](
Particle const &p) {
return p.id(); };
282 auto const trait_pos = [&box_geo](
Particle const &p) {
283 return box_geo.unfolded_position(p.pos(), p.image_box());
291 auto const start_dist2 =
Utils::sqr(r_max + 1.);
292 auto const inv_bin_width =
293 (log_flag) ?
static_cast<double>(r_bins) / std::log(r_max / r_min)
294 :
static_cast<double>(r_bins) / (r_max - r_min);
298 std::vector<double> distribution(r_bins);
300 for (
auto const &[pid1, pos1] : buf1) {
301 auto min_dist2 = start_dist2;
303 for (
auto const &[pid2, pos2] : buf2) {
305 auto const act_dist2 = box_geo.get_mi_vector(pos1, pos2).norm2();
306 if (act_dist2 < min_dist2) {
307 min_dist2 = act_dist2;
311 if (min_dist2 <= r_max2) {
312 if (min_dist2 >= r_min2) {
313 auto const min_dist = std::sqrt(min_dist2);
315 auto const ind =
static_cast<int>(
316 ((log_flag) ? std::log(min_dist / r_min) : (min_dist - r_min)) *
318 if (ind >= 0 and ind < r_bins) {
319 distribution[ind] += 1.0;
330 low /=
static_cast<double>(cnt);
331 for (
int i = 0; i < r_bins; i++) {
332 distribution[i] /=
static_cast<double>(cnt);
337 distribution[0] += low;
338 for (
int i = 0; i < r_bins - 1; i++)
339 distribution[i + 1] += distribution[i];
343 std::vector<double> radii(r_bins);
345 auto const log_fac = std::pow(r_max / r_min, 1. / r_bins);
346 radii[0] = r_min * std::sqrt(log_fac);
347 for (
int i = 1; i < r_bins; ++i) {
348 radii[i] = radii[i - 1] * log_fac;
351 auto const bin_width = (r_max - r_min) /
static_cast<double>(r_bins);
352 for (
int i = 0; i < r_bins; ++i) {
353 radii[i] = r_min + bin_width / 2. +
static_cast<double>(i) * bin_width;
357 return {radii, distribution};
360std::vector<std::vector<double>>
363 auto const &box_geo = *system.
box_geo;
364 auto const trait_pos = [&box_geo](
Particle const &p) {
365 return box_geo.unfolded_position(p.pos(), p.image_box());
368 auto const order_sq =
Utils::sqr(
static_cast<std::size_t
>(order));
369 auto const twoPI_L = 2. * std::numbers::pi * system.
box_geo->length_inv()[0];
370 std::vector<double> ff(2ul * order_sq + 1ul);
371 std::vector<double> wavevectors;
372 std::vector<double> intensities;
375 for (
int i = 0; i <= order; i++) {
376 for (
int j = -order; j <= order; j++) {
377 for (
int k = -order; k <= order; k++) {
378 auto const n = i * i + j * j + k * k;
379 if ((
static_cast<std::size_t
>(n) <= order_sq) && (n >= 1)) {
380 double C_sum = 0.0, S_sum = 0.0;
381 for (
auto const &pos : buf_pos) {
386 ff[2 * n - 2] += C_sum * C_sum + S_sum * S_sum;
393 std::size_t length = 0;
394 for (std::size_t qi = 0; qi < order_sq; qi++) {
395 if (ff[2 * qi + 1] != 0) {
396 ff[2 * qi] /=
static_cast<double>(buf_pos.size()) * ff[2 * qi + 1];
401 wavevectors.resize(length);
402 intensities.resize(length);
405 for (std::size_t i = 0; i < order_sq; i++) {
406 if (ff[2 * i + 1] != 0) {
407 wavevectors[cnt] = twoPI_L * std::sqrt(
static_cast<double>(i + 1));
408 intensities[cnt] = ff[2 * i];
414 return {std::move(wavevectors), std::move(intensities)};
Vector implementation and trait types for boost qvm interoperability.
std::shared_ptr< CellStructure > cell_structure
std::shared_ptr< BoxGeometry > box_geo
boost::mpi::communicator comm_cart
The communicator.
__device__ void vector_product(float const *a, float const *b, float *out)
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
void gather_buffer(std::vector< T, Allocator > &buffer, boost::mpi::communicator const &comm, int root=0)
Gather buffer with different size on each node.
bool contains(InputIt first, InputIt last, T const &value)
Check whether an iterator range contains a value.
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
void serialize(Archive &ar, std::tuple< T... > &pack, unsigned int const)
Serialize std::tuple.
Utils::Vector3d center_of_mass(System::System const &system, int p_type)
Calculate the center of mass of particles of a certain type.
Utils::Vector3d angular_momentum(System::System const &system, int p_type)
Calculate the angular momentum of particles of a certain type.
std::vector< int > nbhood(System::System const &system, Utils::Vector3d const &pos, double dist)
Find all particles within a given radius dist around a position pos.
Utils::Vector9d moment_of_inertia_matrix(System::System const &system, int p_type)
Calculate the moment of inertia of particles of a certain type.
Utils::Vector9d gyration_tensor(System::System const &system, std::vector< int > const &p_types)
Calculate the gyration tensor of particles of certain types.
Utils::Vector3d calc_linear_momentum(System::System const &system, bool include_particles, bool include_lbfluid)
Calculate total momentum of the system (particles & LB fluid).
std::vector< std::vector< double > > structure_factor(System::System const &system, std::vector< int > const &p_types, int order)
Calculate the spherically averaged structure factor.
double mindist(System::System const &system, std::vector< int > const &set1, std::vector< int > const &set2)
Calculate the minimal distance of two particles with types in set1 and set2, respectively.
static auto gather_traits_for_types(System::System const &system, std::vector< int > const &p_types, Trait &&...trait)
Gather particle traits to MPI rank 0.
std::vector< std::vector< double > > calc_part_distribution(System::System const &system, std::vector< int > const &p1_types, std::vector< int > const &p2_types, double r_min, double r_max, int r_bins, bool log_flag, bool int_flag)
Calculate the distribution of particles around others.
Statistical tools to analyze simulations.
std::invoke_result_t< F, Particle const & > type
Decay a tuple of only 1 type to that type.
std::tuple< std::invoke_result_t< F, Particle const & >... > type
bool is_solver_set() const
Return true if a LB solver is active.
auto get_lattice_speed() const
Get the lattice speed (agrid/tau).
Utils::Vector3d get_momentum() const
Struct holding all information for one particle.