From 220aa61c155c42bbcf37eb88df5501ea4304a5ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 21 Feb 2024 20:05:17 +0100 Subject: [PATCH] Make waLBerla dependency private --- src/core/lb/particle_coupling.cpp | 1 + src/core/unit_tests/ek_interface_test.cpp | 4 +- src/script_interface/walberla/EKFFT.hpp | 4 +- src/script_interface/walberla/EKNone.hpp | 3 +- src/script_interface/walberla/EKReaction.hpp | 26 +-- src/script_interface/walberla/EKSpecies.cpp | 3 +- src/walberla_bridge/CMakeLists.txt | 9 +- .../electrokinetics/ek_poisson_fft_init.hpp | 7 +- .../electrokinetics/ek_poisson_none_init.hpp | 7 +- .../electrokinetics/ek_walberla_init.hpp | 16 ++ .../reactions/EKReactionBase.hpp | 15 +- .../reactions/EKReactionBaseIndexed.hpp} | 32 ++- src/walberla_bridge/src/BoundaryHandling.hpp | 16 +- src/walberla_bridge/src/BoundaryPackInfo.hpp | 5 + src/walberla_bridge/src/LatticeWalberla.cpp | 3 +- .../src/electrokinetics/EKinWalberlaImpl.hpp | 13 +- .../electrokinetics/ek_poisson_fft_init.cpp | 4 + .../electrokinetics/ek_poisson_none_init.cpp | 4 + .../src/electrokinetics/ek_walberla_init.cpp | 27 ++- .../electrokinetics/reactions/CMakeLists.txt | 3 - .../reactions/EKReactionImplBulk.hpp | 23 +- .../reactions/EKReactionImplIndexed.cpp | 211 ------------------ .../reactions/EKReactionImplIndexed.hpp | 169 ++++++++++++-- .../src/lattice_boltzmann/LBWalberlaImpl.hpp | 23 +- .../src/lattice_boltzmann/ResetForce.hpp | 6 +- .../utils/boundary.hpp} | 2 +- .../utils/types_conversion.hpp} | 0 .../tests/lb_kernels_unit_tests.cpp | 3 +- 28 files changed, 312 insertions(+), 327 deletions(-) rename src/walberla_bridge/{src/electrokinetics/reactions/EKReactionImplBulk.cpp => include/walberla_bridge/electrokinetics/reactions/EKReactionBaseIndexed.hpp} (51%) delete mode 100644 src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.cpp rename src/walberla_bridge/{include/walberla_bridge/utils/boundary_utils.hpp => src/utils/boundary.hpp} (99%) rename src/walberla_bridge/{include/walberla_bridge/utils/walberla_utils.hpp => src/utils/types_conversion.hpp} (100%) diff --git a/src/core/lb/particle_coupling.cpp b/src/core/lb/particle_coupling.cpp index cf10874730..48b1afa9f1 100644 --- a/src/core/lb/particle_coupling.cpp +++ b/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 * diff --git a/src/core/unit_tests/ek_interface_test.cpp b/src/core/unit_tests/ek_interface_test.cpp index b2f5a9659f..bff055e845 100644 --- a/src/core/unit_tests/ek_interface_test.cpp +++ b/src/core/unit_tests/ek_interface_test.cpp @@ -85,7 +85,7 @@ static auto make_ek_actor() { ek_lattice = std::make_shared( params.grid_dimensions, ::communicator.node_grid, n_ghost_layers); ek_container = std::make_shared( - 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_instance = std::make_shared(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(ek_species, stoich, order); diff --git a/src/script_interface/walberla/EKFFT.hpp b/src/script_interface/walberla/EKFFT.hpp index c02fc1bae6..6ff34dcbb6 100644 --- a/src/script_interface/walberla/EKFFT.hpp +++ b/src/script_interface/walberla/EKFFT.hpp @@ -55,8 +55,8 @@ class EKFFT : public EKPoissonSolver { auto const permittivity = get_value(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", diff --git a/src/script_interface/walberla/EKNone.hpp b/src/script_interface/walberla/EKNone.hpp index 5aa1eb4e70..3c1c38cdc4 100644 --- a/src/script_interface/walberla/EKNone.hpp +++ b/src/script_interface/walberla/EKNone.hpp @@ -45,7 +45,8 @@ class EKNone : public EKPoissonSolver { m_single_precision = get_value_or(args, "single_precision", false); m_lattice = get_value>(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, diff --git a/src/script_interface/walberla/EKReaction.hpp b/src/script_interface/walberla/EKReaction.hpp index 3adc5f8847..75e2d536c1 100644 --- a/src/script_interface/walberla/EKReaction.hpp +++ b/src/script_interface/walberla/EKReaction.hpp @@ -27,9 +27,9 @@ #include "LatticeIndices.hpp" #include "LatticeWalberla.hpp" +#include #include -#include -#include +#include #include #include @@ -80,22 +80,21 @@ class EKReaction : public AutoParameters { return tau / std::pow(Utils::int_pow<3>(agrid), sum_alphas - 1.); } - template - std::shared_ptr make_instance(VariantMap const &args) const { + template + auto make_instance(VariantMap const &args, F &allocator) const { auto lattice = get_value>(args, "lattice"); - auto reactant = get_value>(args, "reactants"); - auto output = - std::vector>(reactant.size()); + auto reactants = get_value>(args, "reactants"); + auto output = ::walberla::EKReactionBase::reactants_type(reactants.size()); auto get_instance = [](Variant const &v) { return get_value>(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(args, "coefficient") * get_conversion_coefficient(); - return std::make_shared(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 diff --git a/src/script_interface/walberla/EKSpecies.cpp b/src/script_interface/walberla/EKSpecies.cpp index 9f908ad15b..c38cb1c5dc 100644 --- a/src/script_interface/walberla/EKSpecies.cpp +++ b/src/script_interface/walberla/EKSpecies.cpp @@ -24,7 +24,6 @@ #include "EKWalberlaNodeState.hpp" #include "WalberlaCheckpoint.hpp" -#include #include #include @@ -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(args, "valency"), ek_ext_efield, ek_density, get_value(args, "advection"), diff --git a/src/walberla_bridge/CMakeLists.txt b/src/walberla_bridge/CMakeLists.txt index f7912560ff..b5595397f3 100644 --- a/src/walberla_bridge/CMakeLists.txt +++ b/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) diff --git a/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_poisson_fft_init.hpp b/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_poisson_fft_init.hpp index f0a7a2db61..10337c34df 100644 --- a/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_poisson_fft_init.hpp +++ b/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_poisson_fft_init.hpp @@ -20,11 +20,14 @@ #pragma once #include - -#include "PoissonSolver/PoissonSolver.hpp" +#include #include +namespace walberla { + std::shared_ptr new_ek_poisson_fft(std::shared_ptr const &lattice, double permittivity, bool single_precision); + +} // namespace walberla diff --git a/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_poisson_none_init.hpp b/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_poisson_none_init.hpp index af7d318989..aa9890d2aa 100644 --- a/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_poisson_none_init.hpp +++ b/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_poisson_none_init.hpp @@ -20,11 +20,14 @@ #pragma once #include - -#include "PoissonSolver/PoissonSolver.hpp" +#include #include +namespace walberla { + std::shared_ptr new_ek_poisson_none(std::shared_ptr const &lattice, bool single_precision); + +} // namespace walberla diff --git a/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_walberla_init.hpp b/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_walberla_init.hpp index fb700df86c..4de3b85435 100644 --- a/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_walberla_init.hpp +++ b/src/walberla_bridge/include/walberla_bridge/electrokinetics/ek_walberla_init.hpp @@ -22,13 +22,29 @@ #include "EKinWalberlaBase.hpp" #include +#include +#include #include #include +namespace walberla { + std::shared_ptr new_ek_walberla(std::shared_ptr const &lattice, double diffusion, double kT, double valency, Utils::Vector3d ext_efield, double density, bool advection, bool friction_coupling, bool single_precision); + +std::shared_ptr +new_ek_reaction_bulk(std::shared_ptr const &lattice, + typename EKReactionBase::reactants_type const &reactants, + double coefficient); + +std::shared_ptr new_ek_reaction_indexed( + std::shared_ptr const &lattice, + typename EKReactionBase::reactants_type const &reactants, + double coefficient); + +} // namespace walberla diff --git a/src/walberla_bridge/include/walberla_bridge/electrokinetics/reactions/EKReactionBase.hpp b/src/walberla_bridge/include/walberla_bridge/electrokinetics/reactions/EKReactionBase.hpp index 049bd3226d..392c515b1a 100644 --- a/src/walberla_bridge/include/walberla_bridge/electrokinetics/reactions/EKReactionBase.hpp +++ b/src/walberla_bridge/include/walberla_bridge/electrokinetics/reactions/EKReactionBase.hpp @@ -29,18 +29,19 @@ namespace walberla { class EKReactionBase { -private: - std::vector> m_reactants; - double m_coefficient; +public: + using reactants_type = std::vector>; +private: std::shared_ptr m_lattice; + reactants_type m_reactants; + double m_coefficient; public: EKReactionBase(std::shared_ptr lattice, - std::vector> reactants, - double coefficient) - : m_reactants(std::move(reactants)), m_coefficient(coefficient), - m_lattice(std::move(lattice)) {} + reactants_type reactants, double coefficient) + : m_lattice(std::move(lattice)), m_reactants(std::move(reactants)), + m_coefficient(coefficient) {} virtual ~EKReactionBase() = default; diff --git a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplBulk.cpp b/src/walberla_bridge/include/walberla_bridge/electrokinetics/reactions/EKReactionBaseIndexed.hpp similarity index 51% rename from src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplBulk.cpp rename to src/walberla_bridge/include/walberla_bridge/electrokinetics/reactions/EKReactionBaseIndexed.hpp index 800a84af57..90d17bedb1 100644 --- a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplBulk.cpp +++ b/src/walberla_bridge/include/walberla_bridge/electrokinetics/reactions/EKReactionBaseIndexed.hpp @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022-2023 The ESPResSo project + * Copyright (C) 2024 The ESPResSo project * * This file is part of ESPResSo. * @@ -17,26 +17,24 @@ * along with this program. If not, see . */ -#include "EKReactionImplBulk.hpp" +#pragma once -#include "generated_kernels/ReactionKernelBulk_all.h" +#include "EKReactionBase.hpp" -#include +#include -namespace walberla { +#include -void EKReactionImplBulk::perform_reaction() { - // TODO: if my understanding is correct: - // the kernels need to either run in the ghost layers and do the - // synchronization before or not run and do a synchronization afterwards. - // The better solution is probably the latter one. Not sure why it fails - // atm. +namespace walberla { - auto kernel = detail::ReactionKernelBulkSelector::get_kernel( - get_reactants(), get_coefficient()); +class EKReactionBaseIndexed : public EKReactionBase { +public: + using EKReactionBase::EKReactionBase; + ~EKReactionBaseIndexed() override = default; + virtual void set_node_is_boundary(Utils::Vector3i const &node, + bool is_boundary) = 0; + virtual std::optional + get_node_is_boundary(Utils::Vector3i const &node) = 0; +}; - for (auto &block : *get_lattice()->get_blocks()) { - kernel(&block); - } -} } // namespace walberla diff --git a/src/walberla_bridge/src/BoundaryHandling.hpp b/src/walberla_bridge/src/BoundaryHandling.hpp index 86c2053888..b5d59cbe87 100644 --- a/src/walberla_bridge/src/BoundaryHandling.hpp +++ b/src/walberla_bridge/src/BoundaryHandling.hpp @@ -20,7 +20,8 @@ #pragma once #include -#include + +#include "utils/types_conversion.hpp" #include #include @@ -38,12 +39,13 @@ namespace walberla { -/// Flag for domain cells, i.e. all cells -FlagUID const Domain_flag("domain"); -/// Flag for boundary cells -FlagUID const Boundary_flag("boundary"); - template class BoundaryHandling { +private: + /** Flag for domain cells, i.e. all cells. */ + FlagUID const Domain_flag{"domain"}; + /** Flag for boundary cells. */ + FlagUID const Boundary_flag{"boundary"}; + /** Container for the map between cells and values. */ class DynamicValueCallback { public: @@ -172,7 +174,7 @@ template class BoundaryHandling { std::shared_ptr m_boundary; bool m_pending_changes; - /** Register flags and set all cells to @ref Domain_flag. */ + /** Register flags and reset all cells. */ void flag_reset_kernel(IBlock *const block) { auto flag_field = block->template getData(m_flag_field_id); // register flags diff --git a/src/walberla_bridge/src/BoundaryPackInfo.hpp b/src/walberla_bridge/src/BoundaryPackInfo.hpp index 3055b3ebe0..baeeb7c385 100644 --- a/src/walberla_bridge/src/BoundaryPackInfo.hpp +++ b/src/walberla_bridge/src/BoundaryPackInfo.hpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -38,6 +39,10 @@ template class BoundaryPackInfo : public PackInfo { protected: using PackInfo::bdId_; + /** Flag for domain cells, i.e. all cells. */ + FlagUID const Domain_flag{"domain"}; + /** Flag for boundary cells. */ + FlagUID const Boundary_flag{"boundary"}; public: using PackInfo::PackInfo; diff --git a/src/walberla_bridge/src/LatticeWalberla.cpp b/src/walberla_bridge/src/LatticeWalberla.cpp index 0559fab379..2dc2943a40 100644 --- a/src/walberla_bridge/src/LatticeWalberla.cpp +++ b/src/walberla_bridge/src/LatticeWalberla.cpp @@ -19,7 +19,8 @@ #include #include -#include + +#include "utils/types_conversion.hpp" #include #include diff --git a/src/walberla_bridge/src/electrokinetics/EKinWalberlaImpl.hpp b/src/walberla_bridge/src/electrokinetics/EKinWalberlaImpl.hpp index 9fc4019aa6..14a7ff9924 100644 --- a/src/walberla_bridge/src/electrokinetics/EKinWalberlaImpl.hpp +++ b/src/walberla_bridge/src/electrokinetics/EKinWalberlaImpl.hpp @@ -22,21 +22,21 @@ #include #include #include +#include #include #include #include #include -#include -#include +#include #include "../BoundaryHandling.hpp" +#include "../utils/boundary.hpp" +#include "../utils/types_conversion.hpp" #include "ek_kernels.hpp" #include #include #include -#include -#include #include @@ -109,6 +109,11 @@ class EKinWalberlaImpl : public EKinWalberlaBase { BlockDataID m_flag_field_density_id; BlockDataID m_flag_field_flux_id; + /** Flag for domain cells, i.e. all cells. */ + FlagUID const Domain_flag{"domain"}; + /** Flag for boundary cells. */ + FlagUID const Boundary_flag{"boundary"}; + /** Block forest */ std::shared_ptr m_lattice; diff --git a/src/walberla_bridge/src/electrokinetics/ek_poisson_fft_init.cpp b/src/walberla_bridge/src/electrokinetics/ek_poisson_fft_init.cpp index 2598c6616e..b67da67b79 100644 --- a/src/walberla_bridge/src/electrokinetics/ek_poisson_fft_init.cpp +++ b/src/walberla_bridge/src/electrokinetics/ek_poisson_fft_init.cpp @@ -23,6 +23,8 @@ #include +namespace walberla { + std::shared_ptr new_ek_poisson_fft(std::shared_ptr const &lattice, double permittivity, bool single_precision) { @@ -31,3 +33,5 @@ new_ek_poisson_fft(std::shared_ptr const &lattice, } return std::make_shared>(lattice, permittivity); } + +} // namespace walberla diff --git a/src/walberla_bridge/src/electrokinetics/ek_poisson_none_init.cpp b/src/walberla_bridge/src/electrokinetics/ek_poisson_none_init.cpp index f912cd4911..2636be3ebd 100644 --- a/src/walberla_bridge/src/electrokinetics/ek_poisson_none_init.cpp +++ b/src/walberla_bridge/src/electrokinetics/ek_poisson_none_init.cpp @@ -23,6 +23,8 @@ #include +namespace walberla { + std::shared_ptr new_ek_poisson_none(std::shared_ptr const &lattice, bool single_precision) { @@ -31,3 +33,5 @@ new_ek_poisson_none(std::shared_ptr const &lattice, } return std::make_shared>(lattice); } + +} // namespace walberla diff --git a/src/walberla_bridge/src/electrokinetics/ek_walberla_init.cpp b/src/walberla_bridge/src/electrokinetics/ek_walberla_init.cpp index a666e8f9ac..03fc1d0fae 100644 --- a/src/walberla_bridge/src/electrokinetics/ek_walberla_init.cpp +++ b/src/walberla_bridge/src/electrokinetics/ek_walberla_init.cpp @@ -18,26 +18,49 @@ */ #include "EKinWalberlaImpl.hpp" +#include "reactions/EKReactionImplBulk.hpp" +#include "reactions/EKReactionImplIndexed.hpp" #include #include +#include +#include #include #include +namespace walberla { + std::shared_ptr new_ek_walberla(std::shared_ptr const &lattice, double diffusion, double kT, double valency, Utils::Vector3d ext_efield, double density, bool advection, bool friction_coupling, bool single_precision) { if (single_precision) { - return std::make_shared>( + return std::make_shared>( lattice, diffusion, kT, valency, ext_efield, density, advection, friction_coupling); } - return std::make_shared>( + return std::make_shared>( lattice, diffusion, kT, valency, ext_efield, density, advection, friction_coupling); } + +std::shared_ptr +new_ek_reaction_bulk(std::shared_ptr const &lattice, + typename EKReactionBase::reactants_type const &reactants, + double coefficient) { + return std::make_shared(lattice, reactants, coefficient); +} + +std::shared_ptr new_ek_reaction_indexed( + std::shared_ptr const &lattice, + typename EKReactionBase::reactants_type const &reactants, + double coefficient) { + return std::make_shared(lattice, reactants, + coefficient); +} + +} // namespace walberla diff --git a/src/walberla_bridge/src/electrokinetics/reactions/CMakeLists.txt b/src/walberla_bridge/src/electrokinetics/reactions/CMakeLists.txt index 6559b5c9ff..4f5e805245 100644 --- a/src/walberla_bridge/src/electrokinetics/reactions/CMakeLists.txt +++ b/src/walberla_bridge/src/electrokinetics/reactions/CMakeLists.txt @@ -18,6 +18,3 @@ # add_subdirectory(generated_kernels) - -target_sources(espresso_walberla PRIVATE EKReactionImplBulk.cpp - EKReactionImplIndexed.cpp) diff --git a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplBulk.hpp b/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplBulk.hpp index 33f7e21770..09f38d319b 100644 --- a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplBulk.hpp +++ b/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplBulk.hpp @@ -19,28 +19,37 @@ #pragma once +#include "generated_kernels/ReactionKernelBulk_all.h" + #include #include #include -#include -#include +#include namespace walberla { class EKReactionImplBulk : public EKReactionBase { public: - EKReactionImplBulk(const std::shared_ptr &lattice, - const std::vector> &reactants, - double coefficient) - : EKReactionBase(lattice, reactants, coefficient) {} ~EKReactionImplBulk() override = default; + using EKReactionBase::EKReactionBase; using EKReactionBase::get_coefficient; using EKReactionBase::get_lattice; using EKReactionBase::get_reactants; - void perform_reaction() override; + void perform_reaction() override { + // TODO: if my understanding is correct: + // the kernels need to either run in the ghost layers and do the + // synchronization before or not run and do a synchronization afterwards. + // The better solution is probably the latter one. Not sure why it fails + // atm. + auto kernel = detail::ReactionKernelBulkSelector::get_kernel( + get_reactants(), get_coefficient()); + for (auto &block : *get_lattice()->get_blocks()) { + kernel(&block); + } + } }; } // namespace walberla diff --git a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.cpp b/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.cpp deleted file mode 100644 index c3f5643fd9..0000000000 --- a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.cpp +++ /dev/null @@ -1,211 +0,0 @@ -/* - * Copyright (C) 2022-2023 The ESPResSo project - * - * This file is part of ESPResSo. - * - * ESPResSo is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * ESPResSo is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#include "EKReactionImplIndexed.hpp" - -#include "generated_kernels/ReactionKernelIndexed_all.h" - -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -namespace walberla { - -/// Flag for domain cells, i.e. all cells -FlagUID const Domain_flag("domain"); -/// Flag for boundary cells -FlagUID const Boundary_flag("boundary"); - -namespace detail { -// FlagField to use -using FlagField = FlagField; - -template -inline auto -get_flag_field_and_flag(IBlock *block, - domain_decomposition::BlockDataID const &flagfield_id) { - auto const flag_field = - block->template uncheckedFastGetData(flagfield_id); - auto const boundary_flag = flag_field->getFlag(Boundary_flag); - return std::make_tuple(flag_field, boundary_flag); -} - -template -void fillFromFlagField(IBlock *block, BlockDataID indexVectorID, - ConstBlockDataID flagFieldID, FlagUID boundaryFlagUID, - FlagUID domainFlagUID) { - auto *indexVectors = block->uncheckedFastGetData(indexVectorID); - auto &indexVectorAll = indexVectors->indexVector(IndexVectors::ALL); - auto &indexVectorInner = indexVectors->indexVector(IndexVectors::INNER); - auto &indexVectorOuter = indexVectors->indexVector(IndexVectors::OUTER); - - auto *flagField = block->getData(flagFieldID); - - assert(flagField->flagExists(boundaryFlagUID) and - flagField->flagExists(domainFlagUID)); - - auto boundaryFlag = flagField->getFlag(boundaryFlagUID); - auto domainFlag = flagField->getFlag(domainFlagUID); - - auto inner = flagField->xyzSize(); - inner.expand(cell_idx_t(-1)); - - indexVectorAll.clear(); - indexVectorInner.clear(); - indexVectorOuter.clear(); - - auto flagWithGLayers = flagField->xyzSizeWithGhostLayer(); - for (auto it = flagField->beginWithGhostLayerXYZ(); it != flagField->end(); - ++it) { - - if (!isFlagSet(it, boundaryFlag)) - continue; - if (flagWithGLayers.contains(it.x(), it.y(), it.z()) && - isFlagSet(it.neighbor(0, 0, 0, 0), domainFlag)) { - - auto element = IndexInfo(it.x(), it.y(), it.z()); - - indexVectorAll.push_back(element); - if (inner.contains(it.x(), it.y(), it.z())) - indexVectorInner.push_back(element); - else - indexVectorOuter.push_back(element); - } - } - - indexVectors->syncGPU(); -} - -template -void fillFromFlagField(const std::shared_ptr &blocks, - BlockDataID indexVectorID, ConstBlockDataID flagFieldID, - FlagUID boundaryFlagUID, FlagUID domainFlagUID) { - for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) - fillFromFlagField( - blockIt.get(), indexVectorID, flagFieldID, boundaryFlagUID, - domainFlagUID); -} -} // namespace detail - -EKReactionImplIndexed::EKReactionImplIndexed( - std::shared_ptr lattice, - std::vector> reactants, double coefficient) - : EKReactionBase(lattice, reactants, coefficient), - m_pending_changes(false) { - m_flagfield_id = - static_cast(field::addFlagFieldToStorage( - get_lattice()->get_blocks(), "flag field reaction", - get_lattice()->get_ghost_layers())); - - // take one IndexVector as a dummy-value - using IndexVectors = detail::ReactionKernelIndexedSelector::KernelTrait<>:: - ReactionKernelIndexed::IndexVectors; - - auto createIdxVector = [](IBlock *const, StructuredBlockStorage *const) { - return new IndexVectors(); - }; - m_indexvector_id = static_cast( - get_lattice() - ->get_blocks() - ->template addStructuredBlockData(createIdxVector, - "IndexField")); - - for (auto &block : *get_lattice()->get_blocks()) { - auto flag_field = - block.template getData(BlockDataID(m_flagfield_id)); - // register flags - flag_field->registerFlag(Domain_flag); - flag_field->registerFlag(Boundary_flag); - // mark all cells as domain cells and fluid cells - auto domain_flag = flag_field->getFlag(Domain_flag); - auto boundary_flag = flag_field->getFlag(Boundary_flag); - for (auto it = flag_field->begin(); it != flag_field->end(); ++it) { - flag_field->addFlag(it.x(), it.y(), it.z(), domain_flag); - flag_field->removeFlag(it.x(), it.y(), it.z(), boundary_flag); - } - } -} - -void EKReactionImplIndexed::perform_reaction() { - boundary_update(); - - auto kernel = detail::ReactionKernelIndexedSelector::get_kernel( - get_reactants(), get_coefficient(), BlockDataID(get_indexvector_id())); - - for (auto &block : *get_lattice()->get_blocks()) { - kernel(&block); - } -} - -void EKReactionImplIndexed::set_node_is_boundary(Utils::Vector3i const &node, - bool is_boundary) { - auto bc = get_block_and_cell(*get_lattice(), node, true); - if (!bc) - return; - - auto [flag_field, boundary_flag] = - detail::get_flag_field_and_flag( - bc->block, BlockDataID(get_flagfield_id())); - if (is_boundary) { - flag_field->addFlag(bc->cell, boundary_flag); - } else { - flag_field->removeFlag(bc->cell, boundary_flag); - } - m_pending_changes = true; -} - -std::optional -EKReactionImplIndexed::get_node_is_boundary(Utils::Vector3i const &node) { - auto bc = get_block_and_cell(*get_lattice(), node, true); - if (!bc) - return std::nullopt; - - auto [flag_field, boundary_flag] = - detail::get_flag_field_and_flag( - bc->block, BlockDataID(get_flagfield_id())); - return {flag_field->isFlagSet(bc->cell, boundary_flag)}; -} - -void EKReactionImplIndexed::boundary_update() { - // take one IndexVector/IndexInfo as a dummy-value - using IndexVectors = detail::ReactionKernelIndexedSelector::KernelTrait<>:: - ReactionKernelIndexed::IndexVectors; - using IndexInfo = detail::ReactionKernelIndexedSelector::KernelTrait<>:: - ReactionKernelIndexed::IndexInfo; - - if (m_pending_changes) { - detail::fillFromFlagField( - get_lattice()->get_blocks(), BlockDataID(get_indexvector_id()), - BlockDataID(get_flagfield_id()), Boundary_flag, Domain_flag); - m_pending_changes = false; - } -} -} // namespace walberla diff --git a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.hpp b/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.hpp index 65686644f7..48125caf36 100644 --- a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.hpp +++ b/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.hpp @@ -19,50 +19,179 @@ #pragma once +#include "generated_kernels/ReactionKernelIndexed_all.h" + +#include #include #include -#include +#include + +#include +#include +#include +#include +#include +#include #include +#include #include #include #include +#include #include namespace walberla { -class EKReactionImplIndexed : public EKReactionBase { +class EKReactionImplIndexed : public EKReactionBaseIndexed { private: - std::size_t m_flagfield_id; - std::size_t m_indexvector_id; - + BlockDataID m_flagfield_id; + BlockDataID m_indexvector_id; bool m_pending_changes; public: - EKReactionImplIndexed(std::shared_ptr lattice, - std::vector> reactants, - double coefficient); + /** Flag for domain cells, i.e. all cells. */ + FlagUID const Domain_flag{"domain"}; + /** Flag for boundary cells. */ + FlagUID const Boundary_flag{"boundary"}; + + using FlagField = field::FlagField; + using IndexVectors = detail::ReactionKernelIndexedSelector::KernelTrait<>:: + ReactionKernelIndexed::IndexVectors; + using IndexInfo = detail::ReactionKernelIndexedSelector::KernelTrait<>:: + ReactionKernelIndexed::IndexInfo; + +private: + auto get_flag_field_and_flag(IBlock *block, BlockDataID const &flagfield_id) { + auto const flag_field = + block->template uncheckedFastGetData(flagfield_id); + auto const boundary_flag = flag_field->getFlag(Boundary_flag); + return std::make_tuple(flag_field, boundary_flag); + } + +public: + EKReactionImplIndexed(std::shared_ptr const &lattice, + reactants_type const &reactants, double coefficient) + : EKReactionBaseIndexed(lattice, reactants, coefficient), + m_pending_changes(false) { + m_flagfield_id = field::addFlagFieldToStorage( + get_lattice()->get_blocks(), "flag field reaction", + get_lattice()->get_ghost_layers()); + + auto createIdxVector = [](IBlock *const, StructuredBlockStorage *const) { + return new IndexVectors(); + }; + m_indexvector_id = get_lattice() + ->get_blocks() + ->template addStructuredBlockData( + createIdxVector, "IndexField"); + + for (auto &block : *get_lattice()->get_blocks()) { + auto flag_field = block.template getData(m_flagfield_id); + // register flags + flag_field->registerFlag(Domain_flag); + flag_field->registerFlag(Boundary_flag); + // mark all cells as domain cells and fluid cells + auto domain_flag = flag_field->getFlag(Domain_flag); + auto boundary_flag = flag_field->getFlag(Boundary_flag); + for (auto it = flag_field->begin(); it != flag_field->end(); ++it) { + flag_field->addFlag(it.x(), it.y(), it.z(), domain_flag); + flag_field->removeFlag(it.x(), it.y(), it.z(), boundary_flag); + } + } + } ~EKReactionImplIndexed() override = default; - using EKReactionBase::get_coefficient; - using EKReactionBase::get_lattice; - using EKReactionBase::get_reactants; + using EKReactionBaseIndexed::get_coefficient; + using EKReactionBaseIndexed::get_lattice; + using EKReactionBaseIndexed::get_reactants; - void perform_reaction() override; + void perform_reaction() override { + boundary_update(); + auto kernel = detail::ReactionKernelIndexedSelector::get_kernel( + get_reactants(), get_coefficient(), m_indexvector_id); + for (auto &block : *get_lattice()->get_blocks()) { + kernel(&block); + } + } - void set_node_is_boundary(Utils::Vector3i const &node, bool is_boundary); - [[nodiscard]] std::optional - get_node_is_boundary(Utils::Vector3i const &node); + void set_node_is_boundary(Utils::Vector3i const &node, + bool is_boundary) override { + if (auto bc = get_block_and_cell(*get_lattice(), node, true)) { + auto const [flag_field, boundary_flag] = + get_flag_field_and_flag(bc->block, m_flagfield_id); + if (is_boundary) { + flag_field->addFlag(bc->cell, boundary_flag); + } else { + flag_field->removeFlag(bc->cell, boundary_flag); + } + m_pending_changes = true; + } + } - [[nodiscard]] auto get_indexvector_id() const noexcept { - return m_indexvector_id; + [[nodiscard]] std::optional + get_node_is_boundary(Utils::Vector3i const &node) override { + if (auto bc = get_block_and_cell(*get_lattice(), node, true)) { + auto const [flag_field, boundary_flag] = + get_flag_field_and_flag(bc->block, m_flagfield_id); + return {flag_field->isFlagSet(bc->cell, boundary_flag)}; + } + return std::nullopt; } - [[nodiscard]] auto get_flagfield_id() const noexcept { - return m_flagfield_id; + + void boundary_update() { + if (m_pending_changes) { + for (auto &block : *get_lattice()->get_blocks()) { + fillFromFlagField(block); + } + m_pending_changes = false; + } } - void boundary_update(); +private: + void fillFromFlagField(IBlock &block) { + auto *indexVectors = + block.uncheckedFastGetData(m_indexvector_id); + auto &indexVectorAll = indexVectors->indexVector(IndexVectors::ALL); + auto &indexVectorInner = indexVectors->indexVector(IndexVectors::INNER); + auto &indexVectorOuter = indexVectors->indexVector(IndexVectors::OUTER); + + auto *flagField = block.getData(m_flagfield_id); + + assert(flagField->flagExists(Boundary_flag) and + flagField->flagExists(Domain_flag)); + + auto boundaryFlag = flagField->getFlag(Boundary_flag); + auto domainFlag = flagField->getFlag(Domain_flag); + + auto inner = flagField->xyzSize(); + inner.expand(cell_idx_t(-1)); + + indexVectorAll.clear(); + indexVectorInner.clear(); + indexVectorOuter.clear(); + + auto flagWithGLayers = flagField->xyzSizeWithGhostLayer(); + for (auto it = flagField->beginWithGhostLayerXYZ(); it != flagField->end(); + ++it) { + if (!isFlagSet(it, boundaryFlag)) + continue; + + if (flagWithGLayers.contains(it.x(), it.y(), it.z()) && + isFlagSet(it.neighbor(0, 0, 0, 0), domainFlag)) { + auto element = IndexInfo(it.x(), it.y(), it.z()); + indexVectorAll.push_back(element); + if (inner.contains(it.x(), it.y(), it.z())) { + indexVectorInner.push_back(element); + } else { + indexVectorOuter.push_back(element); + } + } + } + + indexVectors->syncGPU(); + } }; } // namespace walberla diff --git a/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp b/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp index ecf7e57a64..5270c26aa4 100644 --- a/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp +++ b/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp @@ -29,23 +29,20 @@ #include #include #include +#include #include +#include #include +#include #include #include - -#include -#include -#include -#include -#include -#include - #include #include #include "../BoundaryHandling.hpp" #include "../BoundaryPackInfo.hpp" +#include "../utils/boundary.hpp" +#include "../utils/types_conversion.hpp" #include "InterpolateAndShiftAtBoundary.hpp" #include "ResetForce.hpp" #include "lb_kernels.hpp" @@ -55,8 +52,6 @@ #include #include #include -#include -#include #include #include @@ -71,6 +66,7 @@ #include #include #include +#include #include #include #include @@ -106,8 +102,8 @@ class LBWalberlaImpl : public LBWalberlaBase { protected: template struct FieldTrait { - using PdfField = GhostLayerField; - using VectorField = GhostLayerField; + using PdfField = field::GhostLayerField; + using VectorField = field::GhostLayerField; template using PackInfo = field::communication::PackInfo; }; @@ -217,6 +213,9 @@ class LBWalberlaImpl : public LBWalberlaBase { BlockDataID m_velocity_field_id; BlockDataID m_vec_tmp_field_id; + /** Flag for boundary cells. */ + FlagUID const Boundary_flag{"boundary"}; + /** * @brief Full communicator. * We use the D3Q27 directions to update cells along the diagonals during diff --git a/src/walberla_bridge/src/lattice_boltzmann/ResetForce.hpp b/src/walberla_bridge/src/lattice_boltzmann/ResetForce.hpp index 1d0a154e5b..ce7d19295a 100644 --- a/src/walberla_bridge/src/lattice_boltzmann/ResetForce.hpp +++ b/src/walberla_bridge/src/lattice_boltzmann/ResetForce.hpp @@ -22,11 +22,11 @@ #include "generated_kernels/FieldAccessorsDoublePrecision.h" #include "generated_kernels/FieldAccessorsSinglePrecision.h" -#include +#include "../utils/types_conversion.hpp" #include -#include -#include +#include +#include #include diff --git a/src/walberla_bridge/include/walberla_bridge/utils/boundary_utils.hpp b/src/walberla_bridge/src/utils/boundary.hpp similarity index 99% rename from src/walberla_bridge/include/walberla_bridge/utils/boundary_utils.hpp rename to src/walberla_bridge/src/utils/boundary.hpp index 8c1558ee72..7d3f3cdb07 100644 --- a/src/walberla_bridge/include/walberla_bridge/utils/boundary_utils.hpp +++ b/src/walberla_bridge/src/utils/boundary.hpp @@ -19,7 +19,7 @@ #pragma once -#include "walberla_utils.hpp" +#include "types_conversion.hpp" #include diff --git a/src/walberla_bridge/include/walberla_bridge/utils/walberla_utils.hpp b/src/walberla_bridge/src/utils/types_conversion.hpp similarity index 100% rename from src/walberla_bridge/include/walberla_bridge/utils/walberla_utils.hpp rename to src/walberla_bridge/src/utils/types_conversion.hpp diff --git a/src/walberla_bridge/tests/lb_kernels_unit_tests.cpp b/src/walberla_bridge/tests/lb_kernels_unit_tests.cpp index 2666843c37..3045e9974c 100644 --- a/src/walberla_bridge/tests/lb_kernels_unit_tests.cpp +++ b/src/walberla_bridge/tests/lb_kernels_unit_tests.cpp @@ -29,8 +29,7 @@ #include "../src/lattice_boltzmann/generated_kernels/Dynamic_UBB_single_precision.h" #include "../src/lattice_boltzmann/generated_kernels/FieldAccessorsDoublePrecision.h" #include "../src/lattice_boltzmann/generated_kernels/FieldAccessorsSinglePrecision.h" - -#include +#include "../src/utils/types_conversion.hpp" #include