Fix EKSpecies setters

Author: jngrad

maintainer/walberla_kernels/generate_ek_kernels.py

@@ -56,7 +56,7 @@ def replace_macros(filename: str) -> None:
         content = f.read()
         f.seek(0)
         f.truncate(0)
-        # replace replace_getData_with_uncheckedFastGetData
+        # replace getData with uncheckedFastGetData
         content = content.replace("block->getData<IndexVectors>(indexVectorID);",
                                   "block->uncheckedFastGetData<IndexVectors>(indexVectorID);")
         # remove dummy assignment
@@ -154,7 +154,8 @@ def replace_macros(filename: str) -> None:
             block_offset=block_offsets if fluctuation else None,
             **params)
 
-    # the subsitution for field reads is necessary, because otherwise there are "ResolvedFieldAccess" nodes that fail in the code generation
+    # the substitution for field reads is necessary, because otherwise there are
+    # "ResolvedFieldAccess" nodes that fail in the code generation
     flux_advection = ps.AssignmentCollection(ek.flux_advection())
     flux_advection = ps.simp.add_subexpressions_for_field_reads(flux_advection)
 

src/script_interface/walberla/EKSpecies.hpp

@@ -36,7 +36,9 @@
 #include <utils/math/int_pow.hpp>
 
 #include <memory>
+#include <stdexcept>
 #include <string>
