Refactor waLBerla bridge (#4864)

Fixes #4859, fixes #4855

Description of changes:
- bugfix:
   - LB boundaries are now properly communicated in the ghost layer
   - see details in #4859
- performance:
   - the LB flag field is no longer communicated at every time step
   - the LB UBB field is no longer recalculated at every time step
   - LB boundary setters (node, slice, shape) now always trigger a full ghost communication
- maintainability:
   - the waLBerla header files are no longer visible in the ESPResSo core and script interface
   - the Boost dependency is now checked at the CMake level to prevent building broken ESPResSo shared libraries (see details in #4856)

.github/workflows/push_pull.yml

@@ -10,7 +10,7 @@ permissions:
 jobs:
   macos:
     runs-on: macos-12
-    if: ${{ github.repository == 'espressomd/espresso' }}
+    if: false
     steps:
       - name: Checkout
         uses: actions/checkout@main

CMakeLists.txt

@@ -435,6 +435,9 @@ if(ESPRESSO_BUILD_TESTS)
 endif()
 
 find_package(Boost 1.74.0 REQUIRED ${BOOST_COMPONENTS})
+if(${Boost_VERSION} VERSION_GREATER_EQUAL 1.84.0)
+  message(FATAL_ERROR "Boost version ${Boost_VERSION} is unsupported.")
+endif()
 
 #
 # Paths

maintainer/benchmarks/lb.py

@@ -53,11 +53,11 @@
 
 # process and check arguments
 n_iterations = 30
-assert args.volume_fraction > 0, "volume_fraction must be a positive number"
+assert args.volume_fraction > 0, "--volume_fraction must be a positive number"
 assert args.volume_fraction < np.pi / (3 * np.sqrt(2)), \
-    "volume_fraction exceeds the physical limit of sphere packing (~0.74)"
+    "--volume_fraction exceeds the physical limit of sphere packing (~0.74)"
 assert "box_l" not in args or args.particles_per_core == 0, \
-    "Argument box_l requires particles_per_core=0"
+    "Argument --box_l requires --particles_per_core=0"
 
 required_features = ["LENNARD_JONES", "WALBERLA"]
 if args.gpu:
@@ -85,6 +85,7 @@
 if n_part == 0:
     box_l = args.box_l
     agrid = 1.
+    lb_grid = args.box_l
     measurement_steps = 80
 else:
     # volume of N spheres with radius r: N * (4/3*pi*r^3)

maintainer/walberla_kernels/generate_ek_kernels.py

@@ -186,12 +186,6 @@ def replace_getData_with_uncheckedFastGetData(filename: str) -> None:
         index_shape=density_field.index_shape,
         target=target)
 
-    pystencils_walberla.generate_pack_info_from_kernel(
-        ctx,
-        f"DensityPackInfo_{precision_suffix}",
-        ek_electrostatic.continuity(),
-        target=target)
-
     # ek reactions
     for i in range(1, max_num_reactants + 1):
         assignments = list(reaction_obj.generate_reaction(num_reactants=i))

maintainer/walberla_kernels/templates/Boundary.tmpl.h

@@ -40,8 +40,9 @@
 #include <field/FlagField.h>
 #include <core/debug/Debug.h>
 
+#include <cassert>
 #include <functional>
-#include <set>
+#include <memory>
 #include <vector>
 
 {% for header in interface_spec.headers %}
@@ -122,7 +123,7 @@ class {{class_name}}
         {%- endif %}
     };
 
