diff --git a/doc/tutorials/electrokinetics/electrokinetics.ipynb b/doc/tutorials/electrokinetics/electrokinetics.ipynb index 32e5f952b8..3fe674686e 100644 --- a/doc/tutorials/electrokinetics/electrokinetics.ipynb +++ b/doc/tutorials/electrokinetics/electrokinetics.ipynb @@ -1042,7 +1042,7 @@ " lattice=lattice, density=DENSITY_FLUID, kinematic_viscosity=VISCOSITY_KINEMATIC,\n", " tau=TAU, ext_force_density=EXT_FORCE_DENSITY, kT=KT, seed=42)\n", "system.lb = lbf\n", - "system.thermostat.set_lb(LB_fluid=lbf, seed=42)\n", + "system.thermostat.set_lb(LB_fluid=lbf, seed=42, gamma=0.)\n", "eksolver = espressomd.electrokinetics.EKNone(lattice=lattice)\n", "system.ekcontainer = espressomd.electrokinetics.EKContainer(tau=TAU, solver=eksolver)" ] diff --git a/samples/lb_profile.py b/samples/lb_profile.py index 763bd7f95a..6b47778f29 100644 --- a/samples/lb_profile.py +++ b/samples/lb_profile.py @@ -45,7 +45,6 @@ agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.01, ext_force_density=[0, 0, 0.15], kT=0.0) system.lb = lb_fluid -system.thermostat.set_lb(LB_fluid=lb_fluid, seed=23) ctp = espressomd.math.CylindricalTransformationParameters( center=[5.0, 5.0, 0.0], axis=[0, 0, 1], diff --git a/src/core/lb/particle_coupling.cpp b/src/core/lb/particle_coupling.cpp index 648384f1f6..fad6bb64b8 100644 --- a/src/core/lb/particle_coupling.cpp +++ b/src/core/lb/particle_coupling.cpp @@ -291,18 +291,23 @@ void ParticleCoupling::kernel(std::vector const &particles) { #endif Utils::Vector3d force_on_particle = {}; if (coupling_mode == particle_force) { - auto &v_fluid = *it_interpolated_velocities; - if (m_box_geo.type() == BoxType::LEES_EDWARDS) { - // Account for the case where the interpolated velocity has been read - // from a ghost of the particle across the LE boundary (or vice verssa) - // Then the particle velocity is shifted by +,- the LE shear velocity - auto const vel_correction = lees_edwards_vel_shift( - *it_positions_velocity_coupling, p.pos(), m_box_geo); - v_fluid += vel_correction; +#ifndef THERMOSTAT_PER_PARTICLE + if (m_thermostat.gamma > 0.) +#endif + { + auto &v_fluid = *it_interpolated_velocities; + if (m_box_geo.type() == BoxType::LEES_EDWARDS) { + // Account for the case where the interpolated velocity has been read + // from a ghost of the particle across the LE boundary (or vice versa) + // Then the particle velocity is shifted by +,- the LE shear velocity + auto const vel_correction = lees_edwards_vel_shift( + *it_positions_velocity_coupling, p.pos(), m_box_geo); + v_fluid += vel_correction; + } + auto const drag_force = lb_drag_force(p, m_thermostat.gamma, v_fluid); + auto const random_force = get_noise_term(p); + force_on_particle = drag_force + random_force; } - auto const drag_force = lb_drag_force(p, m_thermostat.gamma, v_fluid); - auto const random_force = get_noise_term(p); - force_on_particle = drag_force + random_force; ++it_interpolated_velocities; ++it_positions_velocity_coupling; } diff --git a/src/core/thermostat.hpp b/src/core/thermostat.hpp index 121c6a4ebd..8d5781257c 100644 --- a/src/core/thermostat.hpp +++ b/src/core/thermostat.hpp @@ -97,6 +97,9 @@ inline bool are_kT_equal(double old_kT, double new_kT) { } auto const large_kT = (old_kT > new_kT) ? old_kT : new_kT; auto const small_kT = (old_kT > new_kT) ? new_kT : old_kT; + if (small_kT == 0.) { + return false; + } return (large_kT / small_kT - 1. < relative_tolerance); } } // namespace Thermostat diff --git a/src/script_interface/thermostat/thermostat.hpp b/src/script_interface/thermostat/thermostat.hpp index 57318778f1..bb63ad38a3 100644 --- a/src/script_interface/thermostat/thermostat.hpp +++ b/src/script_interface/thermostat/thermostat.hpp @@ -33,6 +33,7 @@ #include #include #include +#include #include #include @@ -92,7 +93,7 @@ class Interface : public AutoParameters, System::Leaf> { } virtual bool invalid_rng_state(VariantMap const ¶ms) const { - return (not params.count("seed") or is_none(params.at("seed"))) and + return (not params.contains("seed") or is_none(params.at("seed"))) and is_seed_required(); } @@ -101,14 +102,15 @@ class Interface : public AutoParameters, System::Leaf> { get_member_handle(::Thermostat::Thermostat &thermostat) = 0; void set_new_parameters(VariantMap const ¶ms) { - if (params.count("__check_rng_state") and invalid_rng_state(params)) { - context()->parallel_try_catch([]() { + context()->parallel_try_catch([&]() { + if (params.contains("__check_rng_state") and invalid_rng_state(params)) { throw std::invalid_argument("Parameter 'seed' is needed on first " "activation of the thermostat"); - }); - } + } + check_required_parameters(params); + }); for (auto const &key : get_parameter_insertion_order()) { - if (params.count(key)) { + if (params.contains(key)) { auto const &v = params.at(key); if (key == "is_active") { is_active = get_value(v); @@ -119,6 +121,17 @@ class Interface : public AutoParameters, System::Leaf> { } } + virtual std::span get_required_parameters() const = 0; + + void check_required_parameters(VariantMap const ¶ms) const { + for (auto const &required : get_required_parameters()) { + auto name = std::string(required); + if (not params.contains(name)) { + throw std::runtime_error("Parameter '" + name + "' is missing"); + } + } + } + protected: template auto make_autoparameter(T CoreThermostat::*member, char const *name) { @@ -169,7 +182,7 @@ class Interface : public AutoParameters, System::Leaf> { auto params = parameters; if (not is_seed_required()) { for (auto key : {std::string("seed"), std::string("philox_counter")}) { - if (params.count(key) == 0ul) { + if (not params.contains(key)) { params[key] = get_parameter(key); } } @@ -216,7 +229,7 @@ class Interface : public AutoParameters, System::Leaf> { } virtual std::optional extract_kT(VariantMap const ¶ms) const { - if (params.count("kT")) { + if (params.contains("kT")) { auto const value = get_value(params, "kT"); sanity_checks_positive(value, "kT"); return value; @@ -301,6 +314,11 @@ class Langevin : public Interface<::LangevinThermostat> { return thermostat.langevin; } + std::span get_required_parameters() const override { + static constexpr std::array names{std::string_view("gamma")}; + return names; + } + public: Langevin() { add_parameters({ @@ -319,7 +337,7 @@ class Langevin : public Interface<::LangevinThermostat> { Interface<::LangevinThermostat>::extend_parameters(parameters); #ifdef ROTATION // If gamma_rotation is not set explicitly, use the translational one. - if (params.count("gamma_rotation") == 0ul and params.count("gamma")) { + if (not params.contains("gamma_rotation") and params.contains("gamma")) { params["gamma_rotation"] = params.at("gamma"); } #endif // ROTATION @@ -333,6 +351,11 @@ class Brownian : public Interface<::BrownianThermostat> { return thermostat.brownian; } + std::span get_required_parameters() const override { + static constexpr std::array names{std::string_view("gamma")}; + return names; + } + public: Brownian() { add_parameters({ @@ -351,7 +374,7 @@ class Brownian : public Interface<::BrownianThermostat> { Interface<::BrownianThermostat>::extend_parameters(parameters); #ifdef ROTATION // If gamma_rotation is not set explicitly, use the translational one. - if (params.count("gamma_rotation") == 0ul and params.count("gamma")) { + if (not params.contains("gamma_rotation") and params.contains("gamma")) { params["gamma_rotation"] = params.at("gamma"); } #endif // ROTATION @@ -366,6 +389,12 @@ class IsotropicNpt : public Interface<::IsotropicNptThermostat> { return thermostat.npt_iso; } + std::span get_required_parameters() const override { + static constexpr std::array names{std::string_view("gamma0"), + std::string_view("gammav")}; + return names; + } + public: IsotropicNpt() { add_parameters({ @@ -418,10 +447,15 @@ class LBThermostat : public Interface<::LBThermostat> { protected: bool invalid_rng_state(VariantMap const ¶ms) const override { - return (not params.count("seed") or is_none(params.at("seed"))) and - params.count("__global_kT") and is_seed_required() and + return (not params.contains("seed") or is_none(params.at("seed"))) and + params.contains("__global_kT") and is_seed_required() and get_value(params, "__global_kT") != 0.; } + + std::span get_required_parameters() const override { + static constexpr std::array names{std::string_view("gamma")}; + return names; + } }; #endif // WALBERLA @@ -432,6 +466,10 @@ class DPDThermostat : public Interface<::DPDThermostat> { return thermostat.dpd; } + std::span get_required_parameters() const override { + return {}; + } + public: ::ThermostatFlags get_thermo_flag() const final { return THERMO_DPD; } }; @@ -444,6 +482,10 @@ class Stokesian : public Interface<::StokesianThermostat> { return thermostat.stokesian; } + std::span get_required_parameters() const override { + return {}; + } + public: ::ThermostatFlags get_thermo_flag() const final { return THERMO_SD; } }; @@ -455,6 +497,10 @@ class ThermalizedBond : public Interface<::ThermalizedBondThermostat> { return thermostat.thermalized_bond; } + std::span get_required_parameters() const override { + return {}; + } + public: ::ThermostatFlags get_thermo_flag() const final { return THERMO_BOND; } }; diff --git a/testsuite/python/lb.py b/testsuite/python/lb.py index 4e585d5f08..a16ff6c5a3 100644 --- a/testsuite/python/lb.py +++ b/testsuite/python/lb.py @@ -54,7 +54,6 @@ class LBTest: system.cell_system.skin = 1.0 if espressomd.gpu_available(): system.cuda_init_handle.call_method("set_device_id_per_rank") - interpolation = False n_nodes = system.cell_system.get_state()["n_nodes"] def setUp(self): @@ -75,7 +74,6 @@ def test_properties(self): # activated actor lbf = self.lb_class(kT=1.0, seed=42, **self.params, **self.lb_params) self.system.lb = lbf - self.system.thermostat.set_lb(LB_fluid=lbf, seed=1) self.assertTrue(lbf.is_active) self.check_properties(lbf) self.system.lb = None @@ -83,7 +81,6 @@ def test_properties(self): # deactivated actor lbf = self.lb_class(kT=1.0, seed=42, **self.params, **self.lb_params) self.system.lb = lbf - self.system.thermostat.set_lb(LB_fluid=lbf, seed=1) self.system.lb = None self.assertFalse(lbf.is_active) self.check_properties(lbf) @@ -383,7 +380,7 @@ def test_pressure_tensor_observable(self): lbf = self.lb_class(kT=1., seed=1, ext_force_density=[0, 0, 0], **self.params, **self.lb_params) system.lb = lbf - system.thermostat.set_lb(LB_fluid=lbf, seed=1) + system.thermostat.set_lb(LB_fluid=lbf, seed=1, gamma=1.) system.integrator.run(10) pressure_tensor = np.copy( np.mean(lbf[:, :, :].pressure_tensor, axis=(0, 1, 2))) diff --git a/testsuite/python/lb_interpolation.py b/testsuite/python/lb_interpolation.py index 96e93523b5..3a2bd9a491 100644 --- a/testsuite/python/lb_interpolation.py +++ b/testsuite/python/lb_interpolation.py @@ -128,26 +128,33 @@ def test_mach_limit_check(self): self.assertIsNone(self.lbf[0, 0, 0].boundary) def test_interpolated_force(self): + def sample(kernel): + for i, j, k in itertools.product(range(3), range(3), range(3)): + system.lb[:, :, :].velocity = lb_vel + p.pos = ( + (i + 0.5) * AGRID, + (j + 0.5) * AGRID, + (k + 0.5) * AGRID) + p.v = [0., 0., 0.] + system.integrator.run(1) + np.testing.assert_allclose(np.copy(p.f), kernel(i, j, k)) system = self.system - system.thermostat.set_lb(LB_fluid=system.lb, seed=42, gamma=1.) system.integrator.run(1) lb_vel = np.zeros([12, 12, 12, 3], dtype=float) - for i in range(12): - for j in range(12): - for k in range(12): - lb_vel[i, j, k] = 1e-3 * np.array([i + 1, j + 1, k + 1]) + for i, j, k in itertools.product(range(12), range(12), range(12)): + lb_vel[i, j, k] = 1e-3 * np.array([i + 1, j + 1, k + 1]) p = system.part.add(pos=3 * [1.5]) - for i in range(3): - for j in range(3): - for k in range(3): - system.lb[:, :, :].velocity = lb_vel - p.pos = ( - (i + 0.5) * AGRID, - (j + 0.5) * AGRID, - (k + 0.5) * AGRID) - p.v = [0., 0., 0.] - system.integrator.run(1) - np.testing.assert_allclose(np.copy(p.f), lb_vel[i, j, k]) + # enable particle coupling + system.thermostat.set_lb(LB_fluid=system.lb, seed=42, gamma=1.) + sample(kernel=lambda i, j, k: lb_vel[i, j, k]) + # disable particle coupling + system.thermostat.set_lb(LB_fluid=system.lb, seed=42, gamma=0.) + sample(kernel=lambda i, j, k: [0., 0., 0.]) + if espressomd.has_features("THERMOSTAT_PER_PARTICLE"): + # disable particle coupling globally & enable per-particle coupling + system.thermostat.set_lb(LB_fluid=system.lb, seed=42, gamma=0.) + p.gamma = 1.5 + sample(kernel=lambda i, j, k: 1.5 * lb_vel[i, j, k]) @utx.skipIfMissingFeatures(["WALBERLA"]) diff --git a/testsuite/python/lb_pressure_tensor.py b/testsuite/python/lb_pressure_tensor.py index d4fb48b231..0c0b4df3a6 100644 --- a/testsuite/python/lb_pressure_tensor.py +++ b/testsuite/python/lb_pressure_tensor.py @@ -46,7 +46,6 @@ class TestLBPressureTensor: @classmethod def tearDownClass(cls): cls.system.lb = None - cls.system.thermostat.turn_off() @classmethod def setUpClass(cls): diff --git a/testsuite/python/lb_thermostat.py b/testsuite/python/lb_thermostat.py index 77677bb96e..53de0ea8ca 100644 --- a/testsuite/python/lb_thermostat.py +++ b/testsuite/python/lb_thermostat.py @@ -208,6 +208,8 @@ def test_exceptions(self): with self.assertRaisesRegex(RuntimeError, r"set_lb\(\) got an unexpected keyword argument 'act_on_virtual'"): self.system.thermostat.set_lb( LB_fluid=self.lbf, act_on_virtual=False) + with self.assertRaisesRegex(RuntimeError, "Parameter 'gamma' is missing"): + self.system.thermostat.set_lb(LB_fluid=self.lbf) 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)