+#include <unordered_map>
 
 namespace ScriptInterface::walberla {
 
@@ -66,7 +68,13 @@ class EKSpecies : public LatticeModel<::EKinWalberlaBase, EKVTKHandle> {
           [this]() { return m_instance->get_diffusion() / m_conv_diffusion; }},
          {"kT",
           [this](Variant const &v) {
-            m_instance->set_kT(get_value<double>(v) * m_conv_energy);
+            context()->parallel_try_catch([&]() {
+              auto const kT = get_value<double>(v);
+              if (kT < 0.) {
+                throw std::domain_error("Parameter 'kT' must be >= 0");
+              }
+              m_instance->set_kT(kT * m_conv_energy);
+            });
           },
           [this]() { return m_instance->get_kT() / m_conv_energy; }},
          {"valency",

src/walberla_bridge/include/walberla_bridge/electrokinetics/EKinWalberlaBase.hpp

@@ -140,7 +140,7 @@ class EKinWalberlaBase : public LatticeModel {
   [[nodiscard]] virtual Utils::Vector3d get_ext_efield() const noexcept = 0;
   [[nodiscard]] virtual bool is_double_precision() const noexcept = 0;
   [[nodiscard]] virtual bool is_thermalized() const noexcept = 0;
-  [[nodiscard]] virtual uint get_seed() const noexcept = 0;
+  [[nodiscard]] virtual unsigned int get_seed() const noexcept = 0;
   [[nodiscard]] virtual std::optional<uint64_t> get_rng_state() const = 0;
 
   virtual void set_diffusion(double diffusion) = 0;

src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_walberla_init.hpp

@@ -36,7 +36,7 @@ new_ek_walberla(std::shared_ptr<LatticeWalberla> const &lattice,
                 double diffusion, double kT, double valency,
                 Utils::Vector3d ext_efield, double density, bool advection,
                 bool friction_coupling, bool single_precision, bool thermalized,
-                uint seed);
+                unsigned int seed);
 
 std::shared_ptr<EKReactionBase>
 new_ek_reaction_bulk(std::shared_ptr<LatticeWalberla> const &lattice,

src/walberla_bridge/src/electrokinetics/EKinWalberlaImpl.hpp

@@ -109,7 +109,7 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
   Utils::Vector3d m_ext_efield;
   bool m_advection;
   bool m_friction_coupling;
-  uint m_seed;
+  unsigned int m_seed;
 
 protected:
   // Block data access handles
@@ -177,7 +177,7 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
   EKinWalberlaImpl(std::shared_ptr<LatticeWalberla> lattice, double diffusion,
                    double kT, double valency, Utils::Vector3d const &ext_efield,
                    double density, bool advection, bool friction_coupling,
-                   bool thermalized, uint seed)
+                   bool thermalized, unsigned int seed)
       : m_diffusion(FloatType_c(diffusion)), m_kT(FloatType_c(kT)),
         m_valency(FloatType_c(valency)), m_ext_efield(ext_efield),
         m_advection(advection), m_friction_coupling(friction_coupling),
@@ -247,7 +247,9 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
     return static_cast<bool>(
         std::get_if<DiffusiveFluxKernelThermalized>(&*m_diffusive_flux));
   }
-  [[nodiscard]] uint get_seed() const noexcept override { return m_seed; }
+  [[nodiscard]] unsigned int get_seed() const noexcept override {
+    return m_seed;
+  }
   [[nodiscard]] std::optional<uint64_t> get_rng_state() const override {
     auto const kernel =
         std::get_if<DiffusiveFluxKernelThermalized>(&*m_diffusive_flux);
@@ -259,30 +261,31 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
 
   void set_diffusion(double diffusion) override {
     m_diffusion = FloatType_c(diffusion);
-
-    auto lambda = [m_diffusion = m_diffusion](auto &kernel) {
+    auto visitor = [m_diffusion = m_diffusion](auto &kernel) {
       kernel.D_ = m_diffusion;
     };
-
-    std::visit(lambda, *m_diffusive_flux);
-    std::visit(lambda, *m_diffusive_flux_electrostatic);
+    std::visit(visitor, *m_diffusive_flux);
+    std::visit(visitor, *m_diffusive_flux_electrostatic);
   }
+
   void set_kT(double kT) override {
     m_kT = FloatType_c(kT);
-
     std::visit([m_kT = m_kT](auto &kernel) { kernel.kT_ = m_kT; },
                *m_diffusive_flux_electrostatic);
   }
+
   void set_valency(double valency) override {
     m_valency = FloatType_c(valency);
-
     std::visit([m_valency = m_valency](auto &kernel) { kernel.z_ = m_valency; },
                *m_diffusive_flux_electrostatic);
   }
+
   void set_advection(bool advection) override { m_advection = advection; }
+
   void set_friction_coupling(bool friction_coupling) override {
     m_friction_coupling = friction_coupling;
   }
+
   void set_rng_state(uint64_t counter) override {
     auto const kernel =
         std::get_if<DiffusiveFluxKernelThermalized>(&*m_diffusive_flux);
@@ -298,6 +301,7 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
     kernel->time_step_ = static_cast<uint32_t>(counter);
     kernel_electrostatic->time_step_ = static_cast<uint32_t>(counter);
   }
+
   void set_ext_efield(Utils::Vector3d const &field) override {
     m_ext_efield = field;
 
@@ -330,7 +334,7 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
             std::move(kernel_electrostatic));
   }
 
-  void set_diffusion_kernels(uint seed) {
+  void set_diffusion_kernels(unsigned int seed) {
     auto const grid_dim = get_lattice().get_grid_dimensions();
 
     auto kernel = DiffusiveFluxKernelThermalized(
@@ -412,9 +416,8 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
   }
 
   void kernel_diffusion_electrostatic(const std::size_t &potential_id) {
-    auto const potential_id_blockdata = BlockDataID(potential_id);
-    std::visit([phiID = potential_id_blockdata](
-                   auto &kernel) { kernel.phiID = phiID; },
+    auto const phiID = BlockDataID(potential_id);
+    std::visit([phiID](auto &kernel) { kernel.phiID = phiID; },
                *m_diffusive_flux_electrostatic);
 
     for (auto &block : *m_lattice->get_blocks()) {

src/walberla_bridge/src/electrokinetics/ek_walberla_init.cpp

@@ -37,7 +37,7 @@ new_ek_walberla(std::shared_ptr<LatticeWalberla> const &lattice,
                 double diffusion, double kT, double valency,
                 Utils::Vector3d ext_efield, double density, bool advection,
                 bool friction_coupling, bool single_precision, bool thermalized,
-                uint seed) {
+                unsigned int seed) {
   if (single_precision) {
     return std::make_shared<EKinWalberlaImpl<13, float>>(
         lattice, diffusion, kT, valency, ext_efield, density, advection,

testsuite/python/ek_interface.py

@@ -100,6 +100,9 @@ def test_ek_species(self):
         ek_species = self.make_default_ek_species(thermalized=True, seed=42)
         self.assertTrue(ek_species.thermalized)
         self.assertEqual(ek_species.seed, 42)
+        self.assertEqual(ek_species.rng_state, 0)
+        ek_species.rng_state = 5
+        self.assertEqual(ek_species.rng_state, 5)
 
     def check_ek_species_properties(self, species):
         agrid = self.params["agrid"]
@@ -199,6 +202,15 @@ def test_ek_none_solver(self):
         self.assertIsInstance(self.system.ekcontainer.solver,
                               espressomd.electrokinetics.EKNone)
 
+    def test_ek_species_exceptions(self):
+        ek_species = self.make_default_ek_species()
+        with self.assertRaisesRegex(ValueError, "Parameter 'kT' must be >= 0"):
+            ek_species.kT = -0.4
+        with self.assertRaisesRegex(ValueError, "Parameter 'rng_state' must be >= 0"):
+            ek_species.rng_state = -2
+        with self.assertRaisesRegex(RuntimeError, "This EK instance is unthermalized"):
+            ek_species.rng_state = 5
+
     def test_ek_solver_exceptions(self):
         ek_solver = self.system.ekcontainer.solver
         ek_species = self.make_default_ek_species()