Introduce an event on_lees_edwards_change (#4977)

Fixes #4874

Description of changes:
- introduced a new event function on_lees_edwards_change() to be triggered when the LE BC is changed
- the order of instantiation of LB and LEbc no longer matters: collision kernel is automatically updated

.gitlab-ci.yml

@@ -435,7 +435,8 @@ check_sphinx:
   <<: *global_job_definition
   stage: additional_checks
   needs:
-    - cuda12-maxset
+    - job: cuda12-maxset
+      artifacts: true
   when: on_success
   script:
     - cd ${CI_PROJECT_DIR}/build
@@ -458,7 +459,8 @@ run_tutorials:
   <<: *global_job_definition
   stage: additional_checks
   needs:
-    - cuda12-maxset
+    - job: cuda12-maxset
+      artifacts: true
   when: on_success
   script:
     - cd ${CI_PROJECT_DIR}/build
@@ -486,7 +488,8 @@ run_doxygen:
   <<: *global_job_definition
   stage: additional_checks
   needs:
-    - cuda12-maxset
+    - job: cuda12-maxset
+      artifacts: true
   when: on_success
   only:
     - python
@@ -508,7 +511,8 @@ maxset_no_gpu:
   stage: additional_checks
   when: on_success
   needs:
-    - cuda12-maxset
+    - job: cuda12-maxset
+      artifacts: true
   script:
     - export CUDA_VISIBLE_DEVICES=""
     - cd ${CI_PROJECT_DIR}/build
@@ -524,15 +528,32 @@ maxset_3_cores:
   stage: additional_checks
   when: on_success
   needs:
-    - cuda12-maxset
+    - job: cuda12-maxset
+      artifacts: false
+  variables:
+     CC: 'gcc-12'
+     CXX: 'g++-12'
+     GCOV: 'gcov-12'
+     myconfig: 'maxset'
+     with_cuda: 'true'
+     with_coverage: 'false'
+     with_coverage_python: 'false'
+     srcdir: '${CI_PROJECT_DIR}'
+     with_scafacos: 'true'
+     with_walberla: 'true'
+     with_walberla_avx: 'true'
+     with_stokesian_dynamics: 'true'
+     with_caliper: 'true'
+     check_odd_only: 'true'
+     cmake_params: '-D ESPRESSO_TEST_NP=3'
   script:
-    - cd ${CI_PROJECT_DIR}/build
-    - cmake -D ESPRESSO_TEST_NP=3 .
-    - make -t && make check_unit_tests && make check_python_parallel_odd
+    - bash maintainer/CI/build_cmake.sh
+  timeout: 90m
   tags:
     - espresso
     - cuda
     - numa
+    - avx2
     - reuse-artifacts-same-arch
 
 status_success:

src/core/lb/LBNone.hpp

@@ -67,6 +67,8 @@ struct LBNone {
   void on_node_grid_change() const { throw NoLBActive{}; }
   void on_timestep_change() const { throw NoLBActive{}; }
   void on_temperature_change() const { throw NoLBActive{}; }
+  void on_lees_edwards_change() const { throw NoLBActive{}; }
+  void update_collision_model() const { throw NoLBActive{}; }
 };
 
 } // namespace LB

src/core/lb/LBWalberla.cpp

@@ -30,12 +30,15 @@
 #include "system/System.hpp"
 #include "thermostat.hpp"
 
+#include "lees_edwards/lees_edwards.hpp"
+
 #include <walberla_bridge/lattice_boltzmann/LBWalberlaBase.hpp>
 
 #include <utils/Vector.hpp>
 #include <utils/math/int_pow.hpp>
 
 #include <optional>
+#include <variant>
 
 namespace LB {
 
@@ -112,6 +115,45 @@ void LBWalberla::sanity_checks(System::System const &system) const {
                              system.get_time_step());
 }
 
