From 452a6e62de3d00414e718c6a49b17c16b691a844 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 2 Jun 2023 20:30:10 +0200 Subject: [PATCH 1/2] core: Refactor LB particle coupling --- src/core/forces.cpp | 3 +- .../lb_particle_coupling.cpp | 121 ++++++------------ .../lb_particle_coupling.hpp | 101 +++++++++++---- .../unit_tests/lb_particle_coupling_test.cpp | 59 ++++----- .../VirtualSitesInertialessTracers.cpp | 35 ++--- 5 files changed, 153 insertions(+), 166 deletions(-) diff --git a/src/core/forces.cpp b/src/core/forces.cpp index b2db558000e..5c8e3c83a27 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -225,8 +225,7 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) { immersed_boundaries.volume_conservation(cell_structure); if (lattice_switch != ActiveLB::NONE) { - lb_lbcoupling_calc_particle_lattice_ia(thermo_virtual, particles, - ghost_particles, time_step); + LB::couple_particles(thermo_virtual, particles, ghost_particles, time_step); } #ifdef CUDA diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index a9de5d30392..511785eb3c1 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.cpp @@ -32,17 +32,16 @@ #include #include #include -#include #include -#include -#include -#include +#include +#include +#include #include -#include +#include -LB_Particle_Coupling lb_particle_coupling; +LB::ParticleCouplingConfig lb_particle_coupling; void mpi_bcast_lb_particle_coupling_local() { boost::mpi::broadcast(comm_cart, lb_particle_coupling, 0); @@ -190,36 +189,6 @@ std::vector positions_in_halo(Utils::Vector3d pos, return res; } -/** - * @brief Return if locally there exists a physical particle - * for a given (ghost) particle - */ -bool is_ghost_for_local_particle(const Particle &p) { - return !cell_structure.get_local_particle(p.id())->is_ghost(); -} - -/** - * @brief Determine if a given particle should be coupled. - * In certain cases, there may be more than one ghost for the same particle. - * To make sure, that these are only coupled once, ghosts' ids are stored - * in an unordered_set. - */ -bool should_be_coupled(const Particle &p, - std::unordered_set &coupled_ghost_particles) { - // always couple physical particles - if (not p.is_ghost()) { - return true; - } - // for ghosts check that we don't have the physical particle and - // that a ghost for the same particle has not yet been coupled - if ((not is_ghost_for_local_particle(p)) and - (coupled_ghost_particles.find(p.id()) == coupled_ghost_particles.end())) { - coupled_ghost_particles.insert(p.id()); - return true; - } - return false; -} - #ifdef ENGINE void add_swimmer_force(Particle const &p, double time_step) { if (p.swimming().swimming) { @@ -239,34 +208,32 @@ void add_swimmer_force(Particle const &p, double time_step) { } #endif -Utils::Vector3d lb_particle_coupling_noise(bool enabled, int part_id, - const OptionalCounter &rng_counter) { - if (enabled) { - if (rng_counter) { - return Random::noise_uniform(rng_counter->value(), 0, - part_id); - } +namespace LB { + +Utils::Vector3d ParticleCoupling::get_noise_term(int pid) const { + if (not m_thermalized) { + return Utils::Vector3d{}; + } + auto const &rng_counter = lb_particle_coupling.rng_counter_coupling; + if (not rng_counter) { throw std::runtime_error( "Access to uninitialized LB particle coupling RNG counter"); } - return {}; + auto const counter = rng_counter->value(); + return m_noise * Random::noise_uniform(counter, 0, pid); } -void couple_particle(Particle &p, bool couple_virtual, double noise_amplitude, - const OptionalCounter &rng_counter, double time_step) { - - if (p.is_virtual() and not couple_virtual) +void ParticleCoupling::kernel(Particle &p) { + if (p.is_virtual() and not m_couple_virtual) return; // Calculate coupling force Utils::Vector3d coupling_force = {}; for (auto pos : positions_in_halo(p.pos(), box_geo)) { if (in_local_halo(pos)) { - auto const drag_force = - lb_drag_force(p, pos, lb_particle_coupling_drift_vel_offset(p)); - auto const random_force = - noise_amplitude * lb_particle_coupling_noise(noise_amplitude > 0.0, - p.id(), rng_counter); + auto const vel_offset = lb_particle_coupling_drift_vel_offset(p); + auto const drag_force = lb_drag_force(p, pos, vel_offset); + auto const random_force = get_noise_term(p.id()); coupling_force = drag_force + random_force; break; } @@ -280,52 +247,38 @@ void couple_particle(Particle &p, bool couple_virtual, double noise_amplitude, * is responsible to adding its force */ p.force() += coupling_force; } - add_md_force(pos, coupling_force, time_step); + add_md_force(pos, coupling_force, m_time_step); } #ifdef ENGINE - add_swimmer_force(p, time_step); + add_swimmer_force(p, m_time_step); #endif } -void lb_lbcoupling_calc_particle_lattice_ia(bool couple_virtual, - const ParticleRange &particles, - const ParticleRange &more_particles, - double time_step) { +bool CouplingBookkeeping::is_ghost_for_local_particle(Particle const &p) const { + return not ::cell_structure.get_local_particle(p.id())->is_ghost(); +} + +void couple_particles(bool couple_virtual, ParticleRange const &real_particles, + ParticleRange const &ghost_particles, double time_step) { ESPRESSO_PROFILER_CXX_MARK_FUNCTION; if (lattice_switch == ActiveLB::WALBERLA_LB) { if (lb_particle_coupling.couple_to_md) { - auto const kT = LB::get_kT() * Utils::sqr(LB::get_lattice_speed()); - /* Eq. (16) @cite ahlrichs99a. - * The factor 12 comes from the fact that we use random numbers - * from -0.5 to 0.5 (equally distributed) which have variance 1/12. - * time_step comes from the discretization. - */ - auto const noise_amplitude = - (kT > 0.) - ? std::sqrt(12. * 2. * lb_lbcoupling_get_gamma() * kT / time_step) - : 0.0; - - std::unordered_set coupled_ghost_particles; - - /* Couple particles ranges */ - for (auto &p : particles) { - if (should_be_coupled(p, coupled_ghost_particles)) { - couple_particle(p, couple_virtual, noise_amplitude, - lb_particle_coupling.rng_counter_coupling, time_step); - } - } - - for (auto &p : more_particles) { - if (should_be_coupled(p, coupled_ghost_particles)) { - couple_particle(p, couple_virtual, noise_amplitude, - lb_particle_coupling.rng_counter_coupling, time_step); + ParticleCoupling coupling{couple_virtual, time_step}; + CouplingBookkeeping bookkeeping{}; + for (auto const &particle_range : {real_particles, ghost_particles}) { + for (auto &p : particle_range) { + if (bookkeeping.should_be_coupled(p)) { + coupling.kernel(p); + } } } } } } +} // namespace LB + void lb_lbcoupling_propagate() { if (lattice_switch != ActiveLB::NONE) { if (LB::get_kT() > 0.0) { diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.hpp b/src/core/grid_based_algorithms/lb_particle_coupling.hpp index c93a98ca008..35125aabb23 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.hpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.hpp @@ -25,26 +25,20 @@ #include #include +#include #include #include #include +#include +#include #include #include #include using OptionalCounter = boost::optional>; -/** Calculate particle lattice interactions. - * So far, only viscous coupling with Stokesian friction is implemented. - * Include all particle-lattice forces in this function. - * The function is called from \ref force_calc. - */ -void lb_lbcoupling_calc_particle_lattice_ia(bool couple_virtual, - const ParticleRange &particles, - const ParticleRange &more_particles, - double time_step); void lb_lbcoupling_propagate(); uint64_t lb_lbcoupling_get_rng_state(); void lb_lbcoupling_set_rng_state(uint64_t counter); @@ -75,14 +69,6 @@ void lb_lbcoupling_deactivate(); */ bool in_local_halo(Utils::Vector3d const &pos); -/** @brief Determine if a given particle should be coupled. - * In certain cases, there may be more than one ghost for the same particle. - * To make sure, that these are only coupled once, ghosts' ids are stored - * in an unordered_set. - */ -bool should_be_coupled(const Particle &p, - std::unordered_set &coupled_ghost_particles); - /** * @brief Add a force to the lattice force density. * @param pos Position of the force @@ -92,17 +78,10 @@ bool should_be_coupled(const Particle &p, void add_md_force(Utils::Vector3d const &pos, Utils::Vector3d const &force, double time_step); -Utils::Vector3d lb_particle_coupling_noise(bool enabled, int part_id, - const OptionalCounter &rng_counter); - // internal function exposed for unit testing std::vector positions_in_halo(Utils::Vector3d pos, const BoxGeometry &box); -// internal function exposed for unit testing -void couple_particle(Particle &p, bool couple_virtual, double noise_amplitude, - const OptionalCounter &rng_counter, double time_step); - // internal function exposed for unit testing void add_swimmer_force(Particle const &p, double time_step); @@ -129,7 +108,8 @@ Utils::Vector3d lb_drag_force(Particle const &p, Utils::Vector3d const &shifted_pos, Utils::Vector3d const &vel_offset); -struct LB_Particle_Coupling { +namespace LB { +struct ParticleCouplingConfig { OptionalCounter rng_counter_coupling = {}; /** @brief Friction coefficient for the particle coupling. */ double gamma = 0.0; @@ -144,8 +124,77 @@ struct LB_Particle_Coupling { ar &couple_to_md; } }; +} // namespace LB // internal global exposed for unit testing -extern LB_Particle_Coupling lb_particle_coupling; +extern LB::ParticleCouplingConfig lb_particle_coupling; + +namespace LB { + +/** @brief Calculate particle-lattice interactions. */ +void couple_particles(bool couple_virtual, ParticleRange const &real_particles, + ParticleRange const &ghost_particles, double time_step); + +class ParticleCoupling { + bool m_couple_virtual; + bool m_thermalized; + double m_time_step; + double m_noise; + +public: + auto get_noise(double kT, double gamma) const { + /* Eq. (16) @cite ahlrichs99a. + * The factor 12 comes from the fact that we use random numbers + * from -0.5 to 0.5 (equally distributed) which have variance 1/12. + * The time step comes from the discretization. + */ + auto constexpr variance_correction = 12.; + return std::sqrt(variance_correction * 2. * gamma * kT / m_time_step); + } + + ParticleCoupling(bool couple_virtual, double time_step, double kT) + : m_couple_virtual{couple_virtual}, m_thermalized{kT != 0.}, + m_time_step{time_step} { + assert(kT >= 0.); + m_noise = get_noise(kT, lb_particle_coupling.gamma); + } + + ParticleCoupling(bool couple_virtual, double time_step) + : ParticleCoupling(couple_virtual, time_step, + LB::get_kT() * Utils::sqr(LB::get_lattice_speed())) {} + + Utils::Vector3d get_noise_term(int pid) const; + void kernel(Particle &p); +}; + +/** + * @brief Keep track of ghost particles that have already been coupled once. + * In certain cases, there may be more than one ghost for the same particle. + * To make sure these are only coupled once, ghosts' ids are recorded. + */ +class CouplingBookkeeping { + std::unordered_set m_coupled_ghosts; + + /** @brief Check if there is locally a real particle for the given ghost. */ + bool is_ghost_for_local_particle(Particle const &p) const; + +public: + /** @brief Determine if a given particle should be coupled. */ + bool should_be_coupled(Particle const &p) { + // real particles: always couple + if (not p.is_ghost()) { + return true; + } + // ghosts: check we don't have the corresponding real particle on the same + // node, and that a ghost for the same particle hasn't been coupled already + if (m_coupled_ghosts.count(p.id()) != 0 or is_ghost_for_local_particle(p)) { + return false; + } + m_coupled_ghosts.insert(p.id()); + return true; + } +}; + +} // namespace LB #endif diff --git a/src/core/unit_tests/lb_particle_coupling_test.cpp b/src/core/unit_tests/lb_particle_coupling_test.cpp index 9a0a50dbdc3..bf20400e98e 100644 --- a/src/core/unit_tests/lb_particle_coupling_test.cpp +++ b/src/core/unit_tests/lb_particle_coupling_test.cpp @@ -58,7 +58,6 @@ namespace utf = boost::unit_test; #include #include -#include #include #include #include @@ -178,15 +177,16 @@ BOOST_AUTO_TEST_CASE(rng_initial_state) { BOOST_AUTO_TEST_CASE(rng) { lb_lbcoupling_set_rng_state(17); + lb_particle_coupling.gamma = 0.2; + + LB::ParticleCoupling coupling{true, params.time_step, 1.}; BOOST_REQUIRE(lb_particle_coupling.rng_counter_coupling); BOOST_CHECK_EQUAL(lb_lbcoupling_get_rng_state(), 17); BOOST_CHECK(not lb_lbcoupling_is_seed_required()); - auto const step1_random1 = lb_particle_coupling_noise( - true, 1, lb_particle_coupling.rng_counter_coupling); - auto const step1_random2 = lb_particle_coupling_noise( - true, 4, lb_particle_coupling.rng_counter_coupling); - auto const step1_random2_try2 = lb_particle_coupling_noise( - true, 4, lb_particle_coupling.rng_counter_coupling); + auto const step1_random1 = coupling.get_noise_term(1); + auto const step1_random2 = coupling.get_noise_term(4); + auto const step1_random2_try2 = coupling.get_noise_term(4); + BOOST_REQUIRE(step1_random1.norm() > 1e-10); BOOST_CHECK(step1_random1 != step1_random2); BOOST_CHECK(step1_random2 == step1_random2_try2); @@ -196,15 +196,13 @@ BOOST_AUTO_TEST_CASE(rng) { BOOST_REQUIRE(lb_particle_coupling.rng_counter_coupling); BOOST_CHECK_EQUAL(lb_lbcoupling_get_rng_state(), 18); - auto const step2_random1 = lb_particle_coupling_noise( - true, 1, lb_particle_coupling.rng_counter_coupling); - auto const step2_random2 = lb_particle_coupling_noise( - true, 4, lb_particle_coupling.rng_counter_coupling); + auto const step2_random1 = coupling.get_noise_term(1); + auto const step2_random2 = coupling.get_noise_term(4); BOOST_CHECK(step1_random1 != step2_random1); BOOST_CHECK(step1_random1 != step2_random2); - auto const step3_norandom = lb_particle_coupling_noise( - false, 4, lb_particle_coupling.rng_counter_coupling); + LB::ParticleCoupling coupling_unthermalized{true, params.time_step, 0.}; + auto const step3_norandom = coupling_unthermalized.get_noise_term(4); BOOST_CHECK((step3_norandom == Utils::Vector3d{0., 0., 0.})); } @@ -317,12 +315,10 @@ BOOST_DATA_TEST_CASE(particle_coupling, bdata::make(kTs), kT) { auto const first_lb_node = lb_walberla()->get_lattice().get_local_domain().first; auto const gamma = 0.2; - auto const noise = - (kT > 0.) ? std::sqrt(24. * gamma * kT / params.time_step) : 0.0; - auto &rng = lb_particle_coupling.rng_counter_coupling; + lb_lbcoupling_set_gamma(gamma); Particle p{}; - Utils::Vector3d expected = noise * Random::noise_uniform( - rng->value(), 0, p.id()); + LB::ParticleCoupling coupling{false, params.time_step}; + auto expected = coupling.get_noise_term(p.id()); #ifdef ENGINE p.swimming().swimming = true; p.swimming().v_swim = 2.; @@ -334,12 +330,11 @@ BOOST_DATA_TEST_CASE(particle_coupling, bdata::make(kTs), kT) { expected += gamma * p.mu_E(); #endif p.pos() = first_lb_node + Utils::Vector3d::broadcast(0.5); - lb_lbcoupling_set_gamma(gamma); // coupling { if (in_local_halo(p.pos())) { - couple_particle(p, false, noise, rng, params.time_step); + coupling.kernel(p); BOOST_CHECK_SMALL((p.force() - expected).norm(), eps); auto const interpolated = LB::get_force_to_be_applied(p.pos()); @@ -365,10 +360,6 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, auto const first_lb_node = lb_walberla()->get_lattice().get_local_domain().first; auto const gamma = 0.2; - auto const noise = std::sqrt(24. * gamma * kT / params.time_step * - Utils::sqr(params.agrid / params.tau)); - auto &rng = lb_particle_coupling.rng_counter_coupling; - auto const pid = 0; auto const skin = params.skin; auto const &box_l = params.box_dimensions; @@ -387,8 +378,8 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, set_particle_property(pid, &Particle::mu_E, Utils::Vector3d{-2., 1.5, 1.}); #endif - auto expected = - noise * Random::noise_uniform(rng->value(), 0, pid); + LB::ParticleCoupling coupling{thermo_virtual, params.time_step}; + auto expected = coupling.get_noise_term(pid); auto const p_opt = copy_particle_to_head_node(comm, pid); #if defined(ENGINE) or defined(LB_ELECTROHYDRODYNAMICS) if (rank == 0) { @@ -462,8 +453,8 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, lb_lbcoupling_broadcast(); auto const particles = ::cell_structure.local_particles(); auto const ghost_particles = ::cell_structure.ghost_particles(); - lb_lbcoupling_calc_particle_lattice_ia(thermo_virtual, particles, - ghost_particles, params.time_step); + LB::couple_particles(thermo_virtual, particles, ghost_particles, + params.time_step); auto const p_opt = copy_particle_to_head_node(comm, pid); if (rank == 0) { auto const &p = *p_opt; @@ -487,8 +478,8 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, } } // couple particle to LB - lb_lbcoupling_calc_particle_lattice_ia(thermo_virtual, particles, - ghost_particles, params.time_step); + LB::couple_particles(thermo_virtual, particles, ghost_particles, + params.time_step); { auto const p_opt = copy_particle_to_head_node(comm, pid); if (rank == 0) { @@ -519,6 +510,12 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, for (auto const &error_message : error_messages) { BOOST_CHECK_EQUAL(error_message.what(), error_message_ref); } + lb_particle_coupling.rng_counter_coupling = {}; + if (kT == 0.) { + BOOST_CHECK_EQUAL(coupling.get_noise_term(0).norm(), 0.); + } else { + BOOST_CHECK_THROW(coupling.get_noise_term(0), std::runtime_error); + } } } @@ -565,8 +562,6 @@ BOOST_AUTO_TEST_CASE(exceptions) { // coupling, interpolation, boundaries BOOST_CHECK_THROW(lb_lbcoupling_get_rng_state(), std::runtime_error); BOOST_CHECK_THROW(lb_lbcoupling_set_rng_state(0ul), std::runtime_error); - BOOST_CHECK_THROW(lb_particle_coupling_noise(true, 0, OptionalCounter{}), - std::runtime_error); BOOST_CHECK_THROW(lb_lbinterpolation_get_interpolated_velocity({}), std::runtime_error); BOOST_CHECK_THROW(lb_lbinterpolation_add_force_density({}, {}), diff --git a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp b/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp index 6c06d9d6657..138b55d2ecc 100644 --- a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp +++ b/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp @@ -30,7 +30,7 @@ #include "grid_based_algorithms/lb_particle_coupling.hpp" #include "integrate.hpp" -#include +#include static bool lb_active_check() { if (lattice_switch == ActiveLB::NONE) { @@ -49,30 +49,21 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) { cells_update_ghosts(Cells::DATA_PART_FORCE); // Set to store ghost particles (ids) that have already been coupled - std::unordered_set coupled_ghost_particles; + LB::CouplingBookkeeping bookkeeping{}; // Apply particle forces to the LB fluid at particle positions // For physical particles, also set particle velocity = fluid velocity - for (auto &p : cell_structure.local_particles()) { - if (!p.is_virtual()) - continue; - if (!lb_active_check()) { - return; - } - if (should_be_coupled(p, coupled_ghost_particles)) { - for (auto pos : positions_in_halo(p.pos(), box_geo)) { - add_md_force(pos * to_lb_units, -p.force(), time_step); + for (auto const &particle_range : + {cell_structure.local_particles(), cell_structure.ghost_particles()}) { + for (auto const &p : particle_range) { + if (!p.is_virtual()) + continue; + if (!lb_active_check()) { + return; } - } - } - for (auto const &p : cell_structure.ghost_particles()) { - if (!p.is_virtual()) - continue; - if (!lb_active_check()) { - return; - } - if (should_be_coupled(p, coupled_ghost_particles)) { - for (auto pos : positions_in_halo(p.pos(), box_geo)) { - add_md_force(pos * to_lb_units, -p.force(), time_step); + if (bookkeeping.should_be_coupled(p)) { + for (auto pos : positions_in_halo(p.pos(), box_geo)) { + add_md_force(pos * to_lb_units, -p.force(), time_step); + } } } } From 3c2cb0be68651656f18f95751be1e7a15b62eda7 Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Mon, 12 Jun 2023 10:24:04 +0200 Subject: [PATCH 2/2] LB: per particle gamma --- .../lb_particle_coupling.cpp | 50 ++++- .../lb_particle_coupling.hpp | 22 +-- src/core/thermostat.hpp | 65 +++++-- src/core/thermostats/langevin_inline.hpp | 56 ++---- .../unit_tests/lb_particle_coupling_test.cpp | 29 +-- testsuite/python/lb_thermostat.py | 178 ++++++++++++++---- testsuite/python/thermostats_common.py | 18 +- 7 files changed, 286 insertions(+), 132 deletions(-) diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index 511785eb3c1..6bc67537230 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.cpp @@ -24,6 +24,7 @@ #include "errorhandling.hpp" #include "grid.hpp" #include "random.hpp" +#include "thermostat.hpp" #include "grid_based_algorithms/lb_interface.hpp" #include "grid_based_algorithms/lb_interpolation.hpp" @@ -35,6 +36,7 @@ #include +#include #include #include #include @@ -121,6 +123,21 @@ Utils::Vector3d lb_particle_coupling_drift_vel_offset(const Particle &p) { return vel_offset; } +static Thermostat::GammaType lb_handle_particle_anisotropy(Particle const &p) { + auto const lb_gamma = lb_lbcoupling_get_gamma(); +#ifdef THERMOSTAT_PER_PARTICLE + auto const &partcl_gamma = p.gamma(); +#ifdef PARTICLE_ANISOTROPY + auto const default_gamma = Thermostat::GammaType::broadcast(lb_gamma); +#else + auto const default_gamma = lb_gamma; +#endif // PARTICLE_ANISOTROPY + return handle_particle_gamma(partcl_gamma, default_gamma); +#else + return lb_gamma; +#endif // THERMOSTAT_PER_PARTICLE +} + Utils::Vector3d lb_drag_force(Particle const &p, Utils::Vector3d const &shifted_pos, Utils::Vector3d const &vel_offset) { @@ -131,8 +148,10 @@ Utils::Vector3d lb_drag_force(Particle const &p, LB::get_lattice_speed(); Utils::Vector3d v_drift = interpolated_u + vel_offset; + auto const gamma = lb_handle_particle_anisotropy(p); + /* calculate viscous force (eq. (9) @cite ahlrichs99a) */ - return -lb_lbcoupling_get_gamma() * (p.v() - v_drift); + return Utils::hadamard_product(gamma, v_drift - p.v()); } /** @@ -210,7 +229,7 @@ void add_swimmer_force(Particle const &p, double time_step) { namespace LB { -Utils::Vector3d ParticleCoupling::get_noise_term(int pid) const { +Utils::Vector3d ParticleCoupling::get_noise_term(Particle const &p) const { if (not m_thermalized) { return Utils::Vector3d{}; } @@ -219,8 +238,14 @@ Utils::Vector3d ParticleCoupling::get_noise_term(int pid) const { throw std::runtime_error( "Access to uninitialized LB particle coupling RNG counter"); } + using std::sqrt; + using Utils::sqrt; auto const counter = rng_counter->value(); - return m_noise * Random::noise_uniform(counter, 0, pid); + auto const gamma = lb_handle_particle_anisotropy(p); + return m_noise_pref_wo_gamma * + Utils::hadamard_product( + sqrt(gamma), + Random::noise_uniform(counter, 0, p.id())); } void ParticleCoupling::kernel(Particle &p) { @@ -233,7 +258,7 @@ void ParticleCoupling::kernel(Particle &p) { if (in_local_halo(pos)) { auto const vel_offset = lb_particle_coupling_drift_vel_offset(p); auto const drag_force = lb_drag_force(p, pos, vel_offset); - auto const random_force = get_noise_term(p.id()); + auto const random_force = get_noise_term(p); coupling_force = drag_force + random_force; break; } @@ -259,6 +284,20 @@ bool CouplingBookkeeping::is_ghost_for_local_particle(Particle const &p) const { return not ::cell_structure.get_local_particle(p.id())->is_ghost(); } +#if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY) +static void lb_coupling_sanity_checks(Particle const &p) { + /* + lb does (at the moment) not support rotational particle coupling. + Consequently, anisotropic particles are also not supported. + */ + auto const &p_gamma = p.gamma(); + if (p_gamma[0] != p_gamma[1] or p_gamma[1] != p_gamma[2]) { + runtimeErrorMsg() << "anisotropic particle (id " << p.id() + << ") coupled to LB."; + } +} +#endif + void couple_particles(bool couple_virtual, ParticleRange const &real_particles, ParticleRange const &ghost_particles, double time_step) { ESPRESSO_PROFILER_CXX_MARK_FUNCTION; @@ -269,6 +308,9 @@ void couple_particles(bool couple_virtual, ParticleRange const &real_particles, for (auto const &particle_range : {real_particles, ghost_particles}) { for (auto &p : particle_range) { if (bookkeeping.should_be_coupled(p)) { +#if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY) + lb_coupling_sanity_checks(p); +#endif coupling.kernel(p); } } diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.hpp b/src/core/grid_based_algorithms/lb_particle_coupling.hpp index 35125aabb23..f0b0a6d2821 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.hpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.hpp @@ -139,31 +139,27 @@ class ParticleCoupling { bool m_couple_virtual; bool m_thermalized; double m_time_step; - double m_noise; + double m_noise_pref_wo_gamma; public: - auto get_noise(double kT, double gamma) const { - /* Eq. (16) @cite ahlrichs99a. - * The factor 12 comes from the fact that we use random numbers - * from -0.5 to 0.5 (equally distributed) which have variance 1/12. - * The time step comes from the discretization. - */ - auto constexpr variance_correction = 12.; - return std::sqrt(variance_correction * 2. * gamma * kT / m_time_step); - } - ParticleCoupling(bool couple_virtual, double time_step, double kT) : m_couple_virtual{couple_virtual}, m_thermalized{kT != 0.}, m_time_step{time_step} { assert(kT >= 0.); - m_noise = get_noise(kT, lb_particle_coupling.gamma); + /* Eq. (16) @cite ahlrichs99a, without the gamma term. + * The factor 12 comes from the fact that we use random numbers + * from -0.5 to 0.5 (equally distributed) which have variance 1/12. + * The time step comes from the discretization. + */ + auto constexpr variance_inv = 12.; + m_noise_pref_wo_gamma = std::sqrt(variance_inv * 2. * kT / time_step); } ParticleCoupling(bool couple_virtual, double time_step) : ParticleCoupling(couple_virtual, time_step, LB::get_kT() * Utils::sqr(LB::get_lattice_speed())) {} - Utils::Vector3d get_noise_term(int pid) const; + Utils::Vector3d get_noise_term(Particle const &p) const; void kernel(Particle &p); }; diff --git a/src/core/thermostat.hpp b/src/core/thermostat.hpp index 45e0726b313..dfd0c507c36 100644 --- a/src/core/thermostat.hpp +++ b/src/core/thermostat.hpp @@ -24,6 +24,9 @@ * Implementation in \ref thermostat.cpp. */ +#include "Particle.hpp" +#include "rotation.hpp" + #include "config/config.hpp" #include @@ -55,16 +58,24 @@ using GammaType = double; } // namespace Thermostat namespace { -/** @name Integrators parameters sentinels. - * These functions return the sentinel value for the Langevin/Brownian - * parameters, indicating that they have not been set yet. +/** + * @brief Value for unset friction coefficient. + * Sentinel value for the Langevin/Brownian parameters, + * indicating that they have not been set yet. */ -/**@{*/ -constexpr double sentinel(double) { return -1.0; } -constexpr Utils::Vector3d sentinel(Utils::Vector3d) { - return {-1.0, -1.0, -1.0}; -} -/**@}*/ +#ifdef PARTICLE_ANISOTROPY +constexpr Thermostat::GammaType gamma_sentinel{{-1.0, -1.0, -1.0}}; +#else +constexpr Thermostat::GammaType gamma_sentinel{-1.0}; +#endif +/** + * @brief Value for a null friction coefficient. + */ +#ifdef PARTICLE_ANISOTROPY +constexpr Thermostat::GammaType gamma_null{{0.0, 0.0, 0.0}}; +#else +constexpr Thermostat::GammaType gamma_null{0.0}; +#endif } // namespace /************************************************ @@ -85,6 +96,30 @@ extern double temperature; /** True if the thermostat should act on virtual particles. */ extern bool thermo_virtual; +#ifdef THERMOSTAT_PER_PARTICLE +inline auto const & +handle_particle_gamma(Thermostat::GammaType const &particle_gamma, + Thermostat::GammaType const &default_gamma) { + return particle_gamma >= gamma_null ? particle_gamma : default_gamma; +} +#endif + +inline auto +handle_particle_anisotropy(Particle const &p, + Thermostat::GammaType const &gamma_body) { +#ifdef PARTICLE_ANISOTROPY + auto const aniso_flag = + (gamma_body[0] != gamma_body[1]) || (gamma_body[1] != gamma_body[2]); + const Utils::Matrix gamma_matrix = + boost::qvm::diag_mat(gamma_body); + auto const gamma_space = + aniso_flag ? convert_body_to_space(p, gamma_matrix) : gamma_matrix; + return gamma_space; +#else + return gamma_body; +#endif +} + /************************************************ * parameter structs ************************************************/ @@ -142,9 +177,9 @@ struct LangevinThermostat : public BaseThermostat { /** @name Parameters */ /**@{*/ /** Translational friction coefficient @f$ \gamma_{\text{trans}} @f$. */ - GammaType gamma = sentinel(GammaType{}); + GammaType gamma = gamma_sentinel; /** Rotational friction coefficient @f$ \gamma_{\text{rot}} @f$. */ - GammaType gamma_rotation = sentinel(GammaType{}); + GammaType gamma_rotation = gamma_sentinel; /**@}*/ /** @name Prefactors */ /**@{*/ @@ -219,9 +254,9 @@ struct BrownianThermostat : public BaseThermostat { /** @name Parameters */ /**@{*/ /** Translational friction coefficient @f$ \gamma_{\text{trans}} @f$. */ - GammaType gamma = sentinel(GammaType{}); + GammaType gamma = gamma_sentinel; /** Rotational friction coefficient @f$ \gamma_{\text{rot}} @f$. */ - GammaType gamma_rotation = sentinel(GammaType{}); + GammaType gamma_rotation = gamma_sentinel; /**@}*/ /** @name Prefactors */ /**@{*/ @@ -230,13 +265,13 @@ struct BrownianThermostat : public BaseThermostat { * @f$ D_{\text{trans}} = k_B T/\gamma_{\text{trans}} @f$ * the translational diffusion coefficient. */ - GammaType sigma_pos = sentinel(GammaType{}); + GammaType sigma_pos = gamma_sentinel; /** Rotational noise standard deviation. * Stores @f$ \sqrt{2D_{\text{rot}}} @f$ with * @f$ D_{\text{rot}} = k_B T/\gamma_{\text{rot}} @f$ * the rotational diffusion coefficient. */ - GammaType sigma_pos_rotation = sentinel(GammaType{}); + GammaType sigma_pos_rotation = gamma_sentinel; /** Translational velocity noise standard deviation. * Stores @f$ \sqrt{k_B T} @f$. */ diff --git a/src/core/thermostats/langevin_inline.hpp b/src/core/thermostats/langevin_inline.hpp index bf9b672eafb..8ee57aeacf0 100644 --- a/src/core/thermostats/langevin_inline.hpp +++ b/src/core/thermostats/langevin_inline.hpp @@ -50,16 +50,13 @@ friction_thermo_langevin(LangevinThermostat const &langevin, Particle const &p, } // Determine prefactors for the friction and the noise term - Thermostat::GammaType pref_friction = langevin.pref_friction; - Thermostat::GammaType pref_noise = langevin.pref_noise; #ifdef THERMOSTAT_PER_PARTICLE - // override default if particle-specific gamma - if (p.gamma() >= Thermostat::GammaType{}) { - auto const gamma = - p.gamma() >= Thermostat::GammaType{} ? p.gamma() : langevin.gamma; - pref_friction = -gamma; - pref_noise = LangevinThermostat::sigma(kT, time_step, gamma); - } + auto const gamma = handle_particle_gamma(p.gamma(), langevin.gamma); + auto const pref_friction = -gamma; + auto const pref_noise = LangevinThermostat::sigma(kT, time_step, gamma); +#else + auto const pref_friction = langevin.pref_friction; + auto const pref_noise = langevin.pref_noise; #endif // THERMOSTAT_PER_PARTICLE // Get effective velocity in the thermostatting @@ -70,26 +67,9 @@ friction_thermo_langevin(LangevinThermostat const &langevin, Particle const &p, #else auto const &velocity = p.v(); #endif // ENGINE -#ifdef PARTICLE_ANISOTROPY - // Particle frictional isotropy check - auto const aniso_flag = (pref_friction[0] != pref_friction[1]) || - (pref_friction[1] != pref_friction[2]); - // In case of anisotropic particle: body-fixed reference frame. Otherwise: - // lab-fixed reference frame. - const Utils::Matrix fric_mat = - boost::qvm::diag_mat(pref_friction); - const Utils::Matrix noise_mat = - boost::qvm::diag_mat(pref_noise); - - auto const friction_op = - aniso_flag ? convert_body_to_space(p, fric_mat) : fric_mat; - auto const noise_op = - aniso_flag ? convert_body_to_space(p, noise_mat) : noise_mat; -#else - auto const &friction_op = pref_friction; - auto const &noise_op = pref_noise; -#endif // PARTICLE_ANISOTROPY + auto const friction_op = handle_particle_anisotropy(p, pref_friction); + auto const noise_op = handle_particle_anisotropy(p, pref_noise); return friction_op * velocity + noise_op * Random::noise_uniform( @@ -111,23 +91,19 @@ friction_thermo_langevin_rotation(LangevinThermostat const &langevin, Particle const &p, double time_step, double kT) { - auto pref_friction = -langevin.gamma_rotation; - auto pref_noise = langevin.pref_noise_rotation; - #ifdef THERMOSTAT_PER_PARTICLE - // override default if particle-specific gamma - if (p.gamma_rot() >= Thermostat::GammaType{}) { - auto const gamma = p.gamma_rot() >= Thermostat::GammaType{} - ? p.gamma_rot() - : langevin.gamma_rotation; - pref_friction = -gamma; - pref_noise = LangevinThermostat::sigma(kT, time_step, gamma); - } + auto const gamma = + handle_particle_gamma(p.gamma_rot(), langevin.gamma_rotation); + auto const pref_friction = gamma; + auto const pref_noise = LangevinThermostat::sigma(kT, time_step, gamma); +#else + auto const pref_friction = langevin.gamma_rotation; + auto const pref_noise = langevin.pref_noise_rotation; #endif // THERMOSTAT_PER_PARTICLE auto const noise = Random::noise_uniform( langevin.rng_counter(), langevin.rng_seed(), p.id()); - return hadamard_product(pref_friction, p.omega()) + + return -hadamard_product(pref_friction, p.omega()) + hadamard_product(pref_noise, noise); } diff --git a/src/core/unit_tests/lb_particle_coupling_test.cpp b/src/core/unit_tests/lb_particle_coupling_test.cpp index bf20400e98e..40df3a57118 100644 --- a/src/core/unit_tests/lb_particle_coupling_test.cpp +++ b/src/core/unit_tests/lb_particle_coupling_test.cpp @@ -183,9 +183,13 @@ BOOST_AUTO_TEST_CASE(rng) { BOOST_REQUIRE(lb_particle_coupling.rng_counter_coupling); BOOST_CHECK_EQUAL(lb_lbcoupling_get_rng_state(), 17); BOOST_CHECK(not lb_lbcoupling_is_seed_required()); - auto const step1_random1 = coupling.get_noise_term(1); - auto const step1_random2 = coupling.get_noise_term(4); - auto const step1_random2_try2 = coupling.get_noise_term(4); + Particle test_partcl_1{}; + test_partcl_1.id() = 1; + Particle test_partcl_2{}; + test_partcl_2.id() = 4; + auto const step1_random1 = coupling.get_noise_term(test_partcl_1); + auto const step1_random2 = coupling.get_noise_term(test_partcl_2); + auto const step1_random2_try2 = coupling.get_noise_term(test_partcl_2); BOOST_REQUIRE(step1_random1.norm() > 1e-10); BOOST_CHECK(step1_random1 != step1_random2); BOOST_CHECK(step1_random2 == step1_random2_try2); @@ -196,13 +200,14 @@ BOOST_AUTO_TEST_CASE(rng) { BOOST_REQUIRE(lb_particle_coupling.rng_counter_coupling); BOOST_CHECK_EQUAL(lb_lbcoupling_get_rng_state(), 18); - auto const step2_random1 = coupling.get_noise_term(1); - auto const step2_random2 = coupling.get_noise_term(4); + auto const step2_random1 = coupling.get_noise_term(test_partcl_1); + auto const step2_random2 = coupling.get_noise_term(test_partcl_2); BOOST_CHECK(step1_random1 != step2_random1); BOOST_CHECK(step1_random1 != step2_random2); LB::ParticleCoupling coupling_unthermalized{true, params.time_step, 0.}; - auto const step3_norandom = coupling_unthermalized.get_noise_term(4); + auto const step3_norandom = + coupling_unthermalized.get_noise_term(test_partcl_2); BOOST_CHECK((step3_norandom == Utils::Vector3d{0., 0., 0.})); } @@ -318,7 +323,7 @@ BOOST_DATA_TEST_CASE(particle_coupling, bdata::make(kTs), kT) { lb_lbcoupling_set_gamma(gamma); Particle p{}; LB::ParticleCoupling coupling{false, params.time_step}; - auto expected = coupling.get_noise_term(p.id()); + auto expected = coupling.get_noise_term(p); #ifdef ENGINE p.swimming().swimming = true; p.swimming().v_swim = 2.; @@ -379,11 +384,11 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, #endif LB::ParticleCoupling coupling{thermo_virtual, params.time_step}; - auto expected = coupling.get_noise_term(pid); auto const p_opt = copy_particle_to_head_node(comm, pid); -#if defined(ENGINE) or defined(LB_ELECTROHYDRODYNAMICS) + auto expected = Utils::Vector3d{}; if (rank == 0) { auto const &p = *p_opt; + expected += coupling.get_noise_term(p); #ifdef ENGINE expected += gamma * p.swimming().v_swim * p.calc_director(); #endif @@ -391,7 +396,6 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, expected += gamma * p.mu_E(); #endif } -#endif // defined(ENGINE) or defined(LB_ELECTROHYDRODYNAMICS) boost::mpi::broadcast(comm, expected, 0); auto const p_pos = first_lb_node + Utils::Vector3d::broadcast(0.5); set_particle_pos(pid, p_pos); @@ -512,9 +516,10 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, } lb_particle_coupling.rng_counter_coupling = {}; if (kT == 0.) { - BOOST_CHECK_EQUAL(coupling.get_noise_term(0).norm(), 0.); + BOOST_CHECK_EQUAL(coupling.get_noise_term(Particle{}).norm(), 0.); } else { - BOOST_CHECK_THROW(coupling.get_noise_term(0), std::runtime_error); + BOOST_CHECK_THROW(coupling.get_noise_term(Particle{}), + std::runtime_error); } } } diff --git a/testsuite/python/lb_thermostat.py b/testsuite/python/lb_thermostat.py index 33b3df81e00..029963bdab1 100644 --- a/testsuite/python/lb_thermostat.py +++ b/testsuite/python/lb_thermostat.py @@ -29,16 +29,16 @@ distribution. """ -KT = 0.9 +KT = 0.2 AGRID = 0.8 node_volume = AGRID**3 -VISC = 6 -DENS = 1.7 +KIN_VISC = 0.3 +DENS = 0.8 TIME_STEP = 0.005 -GAMMA = 2 + LB_PARAMS = {'agrid': AGRID, 'density': DENS, - 'kinematic_viscosity': VISC, + 'kinematic_viscosity': KIN_VISC, 'tau': TIME_STEP, 'kT': KT, 'seed': 123} @@ -46,65 +46,163 @@ class LBThermostatCommon(thermostats_common.ThermostatsCommon): - """Check the LB thermostat.""" + """ + Check the LB thermostat. + All temperature checks rely on https://en.wikipedia.org/wiki/Thermal_velocity + with the root-mean-squared-velocity. + The temperature check together with the friction test + ensure that the noise amplitude is correct. + """ system = espressomd.System(box_l=[AGRID * 12] * 3) system.time_step = TIME_STEP system.cell_system.skin = 0.4 * AGRID + # use "small" particles such that radius=hydrodynamic radius + global_gamma = 0.1 * AGRID * 6 * np.pi * KIN_VISC * DENS + partcl_gamma = 0.05 * AGRID * 6 * np.pi * KIN_VISC * DENS + # relaxation time 0.2 sim time units + partcl_mass = 0.2 * partcl_gamma + np.random.seed(41) def setUp(self): self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params) self.system.actors.add(self.lbf) - self.system.thermostat.set_lb(LB_fluid=self.lbf, seed=5, gamma=GAMMA) + self.system.thermostat.set_lb( + LB_fluid=self.lbf, seed=5, gamma=self.global_gamma) def tearDown(self): self.system.actors.clear() self.system.thermostat.turn_off() + self.system.part.clear() + + def get_lb_kT(self, lbf): + nodes_mass = lbf[:, :, :].density * node_volume + nodes_vel_sq = np.sum(np.square(lbf[:, :, :].velocity), axis=3) + return np.mean(nodes_mass * nodes_vel_sq) / 3. - def get_lb_momentum(self, lbf): - nodes_den = lbf[:, :, :].density - nodes_vel = np.sum(np.square(lbf[:, :, :].velocity), axis=3) - return np.multiply(nodes_den, nodes_vel) + def get_lb_velocity(self, lbf): + return np.mean(lbf[:, :, :].velocity, axis=(0, 1, 2)) def test_fluid(self): self.system.integrator.run(100) - fluid_temps = [] + fluid_kTs = [] for _ in range(100): - lb_momentum = self.get_lb_momentum(self.lbf) - fluid_temps.append(np.average(lb_momentum) * node_volume) + fluid_kTs.append(self.get_lb_kT(self.lbf)) self.system.integrator.run(3) - fluid_temp = np.average(fluid_temps) / 3 - self.assertAlmostEqual(fluid_temp, KT, delta=0.05) + fluid_kT = np.mean(fluid_kTs) + np.testing.assert_allclose(fluid_kT, KT, rtol=0.05) - def test_with_particles(self): - self.system.part.add( - pos=np.random.random((100, 3)) * self.system.box_l) - self.system.integrator.run(120) - N = len(self.system.part) - loops = 250 - v_particles = np.zeros((loops, N, 3)) - fluid_temps = [] + def check_partcl_temp(self, partcl_vel): + partcl_vel_rms = np.sqrt( + np.mean(np.linalg.norm(partcl_vel, axis=2)**2)) - for i in range(loops): - self.system.integrator.run(3) - if i % 10 == 0: - lb_momentum = self.get_lb_momentum(self.lbf) - fluid_temps.append(np.average(lb_momentum) * node_volume) - v_particles[i] = self.system.part.all().v - fluid_temp = np.average(fluid_temps) / 3. + np.testing.assert_allclose( + np.mean(partcl_vel), 0, atol=0.05 * partcl_vel_rms) + np.testing.assert_allclose(partcl_vel_rms**2 * self.partcl_mass / 3., + KT, rtol=0.05) - np.testing.assert_allclose(np.average(v_particles), 0, atol=0.033) - np.testing.assert_allclose(np.var(v_particles), KT, atol=0.033) - - minmax = 3 + vel_range = 2 * partcl_vel_rms n_bins = 7 - error_tol = 0.016 self.check_velocity_distribution( - v_particles.reshape((-1, 3)), minmax, n_bins, error_tol, KT) - - np.testing.assert_allclose(fluid_temp, KT, atol=5e-3) + partcl_vel.reshape((-1, 3)), vel_range, n_bins, 0.02, KT, + mass=self.partcl_mass) + + def test_temperature_with_particles(self): + n_partcls_per_type = 50 + partcls_global_gamma = self.system.part.add( + pos=np.random.random((n_partcls_per_type, 3)) * self.system.box_l, + mass=n_partcls_per_type * [self.partcl_mass]) + partcls_per_part_gamma = self.system.part.add( + pos=np.random.random((n_partcls_per_type, 3)) * self.system.box_l, + mass=n_partcls_per_type * [self.partcl_mass], + gamma=n_partcls_per_type * [self.partcl_gamma]) + self.system.integrator.run(50) + n_samples = 400 + vel_global_gamma = np.zeros((n_samples, n_partcls_per_type, 3)) + vel_per_part_gamma = np.zeros_like(vel_global_gamma) + fluid_kTs = [] + + for i in range(n_samples): + self.system.integrator.run(3) + if i % 10 == 0: + fluid_kTs.append(self.get_lb_kT(self.lbf)) + vel_global_gamma[i] = partcls_global_gamma.v + vel_per_part_gamma[i] = partcls_per_part_gamma.v + fluid_kT = np.mean(fluid_kTs) + + self.check_partcl_temp(vel_global_gamma) + self.check_partcl_temp(vel_per_part_gamma) + + np.testing.assert_allclose(fluid_kT, KT, rtol=0.03) + + def test_friction(self): + """apply force and measure if the average velocity matches expectation""" + + # large force to get better signal to noise ratio + ext_force = np.array([1.2, 2.1, 1.1]) + n_partcls_each_type = 50 + partcls_global_gamma = self.system.part.add( + pos=np.random.random((n_partcls_each_type, 3)) * self.system.box_l, + ext_force=n_partcls_each_type * [ext_force], + mass=n_partcls_each_type * [self.partcl_mass]) + partcls_per_partcl_gamma = self.system.part.add( + pos=np.random.random((n_partcls_each_type, 3)) * self.system.box_l, + ext_force=n_partcls_each_type * [ext_force], + mass=n_partcls_each_type * [self.partcl_mass], + gamma=n_partcls_each_type * [self.partcl_gamma]) + + # add counterforce to fluid such that velocities cannot increase + # indefinitely and we get a steady state relative velocity + total_force_partcls = ext_force * 2 * n_partcls_each_type + lb_volume = np.prod(self.system.box_l) + self.lbf.ext_force_density = -total_force_partcls / lb_volume + + # let particles and fluid accelerate towards steady state before + # measuring + self.system.integrator.run(int(0.8 / self.system.time_step)) + startpos_global_gamma = np.copy(partcls_global_gamma.pos) + startpos_per_p_gamma = np.copy(partcls_per_partcl_gamma.pos) + + run_time = 0.3 + n_steps = int(run_time / self.system.time_step) + steps_per_sample = 10 + lb_vels = [] + for _ in range(int(n_steps / steps_per_sample)): + self.system.integrator.run(steps_per_sample) + lb_vels.append(self.get_lb_velocity(self.lbf)) + lb_vel = np.mean(lb_vels, axis=0) + + vel_global_gamma = np.mean( + partcls_global_gamma.pos - startpos_global_gamma, + axis=0) / run_time + vel_per_part_gamma = np.mean( + partcls_per_partcl_gamma.pos - startpos_per_p_gamma, + axis=0) / run_time + # get velocity relative to the fluid + vel_global_gamma -= lb_vel + vel_per_part_gamma -= lb_vel + + # average over 3 spatial directions (isotropic friction) + global_gamma_measured = np.mean(ext_force / vel_global_gamma) + per_partcl_gamma_measured = np.mean(ext_force / vel_per_part_gamma) + + np.testing.assert_allclose( + global_gamma_measured, + self.global_gamma, + rtol=0.05) + np.testing.assert_allclose( + per_partcl_gamma_measured, + self.partcl_gamma, + rtol=0.05) + + @utx.skipIfMissingFeatures(["PARTICLE_ANISOTROPY", + "THERMOSTAT_PER_PARTICLE"]) + def test_exceptions(self): + self.system.part.add(pos=[0., 0., 0.], gamma=[1., 2., 3.], id=2) + with self.assertRaisesRegex(Exception, r"ERROR: anisotropic particle \(id 2\) coupled to LB"): + self.system.integrator.run(1) @utx.skipIfMissingFeatures(["WALBERLA"]) @@ -126,7 +224,7 @@ class LBWalberlaThermostatSinglePrecision(LBThermostatCommon, ut.TestCase): # @utx.skipIfMissingGPU() -# @utx.skipIfMissingFeatures(["WALBERLA"]) +# @utx.skipIfMissingFeatures(["WALBERLA", "CUDA"]) # class LBWalberlaGPUThermostat(LBThermostatCommon, ut.TestCase): # """Test for the GPU implementation of the LB.""" diff --git a/testsuite/python/thermostats_common.py b/testsuite/python/thermostats_common.py index cc2c653a0d5..f0457d4d131 100644 --- a/testsuite/python/thermostats_common.py +++ b/testsuite/python/thermostats_common.py @@ -23,31 +23,33 @@ import espressomd.observables -def single_component_maxwell(x1, x2, kT): +def single_component_maxwell(x1, x2, kT, mass=1.): """Integrate the probability density from x1 to x2 using the trapezoidal rule""" x = np.linspace(x1, x2, 1000) - return np.trapz(np.exp(-x**2 / (2. * kT)), x) / np.sqrt(2. * np.pi * kT) + maxwell_distr = np.exp(- mass * x**2 / (2. * kT)) * \ + np.sqrt(mass / (2. * np.pi * kT)) + return np.trapz(maxwell_distr, x=x) class ThermostatsCommon: """Tests the velocity distribution created by a thermostat.""" - def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT): + def check_velocity_distribution( + self, vel, minmax, n_bins, error_tol, kT, mass=1.): """Check the recorded particle distributions in velocity against a histogram with n_bins bins. Drop velocities outside minmax. Check individual histogram bins up to an accuracy of error_tol against the analytical result for kT.""" for i in range(3): - hist = np.histogram( + hist, bin_edges = np.histogram( vel[:, i], range=(-minmax, minmax), bins=n_bins, density=False) - data = hist[0] / float(vel.shape[0]) - bins = hist[1] + hist = hist / len(vel) expected = [] for j in range(n_bins): expected.append(single_component_maxwell( - bins[j], bins[j + 1], kT)) - np.testing.assert_allclose(data[:n_bins], expected, atol=error_tol) + bin_edges[j], bin_edges[j + 1], kT, mass=mass)) + np.testing.assert_allclose(hist, expected, atol=error_tol) def test_00_verify_single_component_maxwell(self): """Verifies the normalization of the analytical expression."""