From 9f3a9bdb9d44754c69e1fe2be2b7903e904b9ec1 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 24 Oct 2019 21:41:29 +0200 Subject: [PATCH 1/6] core: engine: Removed shear torque calculation --- src/core/Particle.hpp | 6 +--- .../lb_particle_coupling.cpp | 15 -------- src/core/rotation.cpp | 36 ------------------- src/python/espressomd/particle_data.pxd | 3 -- src/python/espressomd/particle_data.pyx | 22 ++---------- testsuite/python/engine_lb.py | 8 ++--- 6 files changed, 8 insertions(+), 82 deletions(-) diff --git a/src/core/Particle.hpp b/src/core/Particle.hpp index 3cb3c9271ba..eb58792621b 100644 --- a/src/core/Particle.hpp +++ b/src/core/Particle.hpp @@ -248,15 +248,11 @@ struct ParticleParametersSwimming { double v_swim = 0.; int push_pull = 0; double dipole_length = 0.; - Utils::Vector3d v_center; - Utils::Vector3d v_source; - double rotational_friction = 0.; #endif template void serialize(Archive &ar, long int) { #ifdef ENGINE - ar &swimming &f_swim &v_swim &push_pull &dipole_length &v_center &v_source - &rotational_friction; + ar &swimming &f_swim &v_swim &push_pull &dipole_length; #endif } }; diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index 9a935f0096f..9ada55db348 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.cpp @@ -223,13 +223,6 @@ bool in_local_halo(Vector3d const &pos) { #ifdef ENGINE void add_swimmer_force(Particle &p) { if (p.swim.swimming) { - if (in_local_domain(p.r.p, local_geo)) { - p.swim.v_center = lb_lbinterpolation_get_interpolated_velocity(p.r.p) * - lb_lbfluid_get_lattice_speed(); - } else { - p.swim.v_center = {}; - } - // calculate source position const double direction = double(p.swim.push_pull) * p.swim.dipole_length; auto const director = p.r.calc_director(); @@ -239,14 +232,6 @@ void add_swimmer_force(Particle &p) { return; } - if (in_local_domain(source_position, local_geo)) { - p.swim.v_source = - lb_lbinterpolation_get_interpolated_velocity(source_position) * - lb_lbfluid_get_lattice_speed(); - } else { - p.swim.v_source = {}; - } - add_md_force(source_position, p.swim.f_swim * director); } } diff --git a/src/core/rotation.cpp b/src/core/rotation.cpp index f2dc1c46112..b480f2fd0a8 100644 --- a/src/core/rotation.cpp +++ b/src/core/rotation.cpp @@ -175,14 +175,6 @@ inline void convert_torque_to_body_frame_apply_fix(Particle &p) { /** convert the torques to the body-fixed frames and propagate angular * velocities */ void convert_torques_propagate_omega(const ParticleRange &particles) { - -#if defined(CUDA) && defined(ENGINE) - if ((lb_lbfluid_get_lattice_switch() == ActiveLB::GPU) && - swimming_particles_exist) { - copy_v_cs_from_GPU(particles); - } -#endif - for (auto &p : particles) { // Skip particle if rotation is turned off entirely for it. if (!p.p.rotation) @@ -190,34 +182,6 @@ void convert_torques_propagate_omega(const ParticleRange &particles) { convert_torque_to_body_frame_apply_fix(p); -#if defined(ENGINE) - if (p.swim.swimming && lb_lbfluid_get_lattice_switch() != ActiveLB::NONE) { - if (lb_lbfluid_get_lattice_switch() == ActiveLB::CPU && n_nodes > 1 && - p.swim.rotational_friction != 0.) { - runtimeErrorMsg() << "ENGINE rotational_friction feature with CPU-LB " - "only implemented for one CPU core"; - } - - auto const dip = p.swim.dipole_length * p.r.calc_director(); - - auto const diff = p.swim.v_center - p.swim.v_source; - - const Utils::Vector3d cross = vector_product(diff, dip); - const double l_diff = diff.norm(); - const double l_cross = cross.norm(); - - if (l_cross > 0 && p.swim.dipole_length > 0) { - auto const omega_swim = - l_diff / (l_cross * p.swim.dipole_length) * cross; - - auto const omega_swim_body = - convert_vector_space_to_body(p, omega_swim); - p.f.torque += - p.swim.rotational_friction * (omega_swim_body - p.m.omega); - } - } -#endif - // Propagation of angular velocities p.m.omega += hadamard_division(time_step_half * p.f.torque, p.p.rinertia); diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd index 8b501b2a2fa..48f052b54e6 100644 --- a/src/python/espressomd/particle_data.pxd +++ b/src/python/espressomd/particle_data.pxd @@ -76,9 +76,6 @@ cdef extern from "particle_data.hpp": double v_swim int push_pull double dipole_length - double v_center[3] - double v_source[3] - double rotational_friction # Setter/getter/modifier functions functions void prefetch_particle_data(vector[int] ids) diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx index aec789bbeaf..7fdfd1903d9 100644 --- a/src/python/espressomd/particle_data.pyx +++ b/src/python/espressomd/particle_data.pyx @@ -1176,19 +1176,11 @@ cdef class ParticleHandle: dipole_length : :obj:`float` This determines the distance of the source of propulsion from the particle's center. - rotational_friction : :obj:`float` - This key can be used to set the friction that causes - the orientation of the particle to change in shear - flow. The torque on the particle is determined by - taking the cross product of the difference between the - fluid velocity at the center of the particle and at - the source point and the vector connecting the center - and source. Notes ----- This needs the feature ``ENGINE``. The keys ``'mode'``, - ``'dipole_length'``, and ``'rotational_friction'`` are only + and ``'dipole_length'`` are only available if ``ENGINE`` is used with LB or ``CUDA``. Examples @@ -1202,7 +1194,7 @@ cdef class ParticleHandle: >>> >>> # Usage with LB >>> system.part.add(id=1, pos=[2,0,0], swimming={'f_swim': 0.01, - ... 'mode': 'pusher', 'dipole_length': 2.0, 'rotational_friction': 20}) + ... 'mode': 'pusher', 'dipole_length': 2.0}) """ @@ -1214,7 +1206,6 @@ cdef class ParticleHandle: swim.f_swim = 0.0 swim.push_pull = 0 swim.dipole_length = 0.0 - swim.rotational_friction = 0.0 if type(_params) == type(True): if _params: @@ -1252,12 +1243,6 @@ cdef class ParticleHandle: _params['dipole_length'], 1, float, "dipole_length has to be a float.") swim.dipole_length = _params['dipole_length'] - if 'rotational_friction' in _params: - check_type_or_throw_except( - _params['rotational_friction'], 1, float, "rotational_friction has to be a float.") - swim.rotational_friction = _params[ - 'rotational_friction'] - if swim.f_swim != 0 or swim.v_swim != 0: swimming_particles_exist = True mpi_bcast_parameter(FIELD_SWIMMING_PARTICLES_EXIST) @@ -1279,8 +1264,7 @@ cdef class ParticleHandle: 'v_swim': _swim.v_swim, 'f_swim': _swim.f_swim, 'mode': mode, - 'dipole_length': _swim.dipole_length, - 'rotational_friction': _swim.rotational_friction + 'dipole_length': _swim.dipole_length } return swim diff --git a/testsuite/python/engine_lb.py b/testsuite/python/engine_lb.py index af30d7a663a..f7a8ac2dfae 100644 --- a/testsuite/python/engine_lb.py +++ b/testsuite/python/engine_lb.py @@ -54,19 +54,19 @@ def add_all_types_of_swimmers( system.part.add(pos=pos0, quat=minus_y, fix=3 * [fix], mass=0.9, rinertia=3 * [7], rotation=3 * [rotation], swimming={"mode": "pusher", "f_swim": 0.10, - "dipole_length": 0.5, "rotational_friction": 0.3}) + "dipole_length": 0.5}) system.part.add(pos=pos1, quat=plus_x, fix=3 * [fix], mass=1.9, rinertia=3 * [8], rotation=3 * [rotation], swimming={"mode": "pusher", "v_swim": 0.02, - "dipole_length": 0.6, "rotational_friction": 0.4}) + "dipole_length": 0.6}) system.part.add(pos=pos2, quat=plus_z, fix=3 * [fix], mass=2.9, rinertia=3 * [9], rotation=3 * [rotation], swimming={"mode": "puller", "f_swim": 0.08, - "dipole_length": 0.7, "rotational_friction": 0.8}) + "dipole_length": 0.7}) system.part.add(pos=pos3, quat=plus_y, fix=3 * [fix], mass=3.9, rinertia=3 * [10], rotation=3 * [rotation], swimming={"mode": "puller", "v_swim": 0.05, - "dipole_length": 0.8, "rotational_friction": 0.3}) + "dipole_length": 0.8}) def tearDown(self): self.system.part.clear() From a1e1925d7b8140351d6f771336cbee714ad63151 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 24 Oct 2019 22:21:47 +0200 Subject: [PATCH 2/6] core: Moved Particle::swim into ParticleProperties and removed custom comm code. --- src/core/Particle.hpp | 40 ++++++--------- src/core/communication.cpp | 1 - src/core/cuda_common_cuda.cu | 26 ---------- src/core/cuda_interface.cpp | 42 ++------------- src/core/cuda_interface.hpp | 18 ------- src/core/forces_inline.hpp | 4 +- src/core/ghosts.cpp | 21 -------- src/core/ghosts.hpp | 4 +- .../lb_particle_coupling.cpp | 15 ++---- src/core/particle_data.cpp | 25 ++------- .../ParticleParametersSwimming.hpp | 51 ------------------- src/core/thermostat.hpp | 4 +- 12 files changed, 36 insertions(+), 215 deletions(-) delete mode 100644 src/core/serialization/ParticleParametersSwimming.hpp diff --git a/src/core/Particle.hpp b/src/core/Particle.hpp index eb58792621b..31ec38c001c 100644 --- a/src/core/Particle.hpp +++ b/src/core/Particle.hpp @@ -34,6 +34,18 @@ enum : uint8_t { ROTATION_Z = 4u }; +struct ParticleParametersSwimming { + bool swimming = false; + double f_swim = 0.; + double v_swim = 0.; + int push_pull = 0; + double dipole_length = 0.; + + template void serialize(Archive &ar, long int /* version */) { + ar &swimming &f_swim &v_swim &push_pull &dipole_length; + } +}; + /** Properties of a particle which are not supposed to * change during the integration, but have to be known * for all ghosts. Ghosts are particles which are @@ -156,6 +168,10 @@ struct ParticleProperties { Utils::Vector3d ext_torque = {0, 0, 0}; #endif #endif + +#ifdef ENGINE + ParticleParametersSwimming swim; +#endif }; /** Positional information on a particle. Information that is @@ -240,23 +256,6 @@ struct ParticleLocal { Utils::Vector3i i = {0, 0, 0}; }; -struct ParticleParametersSwimming { -// ifdef inside because we need this type for some MPI prototypes -#ifdef ENGINE - bool swimming = false; - double f_swim = 0.; - double v_swim = 0.; - int push_pull = 0; - double dipole_length = 0.; -#endif - - template void serialize(Archive &ar, long int) { -#ifdef ENGINE - ar &swimming &f_swim &v_swim &push_pull &dipole_length; -#endif - } -}; - /** Struct holding all information for one particle. */ struct Particle { int &identity() { return p.identity; } @@ -288,9 +287,6 @@ struct Particle { ret.m = m; ret.f = f; ret.l = l; -#ifdef ENGINE - ret.swim = swim; -#endif return ret; } @@ -342,10 +338,6 @@ struct Particle { */ IntList el; #endif - -#ifdef ENGINE - ParticleParametersSwimming swim; -#endif }; #endif diff --git a/src/core/communication.cpp b/src/core/communication.cpp index 6eadbf54039..5b76f60b0b2 100644 --- a/src/core/communication.cpp +++ b/src/core/communication.cpp @@ -67,7 +67,6 @@ #include "serialization/IA_parameters.hpp" #include "serialization/Particle.hpp" -#include "serialization/ParticleParametersSwimming.hpp" #include #include diff --git a/src/core/cuda_common_cuda.cu b/src/core/cuda_common_cuda.cu index bc29fde74f6..3a61f07c742 100644 --- a/src/core/cuda_common_cuda.cu +++ b/src/core/cuda_common_cuda.cu @@ -52,9 +52,6 @@ std::vector particle_forces_host; CUDA_energy energy_host; std::vector particle_torques_host; -#ifdef ENGINE -std::vector host_v_cs; -#endif /**cuda streams for parallel computing on cpu and gpu */ cudaStream_t stream[1]; @@ -187,9 +184,6 @@ void gpu_change_number_of_part_to_comm() { particle_data_device = nullptr; } -#ifdef ENGINE - host_v_cs.clear(); -#endif #ifdef ROTATION particle_torques_host.clear(); #endif @@ -210,9 +204,6 @@ void gpu_change_number_of_part_to_comm() { cudaHostAllocWriteCombined)); particle_forces_host.resize(3 * global_part_vars_host.number_of_particles); -#ifdef ENGINE - host_v_cs.resize(global_part_vars_host.number_of_particles); -#endif #if (defined DIPOLES || defined ROTATION) particle_torques_host.resize(3 * global_part_vars_host.number_of_particles); @@ -352,23 +343,6 @@ void copy_forces_from_GPU(ParticleRange particles) { } } -#if defined(ENGINE) && defined(CUDA) -// setup and call kernel to copy v_cs to host -void copy_v_cs_from_GPU(ParticleRange particles) { - if (global_part_vars_host.communication_enabled == 1 && - global_part_vars_host.number_of_particles) { - // Copy result from device memory to host memory - if (this_node == 0) { - cuda_safe_mem(cudaMemcpy2D( - host_v_cs.data(), sizeof(CUDA_v_cs), particle_data_device, - sizeof(CUDA_particle_data), sizeof(CUDA_v_cs), - global_part_vars_host.number_of_particles, cudaMemcpyDeviceToHost)); - } - cuda_mpi_send_v_cs(particles, host_v_cs); - } -} -#endif - void clear_energy_on_GPU() { if (!global_part_vars_host.communication_enabled) // || !global_part_vars_host.number_of_particles ) diff --git a/src/core/cuda_interface.cpp b/src/core/cuda_interface.cpp index d6ebb5b90d0..e07f68d85cc 100644 --- a/src/core/cuda_interface.cpp +++ b/src/core/cuda_interface.cpp @@ -77,13 +77,13 @@ static void pack_particles(ParticleRange particles, #endif #ifdef ENGINE - buffer[i].swim.v_swim = static_cast(part.swim.v_swim); - buffer[i].swim.f_swim = static_cast(part.swim.f_swim); + buffer[i].swim.v_swim = static_cast(part.p.swim.v_swim); + buffer[i].swim.f_swim = static_cast(part.p.swim.f_swim); buffer[i].swim.director = buffer[i].director; - buffer[i].swim.push_pull = part.swim.push_pull; - buffer[i].swim.dipole_length = static_cast(part.swim.dipole_length); - buffer[i].swim.swimming = part.swim.swimming; + buffer[i].swim.push_pull = part.p.swim.push_pull; + buffer[i].swim.dipole_length = static_cast(part.p.swim.dipole_length); + buffer[i].swim.swimming = part.p.swim.swimming; #endif i++; } @@ -162,38 +162,6 @@ void cuda_mpi_send_forces(ParticleRange particles, } } -#if defined(ENGINE) && defined(CUDA) -namespace { -void set_v_cs(ParticleRange particles, const std::vector &v_cs) { - int ind = 0; - for (auto &p : particles) { - for (int i = 0; i < 3; i++) { - p.swim.v_center[i] = v_cs[ind].v_cs[0 + i]; - p.swim.v_source[i] = v_cs[ind].v_cs[3 + i]; - } - ind++; - } -} -} // namespace - -void cuda_mpi_send_v_cs(ParticleRange particles, - std::vector host_v_cs) { - // first collect number of particles on each node - auto const n_part = particles.size(); - - // call slave functions to provide the slave data - if (this_node > 0) { - std::vector buffer(n_part); - - Utils::Mpi::scatter_buffer(buffer.data(), n_part, comm_cart); - set_v_cs(particles, buffer); - } else { - Utils::Mpi::scatter_buffer(host_v_cs.data(), n_part, comm_cart); - set_v_cs(particles, host_v_cs); - } -} -#endif // ifdef ENGINE - /** Takes a CUDA_energy struct and adds it to the core energy struct. This cannot be done from inside cuda_common_cuda.cu:copy_energy_from_GPU() because energy.hpp indirectly includes on mpi.h while .cu files may not depend diff --git a/src/core/cuda_interface.hpp b/src/core/cuda_interface.hpp index 6bf505c1b8b..1f8142dcf54 100644 --- a/src/core/cuda_interface.hpp +++ b/src/core/cuda_interface.hpp @@ -28,17 +28,6 @@ #include #include -#ifdef ENGINE -// velocities which need to be copied from the GPU to the CPU to calculate a -// torque -typedef struct { - - // center and source velocity of the md part - float v_cs[6]; - -} CUDA_v_cs; -#endif - // Parameters for swimmers #ifdef ENGINE struct CUDA_ParticleParametersSwimming { @@ -156,13 +145,6 @@ void cuda_mpi_send_forces(ParticleRange particles, void cuda_bcast_global_part_params(); void cuda_copy_to_device(void *host_data, void *device_data, size_t n); void cuda_copy_to_host(void *host_device, void *device_host, size_t n); - -#ifdef ENGINE -void copy_v_cs_from_GPU(ParticleRange particles); -void cuda_mpi_send_v_cs(ParticleRange particles, - std::vector host_v_cs); -#endif - #endif /* ifdef CUDA */ #endif /* ifdef CUDA_INTERFACE_HPP */ diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index bee367c18cc..7c7c9965428 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -94,8 +94,8 @@ inline ParticleForce external_force(Particle const &p) { #ifdef ENGINE // apply a swimming force in the direction of // the particle's orientation axis - if (p.swim.swimming) { - f.f += p.swim.f_swim * p.r.calc_director(); + if (p.p.swim.swimming) { + f.f += p.p.swim.f_swim * p.r.calc_director(); } #endif diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp index b7fbfa11e4c..6710b76e020 100644 --- a/src/core/ghosts.cpp +++ b/src/core/ghosts.cpp @@ -183,10 +183,6 @@ static int calc_transmit_size(GhostCommunication &ghost_comm, if (data_parts & GHOSTTRANS_FORCE) n_buffer_new += sizeof(ParticleForce); -#ifdef ENGINE - if (data_parts & GHOSTTRANS_SWIMMING) - n_buffer_new += sizeof(ParticleParametersSwimming); -#endif int count = 0; for (auto const &pl : ghost_comm.part_lists) count += pl->n; @@ -231,12 +227,6 @@ static void prepare_send_buffer(CommBuf &send_buffer, if (data_parts & GHOSTTRANS_FORCE) { archiver << part.f; } - -#ifdef ENGINE - if (data_parts & GHOSTTRANS_SWIMMING) { - archiver << part.swim; - } -#endif } } } @@ -311,12 +301,6 @@ static void put_recv_buffer(CommBuf &recv_buffer, if (data_parts & GHOSTTRANS_FORCE) { archiver >> part.f; } - -#ifdef ENGINE - if (data_parts & GHOSTTRANS_SWIMMING) { - archiver >> part.swim; - } -#endif } } } @@ -368,11 +352,6 @@ static void cell_cell_transfer(GhostCommunication &ghost_comm, } if (data_parts & GHOSTTRANS_FORCE) part2.f += part1.f; - -#ifdef ENGINE - if (data_parts & GHOSTTRANS_SWIMMING) - part2.swim = part1.swim; -#endif } } } diff --git a/src/core/ghosts.hpp b/src/core/ghosts.hpp index c4f5f8f1137..dc2c4f58d6c 100644 --- a/src/core/ghosts.hpp +++ b/src/core/ghosts.hpp @@ -134,9 +134,7 @@ enum : unsigned { /// transfer \ref ParticleForce GHOSTTRANS_FORCE = 16u, /// resize the receiver particle arrays to the size of the senders - GHOSTTRANS_PARTNUM = 64u, - /// transfer \ref ParticleParametersSwimming - GHOSTTRANS_SWIMMING = 128u + GHOSTTRANS_PARTNUM = 64u }; /** \name Data Types */ diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index 9ada55db348..853a4a14be0 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.cpp @@ -146,8 +146,8 @@ Utils::Vector3d lb_viscous_coupling(Particle const &p, Utils::Vector3d v_drift = interpolated_u; #ifdef ENGINE - if (p.swim.swimming) { - v_drift += p.swim.v_swim * p.r.calc_director(); + if (p.p.swim.swimming) { + v_drift += p.p.swim.v_swim * p.r.calc_director(); } #endif @@ -222,9 +222,9 @@ bool in_local_halo(Vector3d const &pos) { #ifdef ENGINE void add_swimmer_force(Particle &p) { - if (p.swim.swimming) { + if (p.p.swim.swimming) { // calculate source position - const double direction = double(p.swim.push_pull) * p.swim.dipole_length; + const double direction = double(p.p.swim.push_pull) * p.p.swim.dipole_length; auto const director = p.r.calc_director(); auto const source_position = p.r.p + direction * director; @@ -232,7 +232,7 @@ void add_swimmer_force(Particle &p) { return; } - add_md_force(source_position, p.swim.f_swim * director); + add_md_force(source_position, p.p.swim.f_swim * director); } } #endif @@ -265,11 +265,6 @@ void lb_lbcoupling_calc_particle_lattice_ia( throw std::runtime_error("The non-linear interpolation scheme is not " "implemented for the CPU LB."); case (InterpolationOrder::linear): { -#ifdef ENGINE - ghost_communicator(&cell_structure.exchange_ghosts_comm, - GHOSTTRANS_SWIMMING); -#endif - using rng_type = r123::Philox4x64; using ctr_type = rng_type::ctr_type; using key_type = rng_type::key_type; diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index 79d67588d2d..fe8b4aa7e27 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -120,6 +120,9 @@ using UpdatePropertyMessage = boost::variant #ifdef LB_ELECTROHYDRODYNAMICS , UpdateProperty #endif +#ifdef ENGINE + , UpdateProperty +#endif #ifdef DIPOLES , UpdateProperty #endif @@ -224,21 +227,6 @@ using UpdateBondMessage = boost::variant , AddBond >; -#ifdef ENGINE -struct UpdateSwim { - ParticleParametersSwimming swim; - - void operator()(Particle &p) const { - p.swim = swim; - } - - template - void serialize(Archive &ar, long int) { - ar & swim; - } -}; -#endif - #ifdef ROTATION struct UpdateOrientation { Utils::Vector3d axis; @@ -273,9 +261,6 @@ using UpdateMessage = boost::variant , UpdateMomentumMessage , UpdateForceMessage , UpdateBondMessage -#ifdef ENGINE - , UpdateSwim -#endif #ifdef ROTATION , UpdateOrientation #endif @@ -791,7 +776,7 @@ void set_particle_v(int part, double *v) { #ifdef ENGINE void set_particle_swimming(int part, ParticleParametersSwimming swim) { - mpi_send_update_message(part, UpdateSwim{swim}); + mpi_update_particle_property(part, swim); } #endif @@ -1512,7 +1497,7 @@ void pointer_to_temperature(Particle const *p, double const *&res) { #ifdef ENGINE void pointer_to_swimming(Particle const *p, ParticleParametersSwimming const *&swim) { - swim = &(p->swim); + swim = &(p->p.swim); } #endif diff --git a/src/core/serialization/ParticleParametersSwimming.hpp b/src/core/serialization/ParticleParametersSwimming.hpp deleted file mode 100644 index 7a20f27344c..00000000000 --- a/src/core/serialization/ParticleParametersSwimming.hpp +++ /dev/null @@ -1,51 +0,0 @@ -/* - * Copyright (C) 2010-2019 The ESPResSo project - * - * This file is part of ESPResSo. - * - * ESPResSo is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * ESPResSo is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ -#ifndef CORE_UTILS_SERIALIZATION_PARTICLE_SWIM_HPP -#define CORE_UTILS_SERIALIZATION_PARTICLE_SWIM_HPP - -#include "Particle.hpp" - -#include - -namespace boost { -namespace serialization { -/* Pod serialize for ParticleParametersSwimming */ -template -void load(Archive &ar, ParticleParametersSwimming &swim, - const unsigned int /* file_version */) { - ar >> make_array(reinterpret_cast(&swim), - sizeof(ParticleParametersSwimming)); -} - -template -void save(Archive &ar, ParticleParametersSwimming const &swim, - const unsigned int /* file_version */) { - ar << make_array(reinterpret_cast(&swim), - sizeof(ParticleParametersSwimming)); -} - -template -void serialize(Archive &ar, ParticleParametersSwimming &swim, - const unsigned int file_version) { - split_free(ar, swim, file_version); -} -} // namespace serialization -} // namespace boost - -#endif diff --git a/src/core/thermostat.hpp b/src/core/thermostat.hpp index 82f42a0e202..39fdf6797d7 100644 --- a/src/core/thermostat.hpp +++ b/src/core/thermostat.hpp @@ -224,8 +224,8 @@ inline Utils::Vector3d friction_thermo_langevin(Particle const &p) { // Get velocity effective in the thermostatting #ifdef ENGINE - auto const velocity = (p.swim.v_swim != 0) - ? p.m.v - p.swim.v_swim * p.r.calc_director() + auto const velocity = (p.p.swim.v_swim != 0) + ? p.m.v - p.p.swim.v_swim * p.r.calc_director() : p.m.v; #else auto const &velocity = p.m.v; From c7e80240a5510086b9fd848e4c4cb44c5732b981 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 24 Oct 2019 23:08:50 +0200 Subject: [PATCH 3/6] doc: Removed mentions of shear torque from docs and tutorial --- doc/sphinx/particles.rst | 15 ++------------- .../06-active_matter/06-active_matter.tex | 8 +------- 2 files changed, 3 insertions(+), 20 deletions(-) diff --git a/doc/sphinx/particles.rst b/doc/sphinx/particles.rst index 8af73159866..1a0b7052134 100644 --- a/doc/sphinx/particles.rst +++ b/doc/sphinx/particles.rst @@ -535,7 +535,7 @@ Lattice Boltzmann (LB) swimmers system = espressomd.System() system.part.add(id=1, pos=[2, 0, 0], rotation=[1, 1, 1], swimming={ - 'f_swim': 0.01, 'mode': 'pusher', 'dipole_length': 2.0, 'rotational_friction': 20}) + 'f_swim': 0.01, 'mode': 'pusher', 'dipole_length': 2.0}) For an explanation of the parameters ``v_swim`` and ``f_swim`` see the previous item. In lattice Boltzmann self-propulsion is less trivial than for regular MD, @@ -549,11 +549,7 @@ nature of the particle's flow field by using one of the modes: ``pusher`` or ``puller``. You will also need to specify a ``dipole_length`` which determines the distance of the source of propulsion from the particle's center. Note that you should not put this distance to zero; |es| (currently) does not support -mathematical dipole flow fields. The key ``rotational_friction`` can be used to -set the friction that causes the orientation of the particle to change in shear -flow. The torque on the particle is determined by taking the cross product of -the difference between the fluid velocity at the center of the particle and at -the source point and the vector connecting the center and source. +mathematical dipole flow fields. You may ask: "Why are there two methods ``v_swim`` and ``f_swim`` for the self-propulsion using the lattice Boltzmann algorithm?" The answer is @@ -566,10 +562,3 @@ reaches a constant speed (given by ``v_swim``) this monopolar moment is gone and the flow field is zero! In contrast, ``f_swim`` always, i.e., while accelerating *and* while swimming at constant force possesses a dipolar flow field. - -.. warning:: - - Please note that even though swimming is interoperable with the - CPU version of LB it is only supported on *one* MPI - rank, i.e. ``n_nodes`` = 1. - diff --git a/doc/tutorials/06-active_matter/06-active_matter.tex b/doc/tutorials/06-active_matter/06-active_matter.tex index 951719c1b8c..2ac52096ded 100644 --- a/doc/tutorials/06-active_matter/06-active_matter.tex +++ b/doc/tutorials/06-active_matter/06-active_matter.tex @@ -230,13 +230,7 @@ \subsection{\label{sub:lattice}Self-Propulsion with Hydrodynamics} flow,~\textit{e.g.}, a swimmer in a Posseuille flow. If the thermostat is switched on, the rotational degrees of freedom will also be thermalized, but there is still no contribution of rotation due to `external' flow fields. A -rather unsatisfying solution to this problem is offered by the -\codees{rotational_friction} option, which allows one to include the effect of -fluid flow on the rotation. The algorithm computes the difference in fluid flow -velocity between the center of the particle and the coupling point, and takes -the cross product of that with the dipole length. This quantity is then -converted into a torque using the \codees{rotational_friction}. It is -recommended to use an alternative means of obtaining rotations in your LB +It is recommended to use an alternative means of obtaining rotations in your LB swimming simulations. For example, by constructing a raspberry particle~\cite{lobaskin04,chatterji05,fischer15,degraaf15,degraaf16b}. From 73f76a1b7f25286519520d993e45e933a106f027 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 24 Oct 2019 23:28:55 +0200 Subject: [PATCH 4/6] Formatting --- src/core/cuda_interface.cpp | 3 ++- src/core/grid_based_algorithms/lb_particle_coupling.cpp | 3 ++- src/core/particle_data.cpp | 3 ++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/core/cuda_interface.cpp b/src/core/cuda_interface.cpp index e07f68d85cc..2d8171a68e9 100644 --- a/src/core/cuda_interface.cpp +++ b/src/core/cuda_interface.cpp @@ -82,7 +82,8 @@ static void pack_particles(ParticleRange particles, buffer[i].swim.director = buffer[i].director; buffer[i].swim.push_pull = part.p.swim.push_pull; - buffer[i].swim.dipole_length = static_cast(part.p.swim.dipole_length); + buffer[i].swim.dipole_length = + static_cast(part.p.swim.dipole_length); buffer[i].swim.swimming = part.p.swim.swimming; #endif i++; diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index 853a4a14be0..c97e4e39eb0 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.cpp @@ -224,7 +224,8 @@ bool in_local_halo(Vector3d const &pos) { void add_swimmer_force(Particle &p) { if (p.p.swim.swimming) { // calculate source position - const double direction = double(p.p.swim.push_pull) * p.p.swim.dipole_length; + const double direction = + double(p.p.swim.push_pull) * p.p.swim.dipole_length; auto const director = p.r.calc_director(); auto const source_position = p.r.p + direction * director; diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index fe8b4aa7e27..d91a1646e19 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -776,7 +776,8 @@ void set_particle_v(int part, double *v) { #ifdef ENGINE void set_particle_swimming(int part, ParticleParametersSwimming swim) { - mpi_update_particle_property(part, swim); + mpi_update_particle_property(part, swim); } #endif From 621aaab287ca84c3ecbc500ddf34caecf6773f68 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 25 Oct 2019 00:19:36 +0200 Subject: [PATCH 5/6] core: Removed unused headers vom rotation --- src/core/rotation.cpp | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/src/core/rotation.cpp b/src/core/rotation.cpp index b480f2fd0a8..73fad9d5158 100644 --- a/src/core/rotation.cpp +++ b/src/core/rotation.cpp @@ -35,23 +35,12 @@ #include "rotation.hpp" #ifdef ROTATION -#include "Particle.hpp" -#include "cells.hpp" -#include "communication.hpp" -#include "cuda_interface.hpp" -#include "errorhandling.hpp" -#include "forces.hpp" -#include "global.hpp" -#include "grid_based_algorithms/lb_interface.hpp" #include "integrate.hpp" -#include "particle_data.hpp" #include -#include -#include +#include #include -#include /** Calculate the derivatives of the quaternion and angular acceleration * for a given particle @@ -63,10 +52,7 @@ * @param[out] Wd Angular acceleration of the particle */ static void define_Qdd(Particle const &p, double Qd[4], double Qdd[4], - double S[3], double Wd[3]); - -void define_Qdd(Particle const &p, double Qd[4], double Qdd[4], double S[3], - double Wd[3]) { + double S[3], double Wd[3]) { /* calculate the first derivative of the quaternion */ /* Taken from "An improved algorithm for molecular dynamics simulation of * rigid molecules", Sonnenschein, Roland (1985), Eq. 4.*/ From dc9202d7b6acc50aac36ec01e75f5836958c3097 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 25 Oct 2019 21:39:56 +0200 Subject: [PATCH 6/6] testsuite: Fixed engine_lb test return code core/testsuite: engine_lb test fix from christophlohrmann --- src/core/integrate.cpp | 6 ++- testsuite/python/engine_lb.py | 75 ++++------------------------------- 2 files changed, 11 insertions(+), 70 deletions(-) diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 517bb524644..cd036e68074 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -200,8 +200,10 @@ void integrate_vv(int n_steps, int reuse_forces) { thermo_cool_down(); ESPRESSO_PROFILER_MARK_END("Initial Force Calculation"); - if (n_part > 0) - lb_lbcoupling_activate(); + } + + if (n_part > 0) { + lb_lbcoupling_activate(); } if (check_runtime_errors(comm_cart)) diff --git a/testsuite/python/engine_lb.py b/testsuite/python/engine_lb.py index f7a8ac2dfae..fbdb83a26ff 100644 --- a/testsuite/python/engine_lb.py +++ b/testsuite/python/engine_lb.py @@ -90,6 +90,10 @@ def test_momentum_conservation(self): self.system.integrator.run(20, reuse_forces=True) tot_mom = self.system.analysis.linear_momentum(include_particles=True, include_lbfluid=True) + # compensate half-step offset between force calculation and LB-update + for part in self.system.part: + tot_mom += part.f * self.system.time_step / 2. + np.testing.assert_allclose(tot_mom, 3 * [0.], atol=self.tol) def test_particle_forces(self): @@ -114,15 +118,14 @@ def test_particle_forces(self): f_swim * director + self.gamma * v_swim * director self.system.integrator.run(1, reuse_forces=True) - np.testing.assert_allclose(swimmer.f, force, atol=self.tol) + np.testing.assert_allclose( + np.copy(swimmer.f), force, atol=self.tol) def check_fluid_force(self, swimmer): pass # forces on particles are checked # total force on the fluid matches (momentum conservation) # TODO: only thing left to check is the location of the fluid force. - # There is no counter torque when using the rotational_friction-feature - # so there is nothing to be tested @utx.skipIfMissingFeatures(["ENGINE", "ROTATION", "MASS"]) @@ -134,22 +137,6 @@ def setUp(self): self.system.actors.add(self.lbf) self.system.thermostat.set_lb(LB_fluid=self.lbf, gamma=self.gamma) - def test_rotfric_exception(self): - """rotational_friction feature is disabled on CPU for more than one core - """ - if self.system.cell_system.get_state()["n_nodes"] > 1: - # swimming without rot_fric is fine - self.system.part.add(pos=[0, 0, 0], rotation=3 * [True], - swimming={"f_swim": 0.1, "mode": "pusher"}) - self.system.integrator.run(3) - # with rot_fric it is not - with self.assertRaises(Exception): - self.system.part.add(pos=[0, 0, 0], rotation=3 * [True], - swimming={"f_swim": 0.1, - "mode": "pusher", - "rotational_friction": 0.3}) - self.system.integrator.run(3) - @utx.skipIfMissingGPU() @utx.skipIfMissingFeatures(["ENGINE", "ROTATION", "MASS"]) @@ -161,54 +148,6 @@ def setUp(self): self.system.actors.add(self.lbf) self.system.thermostat.set_lb(LB_fluid=self.lbf, gamma=self.gamma) - def test_particle_torques(self): - """setup shear flow and check if resulting torques match - the formulae in the core - """ - - bottom = shapes.Wall(normal=[0, 0, 1], - dist=self.LB_params['agrid']) - top = shapes.Wall(normal=[0, 0, -1], - dist=-self.system.box_l[2] + self.LB_params['agrid']) - self.system.lbboundaries.add(lbboundaries.LBBoundary(shape=bottom)) - self.system.lbboundaries.add( - lbboundaries.LBBoundary(shape=top, velocity=[1e-3, 1e-3, 0])) - self.system.integrator.run(100) - - # fix the particles so inaccuracies from position updates - # before torque calculation don't matter - self.add_all_types_of_swimmers(fix=True, rotation=True, - put_in_corners=False) - - self.system.integrator.run(20) - for swimmer in self.system.part: - - director = swimmer.director - dip_len = swimmer.swimming["dipole_length"] - mode_fac = 1. if swimmer.swimming["mode"] == "puller" else -1. - source_pos = swimmer.pos + mode_fac * dip_len * director - - v_center = self.lbf.get_interpolated_velocity(swimmer.pos) - v_source = self.lbf.get_interpolated_velocity(source_pos) - diff = v_center - v_source - cross = np.cross(diff, director) - # half-step omega with isotropic rinertia - omega_part = swimmer.omega_lab + 0.5 * self.system.time_step * \ - swimmer.torque_lab / swimmer.rinertia[0] - omega_swim = cross / np.linalg.norm(cross) * \ - np.linalg.norm(diff) / dip_len - torque = swimmer.swimming["rotational_friction"] * \ - (omega_swim - omega_part) - - self.system.integrator.run(1, reuse_forces=True) - np.testing.assert_allclose( - swimmer.torque_lab, - torque, - atol=self.tol) - if __name__ == "__main__": - suite = ut.TestSuite() - suite.addTests(ut.TestLoader().loadTestsFromTestCase(SwimmerTestGPU)) - suite.addTests(ut.TestLoader().loadTestsFromTestCase(SwimmerTestCPU)) - result = ut.TextTestRunner(verbosity=4).run(suite) + ut.main()