+void LBWalberla::on_lees_edwards_change() { update_collision_model(); }
+
+void LBWalberla::update_collision_model() {
+  auto const energy_conversion =
+      Utils::int_pow<2>(lb_params->get_agrid() / lb_params->get_tau());
+  auto const kT = lb_fluid->get_kT() * energy_conversion;
+  auto const seed = lb_fluid->get_seed();
+  update_collision_model(*lb_fluid, *lb_params, kT, seed);
+}
+
+void LBWalberla::update_collision_model(LBWalberlaBase &lb,
+                                        LBWalberlaParams &params, double kT,
+                                        unsigned int seed) {
+  auto const &system = ::System::get_system();
+  auto le_protocol = system.lees_edwards->get_protocol();
+  if (le_protocol and
+      not std::holds_alternative<LeesEdwards::Off>(*le_protocol)) {
+    if (kT != 0.) {
+      throw std::runtime_error(
+          "Lees-Edwards LB doesn't support thermalization");
+    }
+    auto const &le_bc = system.box_geo->lees_edwards_bc();
+    auto lees_edwards_object = std::make_unique<LeesEdwardsPack>(
+        le_bc.shear_direction, le_bc.shear_plane_normal,
+        [&params, le_protocol, &system]() {
+          return get_pos_offset(system.get_sim_time(), *le_protocol) /
+                 params.get_agrid();
+        },
+        [&params, le_protocol, &system]() {
+          return get_shear_velocity(system.get_sim_time(), *le_protocol) *
+                 (params.get_tau() / params.get_agrid());
+        });
+    lb.set_collision_model(std::move(lees_edwards_object));
+  } else {
+    lb.set_collision_model(kT, seed);
+  }
+  lb.ghost_communication();
+}
+
 } // namespace LB
 
 #endif // WALBERLA

src/core/lb/LBWalberla.hpp

@@ -88,6 +88,11 @@ struct LBWalberla {
   }
   void on_timestep_change() const {}
   void on_temperature_change() const {}
+  void on_lees_edwards_change();
+  void update_collision_model();
+  static void update_collision_model(LBWalberlaBase &instance,
+                                     LBWalberlaParams &params, double kT,
+                                     unsigned int seed);
 };
 
 } // namespace LB

src/core/lb/Solver.cpp

@@ -140,6 +140,18 @@ void Solver::on_temperature_change() {
   }
 }
 
