Introduce an event on_lees_edwards_change

.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);