-    {{class_name}}( const shared_ptr<StructuredBlockForest> & blocks,
+    {{class_name}}( const std::shared_ptr<StructuredBlockForest> & blocks,
                    {{kernel|generate_constructor_parameters(['indexVector', 'indexVectorSize'])}}{{additional_data_handler.constructor_arguments}})
         :{{additional_data_handler.initialiser_list}} {{ kernel|generate_constructor_initializer_list(['indexVector', 'indexVectorSize']) }}
     {
@@ -177,7 +178,7 @@ class {{class_name}}
     }
 
     template<typename FlagField_T>
-    void fillFromFlagField( const shared_ptr<StructuredBlockForest> & blocks, ConstBlockDataID flagFieldID,
+    void fillFromFlagField( const std::shared_ptr<StructuredBlockForest> & blocks, ConstBlockDataID flagFieldID,
                             FlagUID boundaryFlagUID, FlagUID domainFlagUID)
     {
         for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
@@ -197,8 +198,8 @@ class {{class_name}}
         auto * flagField = block->getData< FlagField_T > ( flagFieldID );
         {{additional_data_handler.additional_field_data|indent(4)}}
 
-        if( !(flagField->flagExists(boundaryFlagUID) && flagField->flagExists(domainFlagUID) ))
-            return;
+        assert(flagField->flagExists(boundaryFlagUID and
+               flagField->flagExists(domainFlagUID));
 
         auto boundaryFlag = flagField->getFlag(boundaryFlagUID);
         auto domainFlag = flagField->getFlag(domainFlagUID);

samples/lb_circular_couette.py

@@ -57,7 +57,7 @@
 cyl_center = agrid * (grid_size // 2 + 0.5) * [1, 1, 0]
 cylinder_in = espressomd.shapes.Cylinder(
     center=cyl_center, axis=[0, 0, 1], length=3 * system.box_l[2],
-    radius=8.1 * agrid, direction=1)
+    radius=8.6 * agrid, direction=1)
 cylinder_out = espressomd.shapes.Cylinder(
     center=cyl_center, axis=[0, 0, 1], length=3 * system.box_l[2],
     radius=14.5 * agrid, direction=-1)

src/core/electrostatics/p3m.cpp

@@ -599,9 +599,11 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
       }
     }
     energy *= prefactor;
+#ifdef NPT
     if (npt_flag) {
       npt_add_virial_contribution(energy);
     }
+#endif
     if (not energy_flag) {
       energy = 0.;
     }

src/core/lb/particle_coupling.cpp

@@ -86,6 +86,7 @@ Utils::Vector3d lb_drag_force(LB::Solver const &lb, double lb_gamma,
 /**
  * @brief Check if a position is within the local box + halo.
  *
+ * @param local_box Local geometry
  * @param pos Position to check
  * @param halo Halo
  *

src/core/unit_tests/ek_interface_test.cpp

@@ -85,7 +85,7 @@ static auto make_ek_actor() {
   ek_lattice = std::make_shared<LatticeWalberla>(
       params.grid_dimensions, ::communicator.node_grid, n_ghost_layers);
   ek_container = std::make_shared<EK::EKWalberla::ek_container_type>(
-      params.tau, new_ek_poisson_none(ek_lattice, single_precision));
+      params.tau, walberla::new_ek_poisson_none(ek_lattice, single_precision));
   ek_reactions = std::make_shared<EK::EKWalberla::ek_reactions_type>();
   ek_instance = std::make_shared<EK::EKWalberla>(ek_container, ek_reactions);
 #endif
@@ -146,7 +146,7 @@ BOOST_AUTO_TEST_CASE(ek_interface_walberla) {
     auto constexpr single_precision = true;
     auto constexpr stoich = 1.;
     auto constexpr order = 2.;
-    auto ek_species = new_ek_walberla(
+    auto ek_species = walberla::new_ek_walberla(
         espresso::ek_lattice, params.diffusion, params.kT, params.valency,
         params.ext_efield, params.density, false, false, single_precision);
     auto ek_reactant = std::make_shared<EKReactant>(ek_species, stoich, order);

src/python/espressomd/lb.py

@@ -53,20 +53,6 @@ def __getitem__(self, key):
     def __str__(self):
         return f"{self.__class__.__name__}({self.get_params()})"
 
-    def _activate(self):
-        self._activate_method()
-
-    def _deactivate(self):
-        self._deactivate_method()
-
-    def _activate_method(self):
-        self.call_method("activate")
-        utils.handle_errors("HydrodynamicInteraction activation failed")
-
-    def _deactivate_method(self):
-        self.call_method("deactivate")
-        utils.handle_errors("HydrodynamicInteraction deactivation failed")
-
     def validate_params(self, params):
         pass
 
@@ -342,13 +328,6 @@ class LBFluidNodeWalberla(ScriptInterfaceHelper):
     def required_keys(self):
         return {"parent_sip", "index"}
 
-    def __init__(self, *args, **kwargs):
-        if "sip" not in kwargs:
-            super().__init__(*args, **kwargs)
-            utils.handle_errors("LBFluidNode instantiation failed")
-        else:
-            super().__init__(**kwargs)
-
     def __reduce__(self):
         raise NotImplementedError("Cannot serialize LB fluid node objects")
 
@@ -494,7 +473,6 @@ def __init__(self, *args, **kwargs):
                 slice_range, grid_size)
             node = LBFluidNodeWalberla(index=np.array([0, 0, 0]), **kwargs)
             super().__init__(*args, node_sip=node, **kwargs, **extra_kwargs)
-            utils.handle_errors("LBFluidSliceWalberla instantiation failed")
 
     def __iter__(self):
         lower, upper = self.call_method("get_slice_ranges")

src/script_interface/electrostatics/CoulombScafacos.hpp

@@ -31,6 +31,7 @@
 #include "script_interface/get_value.hpp"
 #include "script_interface/scafacos/scafacos.hpp"
 
+#include <iomanip>
 #include <regex>
 #include <set>
 #include <sstream>

src/script_interface/walberla/EKFFT.hpp

@@ -55,8 +55,8 @@ class EKFFT : public EKPoissonSolver {
     auto const permittivity =
         get_value<double>(args, "permittivity") * m_conv_permittivity;
 
-    m_instance = new_ek_poisson_fft(m_lattice->lattice(), permittivity,
-                                    m_single_precision);
+    m_instance = ::walberla::new_ek_poisson_fft(
+        m_lattice->lattice(), permittivity, m_single_precision);
 
     add_parameters({
         {"permittivity",

src/script_interface/walberla/EKNone.hpp

@@ -45,7 +45,8 @@ class EKNone : public EKPoissonSolver {
     m_single_precision = get_value_or<bool>(args, "single_precision", false);
     m_lattice = get_value<std::shared_ptr<LatticeWalberla>>(args, "lattice");
 
-    m_instance = new_ek_poisson_none(m_lattice->lattice(), m_single_precision);
+    m_instance = ::walberla::new_ek_poisson_none(m_lattice->lattice(),
+                                                 m_single_precision);
 
     add_parameters({
         {"single_precision", AutoParameter::read_only,

src/script_interface/walberla/EKReaction.hpp

@@ -27,9 +27,9 @@
 #include "LatticeIndices.hpp"
 #include "LatticeWalberla.hpp"
 
+#include <walberla_bridge/electrokinetics/ek_walberla_init.hpp>
 #include <walberla_bridge/electrokinetics/reactions/EKReactionBase.hpp>
-#include <walberla_bridge/src/electrokinetics/reactions/EKReactionImplBulk.hpp>
-#include <walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.hpp>
+#include <walberla_bridge/electrokinetics/reactions/EKReactionBaseIndexed.hpp>
 
 #include <script_interface/ScriptInterface.hpp>
 #include <script_interface/auto_parameters/AutoParameters.hpp>
@@ -80,22 +80,21 @@ class EKReaction : public AutoParameters<EKReaction, LatticeIndices> {
     return tau / std::pow(Utils::int_pow<3>(agrid), sum_alphas - 1.);
   }
 
-  template <typename T>
-  std::shared_ptr<T> make_instance(VariantMap const &args) const {
+  template <typename F>
+  auto make_instance(VariantMap const &args, F &allocator) const {
     auto lattice = get_value<std::shared_ptr<LatticeWalberla>>(args, "lattice");
-    auto reactant = get_value<std::vector<Variant>>(args, "reactants");
-    auto output =
-        std::vector<std::shared_ptr<::walberla::EKReactant>>(reactant.size());
+    auto reactants = get_value<std::vector<Variant>>(args, "reactants");
+    auto output = ::walberla::EKReactionBase::reactants_type(reactants.size());
     auto get_instance = [](Variant const &v) {
       return get_value<std::shared_ptr<EKReactant>>(v)->get_instance();
     };
-    std::transform(reactant.begin(), reactant.end(), output.begin(),
+    std::transform(reactants.begin(), reactants.end(), output.begin(),
                    get_instance);
 
     auto const coefficient =
         get_value<double>(args, "coefficient") * get_conversion_coefficient();
 
-    return std::make_shared<T>(lattice->lattice(), output, coefficient);
+    return allocator(lattice->lattice(), output, coefficient);
   }
 
   std::shared_ptr<::walberla::EKReactionBase> m_ekreaction;
@@ -118,7 +117,7 @@ class EKBulkReaction : public EKReaction {
 
   void do_construct(VariantMap const &args) override {
     m_conv_coefficient = calculate_bulk_conversion_factor(args);
-    m_ekreaction = make_instance<::walberla::EKReactionImplBulk>(args);
+    m_ekreaction = make_instance(args, ::walberla::new_ek_reaction_bulk);
   }
 };
 
@@ -143,10 +142,9 @@ class EKIndexedReaction : public EKReaction {
   void do_construct(VariantMap const &args) override {
     auto const agrid = get_agrid(args);
     m_conv_coefficient = calculate_bulk_conversion_factor(args) / agrid;
-    m_ekreaction = make_instance<::walberla::EKReactionImplIndexed>(args);
     m_ekreaction_impl =
-        std::dynamic_pointer_cast<::walberla::EKReactionImplIndexed>(
-            get_instance());
+        make_instance(args, ::walberla::new_ek_reaction_indexed);
+    m_ekreaction = m_ekreaction_impl;
   }
 
   [[nodiscard]] Variant do_call_method(std::string const &method,
@@ -170,7 +168,7 @@ class EKIndexedReaction : public EKReaction {
   }
 
 private:
-  std::shared_ptr<::walberla::EKReactionImplIndexed> m_ekreaction_impl;
+  std::shared_ptr<::walberla::EKReactionBaseIndexed> m_ekreaction_impl;
 };
 
 } // namespace ScriptInterface::walberla

src/script_interface/walberla/EKSpecies.cpp

@@ -24,7 +24,6 @@
 #include "EKWalberlaNodeState.hpp"
 #include "WalberlaCheckpoint.hpp"
 
-#include <walberla_bridge/LatticeWalberla.hpp>
 #include <walberla_bridge/electrokinetics/ek_walberla_init.hpp>
 
 #include <boost/mpi.hpp>
@@ -119,7 +118,7 @@ void EKSpecies::do_construct(VariantMap const &args) {
     auto const ek_ext_efield = ext_efield * m_conv_ext_efield;
     auto const ek_density = m_density = density * m_conv_density;
     auto const ek_kT = kT * m_conv_energy;
-    m_instance = new_ek_walberla(
+    m_instance = ::walberla::new_ek_walberla(
         m_lattice->lattice(), ek_diffusion, ek_kT,
         get_value<double>(args, "valency"), ek_ext_efield, ek_density,
         get_value<bool>(args, "advection"),

src/script_interface/walberla/LBFluid.cpp

@@ -108,7 +108,6 @@ Variant LBFluid::do_call_method(std::string const &name,
   }
   if (name == "clear_boundaries") {
     m_instance->clear_boundaries();
-    m_instance->ghost_communication();
     ::System::get_system().on_lb_boundary_conditions_change();
     return {};
   }
@@ -269,8 +268,8 @@ void LBFluid::load_checkpoint(std::string const &filename, int mode) {
   };
 
   auto const on_success = [&lb_obj]() {
-    lb_obj.reallocate_ubb_field();
     lb_obj.ghost_communication();
+    lb_obj.reallocate_ubb_field();
   };
 
   load_checkpoint_common(*context(), "LB", filename, mode, read_metadata,

src/script_interface/walberla/LBFluidNode.cpp

@@ -47,13 +47,13 @@ Variant LBFluidNode::do_call_method(std::string const &name,
   if (name == "set_velocity_at_boundary") {
     if (is_none(params.at("value"))) {
       m_lb_fluid->remove_node_from_boundary(m_index);
-      m_lb_fluid->ghost_communication();
     } else {
       auto const u =
           get_value<Utils::Vector3d>(params, "value") * m_conv_velocity;
       m_lb_fluid->set_node_velocity_at_boundary(m_index, u);
-      m_lb_fluid->ghost_communication();
     }
+    m_lb_fluid->ghost_communication();
+    m_lb_fluid->reallocate_ubb_field();
     return {};
   }
   if (name == "get_velocity_at_boundary") {

src/script_interface/walberla/LBFluidSlice.cpp

@@ -99,8 +99,10 @@ Variant LBFluidSlice::do_call_method(std::string const &name,
                 1. / m_conv_velocity);
   }
   if (name == "set_velocity_at_boundary") {
-    return call(&LatticeModel::set_slice_velocity_at_boundary, {1},
-                m_conv_velocity);
+    auto const retval = call(&LatticeModel::set_slice_velocity_at_boundary, {1},
+                             m_conv_velocity);
+    m_lb_fluid->reallocate_ubb_field();
+    return retval;
   }
   if (name == "get_pressure_tensor") {
     return call(&LatticeModel::get_slice_pressure_tensor, {3, 3},

src/walberla_bridge/CMakeLists.txt

@@ -42,8 +42,8 @@ if(ESPRESSO_BUILD_WITH_CUDA AND WALBERLA_BUILD_WITH_CUDA)
                         PRIVATE ${WALBERLA_LIBS})
   target_include_directories(espresso_walberla_cuda PUBLIC include)
   target_include_directories(
-    espresso_walberla_cuda PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
-    PRIVATE ${WALBERLA_INCLUDE_DIRS} ${walberla_BINARY_DIR}/src)
+    espresso_walberla_cuda PRIVATE ${WALBERLA_INCLUDE_DIRS}
+                                   ${walberla_BINARY_DIR}/src)
   install(TARGETS espresso_walberla_cuda
           LIBRARY DESTINATION ${ESPRESSO_INSTALL_PYTHON}/espressomd)
   target_link_libraries(espresso_walberla PUBLIC espresso::walberla_cuda)
@@ -52,9 +52,8 @@ endif()
 target_link_libraries(
   espresso_walberla PUBLIC MPI::MPI_CXX espresso::utils
   PRIVATE espresso::cpp_flags espresso::walberla::cpp_flags ${WALBERLA_LIBS})
-target_include_directories(
-  espresso_walberla PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
-  PRIVATE ${WALBERLA_INCLUDE_DIRS} ${walberla_BINARY_DIR}/src)
+target_include_directories(espresso_walberla PRIVATE ${WALBERLA_INCLUDE_DIRS}
+                                                     ${walberla_BINARY_DIR}/src)
 
 add_subdirectory(src)
 

src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_poisson_fft_init.hpp

@@ -20,11 +20,14 @@
 #pragma once
 
 #include <walberla_bridge/LatticeWalberla.hpp>
-
-#include "PoissonSolver/PoissonSolver.hpp"
+#include <walberla_bridge/electrokinetics/PoissonSolver/PoissonSolver.hpp>
 
 #include <memory>
 
+namespace walberla {
+
 std::shared_ptr<walberla::PoissonSolver>
 new_ek_poisson_fft(std::shared_ptr<LatticeWalberla> const &lattice,
                    double permittivity, bool single_precision);
+
+} // namespace walberla