+void Solver::on_lees_edwards_change() {
+  if (impl->solver) {
+    std::visit([](auto &ptr) { ptr->on_lees_edwards_change(); }, *impl->solver);
+  }
+}
+
+void Solver::update_collision_model() {
+  if (impl->solver) {
+    std::visit([](auto &ptr) { ptr->update_collision_model(); }, *impl->solver);
+  }
+}
+
 bool Solver::is_gpu() const {
   check_solver(impl);
   return std::visit([](auto &ptr) { return ptr->is_gpu(); }, *impl->solver);

src/core/lb/Solver.hpp

@@ -168,6 +168,8 @@ struct Solver : public System::Leaf<Solver> {
   void on_cell_structure_change();
   void on_timestep_change();
   void on_temperature_change();
+  void on_lees_edwards_change();
+  void update_collision_model();
   void veto_boxl_change() const;
 
 private:

src/core/system/System.cpp

@@ -378,6 +378,8 @@ void System::on_observable_calc() {
   clear_particle_node();
 }
 
+void System::on_lees_edwards_change() { lb.on_lees_edwards_change(); }
+
 void System::update_local_geo() {
   *local_geo = LocalBox::make_regular_decomposition(
       box_geo->length(), ::communicator.calc_node_index(),

src/core/system/System.hpp

@@ -253,6 +253,7 @@ class System : public std::enable_shared_from_this<System> {
    *  initialized immediately (P3M etc.).
    */
   void on_observable_calc();
+  void on_lees_edwards_change();
   void veto_boxl_change(bool skip_particle_checks = false) const;
   /**@}*/
 

src/core/unit_tests/lb_particle_coupling_test.cpp

@@ -608,11 +608,13 @@ BOOST_AUTO_TEST_CASE(lb_exceptions) {
     BOOST_CHECK_THROW(lb.veto_kT(0.), NoLBActive);
     BOOST_CHECK_THROW(lb.lebc_sanity_checks(0u, 1u), NoLBActive);
     BOOST_CHECK_THROW(lb.propagate(), NoLBActive);
+    BOOST_CHECK_THROW(lb.update_collision_model(), NoLBActive);
     BOOST_CHECK_THROW(lb.on_cell_structure_change(), NoLBActive);
     BOOST_CHECK_THROW(lb.on_boxl_change(), NoLBActive);
     BOOST_CHECK_THROW(lb.on_node_grid_change(), NoLBActive);
     BOOST_CHECK_THROW(lb.on_timestep_change(), NoLBActive);
     BOOST_CHECK_THROW(lb.on_temperature_change(), NoLBActive);
+    BOOST_CHECK_THROW(lb.on_lees_edwards_change(), NoLBActive);
     BOOST_CHECK_THROW(lb_impl->get_density_at_pos(vec, true), NoLBActive);
     BOOST_CHECK_THROW(lb_impl->get_velocity_at_pos(vec, true), NoLBActive);
     BOOST_CHECK_THROW(lb_impl->add_force_at_pos(vec, vec), NoLBActive);

src/script_interface/lees_edwards/LeesEdwards.hpp

@@ -45,15 +45,16 @@ class LeesEdwards : public AutoParameters<LeesEdwards, System::Leaf> {
     add_parameters(
         {{"protocol",
           [this](Variant const &value) {
+            auto &system = get_system();
             if (is_none(value)) {
-              auto const &system = get_system();
-              context()->parallel_try_catch([&system]() {
+              context()->parallel_try_catch([&]() {
                 auto constexpr invalid_dir = LeesEdwardsBC::invalid_dir;
                 system.lb.lebc_sanity_checks(invalid_dir, invalid_dir);
               });
               m_protocol = nullptr;
               system.box_geo->set_lees_edwards_bc(LeesEdwardsBC{});
               m_lees_edwards->unset_protocol();
+              system.on_lees_edwards_change();
               return;
             }
             context()->parallel_try_catch([this]() {
@@ -67,8 +68,15 @@ class LeesEdwards : public AutoParameters<LeesEdwards, System::Leaf> {
                     "activation via set_boundary_conditions()");
               }
             });
-            m_protocol = get_value<std::shared_ptr<Protocol>>(value);
-            m_lees_edwards->set_protocol(m_protocol->protocol());
+            auto const protocol = get_value<std::shared_ptr<Protocol>>(value);
+            context()->parallel_try_catch([&]() {
+              try {
+                set_protocol(protocol);
+              } catch (...) {
+                set_protocol(m_protocol);
+                throw;
+              }
+            });
           },
           [this]() {
             if (m_protocol)
@@ -89,26 +97,35 @@ class LeesEdwards : public AutoParameters<LeesEdwards, System::Leaf> {
                          VariantMap const &params) override {
     if (name == "set_boundary_conditions") {
       context()->parallel_try_catch([this, &params]() {
-        auto const protocol = params.at("protocol");
-        if (is_none(protocol)) {
-          do_set_parameter("protocol", protocol);
+        auto const variant = params.at("protocol");
+        if (is_none(variant)) {
+          do_set_parameter("protocol", variant);
           return;
         }
         // check input arguments
-        m_protocol = get_value<std::shared_ptr<Protocol>>(protocol);
+        auto const protocol = get_value<std::shared_ptr<Protocol>>(variant);
         auto const shear_direction = get_shear_axis(params, "shear_direction");
         auto const shear_plane_normal =
             get_shear_axis(params, "shear_plane_normal");
         if (shear_plane_normal == shear_direction) {
           throw std::invalid_argument("Parameters 'shear_direction' and "
                                       "'shear_plane_normal' must differ");
         }
-        auto const &system = get_system();
+        auto &system = get_system();
         system.lb.lebc_sanity_checks(shear_direction, shear_plane_normal);
         // update box geometry and cell structure
-        system.box_geo->set_lees_edwards_bc(
-            LeesEdwardsBC{0., 0., shear_direction, shear_plane_normal});
-        m_lees_edwards->set_protocol(m_protocol->protocol());
+        auto const old_shear_direction = get_lebc().shear_direction;
+        auto const old_shear_plane_normal = get_lebc().shear_plane_normal;
+        try {
+          system.box_geo->set_lees_edwards_bc(
+              LeesEdwardsBC{0., 0., shear_direction, shear_plane_normal});
+          set_protocol(protocol);
+        } catch (...) {
+          system.box_geo->set_lees_edwards_bc(LeesEdwardsBC{
+              0., 0., old_shear_direction, old_shear_plane_normal});
+          set_protocol(m_protocol);
+          throw;
+        }
       });
     }
     return {};
@@ -158,6 +175,18 @@ class LeesEdwards : public AutoParameters<LeesEdwards, System::Leaf> {
     }
     m_params.reset();
   }
+
+  void set_protocol(std::shared_ptr<Protocol> const &protocol) {
+    auto &system = get_system();
+    if (protocol) {
+      m_lees_edwards->set_protocol(protocol->protocol());
+    } else {
+      system.box_geo->set_lees_edwards_bc(LeesEdwardsBC{});
+      m_lees_edwards->unset_protocol();
+    }
+    system.on_lees_edwards_change();
+    m_protocol = protocol;
+  }
 };
 
 } // namespace LeesEdwards

src/script_interface/walberla/LBFluid.cpp

@@ -158,7 +158,7 @@ void LBFluid::do_construct(VariantMap const &params) {
   auto const ext_f = get_value<Utils::Vector3d>(params, "ext_force_density");
   m_lb_params = std::make_shared<::LB::LBWalberlaParams>(agrid, tau);
   m_is_active = false;
-  m_seed = get_value<int>(params, "seed");
+  auto const seed = get_value<int>(params, "seed");
   context()->parallel_try_catch([&]() {
     if (tau <= 0.) {
       throw std::domain_error("Parameter 'tau' must be > 0");
@@ -175,7 +175,7 @@ void LBFluid::do_construct(VariantMap const &params) {
     auto const lb_dens = m_conv_dens * dens;
     auto const lb_kT = m_conv_energy * kT;
     auto const lb_ext_f = m_conv_force_dens * ext_f;
-    if (m_seed < 0) {
+    if (seed < 0) {
       throw std::domain_error("Parameter 'seed' must be >= 0");
     }
     if (lb_kT < 0.) {
@@ -188,28 +188,8 @@ void LBFluid::do_construct(VariantMap const &params) {
       throw std::domain_error("Parameter 'kinematic_viscosity' must be >= 0");
     }
     make_instance(params);
-    auto const &system = ::System::get_system();
-    if (auto le_protocol = system.lees_edwards->get_protocol()) {
-      if (lb_kT != 0.) {
-        throw std::runtime_error(
-            "Lees-Edwards LB doesn't support thermalization");
-      }
-      auto const &le_bc = system.box_geo->lees_edwards_bc();
-      auto lees_edwards_object = std::make_unique<LeesEdwardsPack>(
-          le_bc.shear_direction, le_bc.shear_plane_normal,
-          [this, le_protocol, &system]() {
-            return get_pos_offset(system.get_sim_time(), *le_protocol) /
-                   m_lb_params->get_agrid();
-          },
-          [this, le_protocol, &system]() {
-            return get_shear_velocity(system.get_sim_time(), *le_protocol) *
-                   (m_lb_params->get_tau() / m_lb_params->get_agrid());
-          });
-      m_instance->set_collision_model(std::move(lees_edwards_object));
-    } else {
-      m_instance->set_collision_model(lb_kT, m_seed);
-    }
-    m_instance->ghost_communication(); // synchronize ghost layers
+    ::LB::LBWalberla::update_collision_model(*m_instance, *m_lb_params, lb_kT,
+                                             static_cast<unsigned int>(seed));
     m_instance->set_external_force(lb_ext_f);
     for (auto &vtk : m_vtk_writers) {
       vtk->attach_to_lattice(m_instance, get_latice_to_md_units_conversion());

src/script_interface/walberla/LBFluid.hpp

@@ -51,7 +51,6 @@ class LBFluid : public LatticeModel<::LBWalberlaBase, LBVTKHandle> {
   using Base = LatticeModel<::LBWalberlaBase, LBVTKHandle>;
   std::shared_ptr<::LB::LBWalberlaParams> m_lb_params;
   bool m_is_active;
-  int m_seed;
   double m_conv_dist;
   double m_conv_visc;
   double m_conv_dens;
@@ -77,7 +76,8 @@ class LBFluid : public LatticeModel<::LBWalberlaBase, LBVTKHandle> {
          [this]() { return m_instance->get_lattice().get_grid_dimensions(); }},
         {"kT", AutoParameter::read_only,
          [this]() { return m_instance->get_kT() / m_conv_energy; }},
-        {"seed", AutoParameter::read_only, [this]() { return m_seed; }},
+        {"seed", AutoParameter::read_only,
+         [this]() { return static_cast<int>(m_instance->get_seed()); }},
         {"rng_state",
          [this](Variant const &v) {
            auto const rng_state = get_value<int>(v);