From f3f24616d679bc0ff0f19955fa4a0b56fc2bb843 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Mon, 19 Aug 2019 18:28:00 +0200 Subject: [PATCH] WIP: Walberla boundary force --- src/core/CMakeLists.txt | 1 + src/core/grid_based_algorithms/LbWalberla.cpp | 22 +++++++++- src/core/grid_based_algorithms/LbWalberla.hpp | 5 ++- .../lbboundaries/LBBoundary.cpp | 41 +++++++++++++++++++ .../lbboundaries/LBBoundary.hpp | 14 +++---- .../lbboundaries/LBBoundary.hpp | 3 +- testsuite/python/lb_shear.py | 31 +++++++++----- testsuite/python/lb_stokes_sphere.py | 7 ++++ 8 files changed, 103 insertions(+), 21 deletions(-) create mode 100644 src/core/grid_based_algorithms/lbboundaries/LBBoundary.cpp diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 9105639a619..9f0827b0bb4 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -91,6 +91,7 @@ set(EspressoCore_SRC grid_based_algorithms/halo.cpp grid_based_algorithms/lattice.cpp grid_based_algorithms/lb_boundaries.cpp + grid_based_algorithms/lbboundaries/LBBoundary.cpp grid_based_algorithms/lb.cpp grid_based_algorithms/LbWalberla.cpp grid_based_algorithms/lb_walberla_instance.cpp diff --git a/src/core/grid_based_algorithms/LbWalberla.cpp b/src/core/grid_based_algorithms/LbWalberla.cpp index 353ce23bab2..d30868370d9 100644 --- a/src/core/grid_based_algorithms/LbWalberla.cpp +++ b/src/core/grid_based_algorithms/LbWalberla.cpp @@ -55,8 +55,8 @@ LbWalberla::LbWalberla(double viscosity, double density, double agrid, Utils::Vector3i grid_dimensions; for (int i = 0; i < 3; i++) { - if (fabs(floor(box_dimensions[i] / agrid) * agrid - box_dimensions[i]) > - std::numeric_limits::epsilon()) { + if (fabs((box_dimensions[i] / agrid) * agrid - box_dimensions[i]) > + 2 * std::numeric_limits::epsilon()) { throw std::runtime_error( "Box length not commensurate with agrid in direction " + std::to_string(i)); @@ -293,6 +293,24 @@ bool LbWalberla::add_force_at_pos(const Utils::Vector3d &pos, return true; } +boost::optional +LbWalberla::get_node_boundary_force(const Utils::Vector3i node) const { + auto bc = get_block_and_cell(node, true); // including ghosts + if (!bc) + return {boost::none}; + // Get boundary hanlindg + auto const &bh = + (*bc).block->getData(m_boundary_handling_id); + auto const &ff = (*bc).block->getData(m_flag_field_id); + if (!ff->isFlagSet((*bc).cell, ff->getFlag(UBB_flag))) + return {boost::none}; + + auto const uid = bh->getBoundaryUID(UBB_flag); + auto const &ubb = bh->getBoundaryCondition(uid); + return {to_vector3d( + ubb.getForce((*bc).cell.x(), (*bc).cell.y(), (*bc).cell.z()))}; +} + boost::optional LbWalberla::get_node_velocity(const Utils::Vector3i node) const { auto bc = get_block_and_cell(node); diff --git a/src/core/grid_based_algorithms/LbWalberla.hpp b/src/core/grid_based_algorithms/LbWalberla.hpp index 5b43e7443e3..5e7a3f809b6 100644 --- a/src/core/grid_based_algorithms/LbWalberla.hpp +++ b/src/core/grid_based_algorithms/LbWalberla.hpp @@ -112,7 +112,8 @@ class LbWalberla { using Flag_field_t = walberla::FlagField; using Pdf_field_t = walberla::lbm::PdfField; - using UBB_t = walberla::lbm::UBB; + using UBB_t = + walberla::lbm::UBB; using Boundary_handling_t = walberla::BoundaryHandling; @@ -173,6 +174,8 @@ class LbWalberla { bool add_force_at_pos(const Utils::Vector3d &position, const Utils::Vector3d &force); boost::optional + get_node_boundary_force(const Utils::Vector3i node) const; + boost::optional get_velocity_at_pos(const Utils::Vector3d &position) const; boost::optional get_density_at_pos(const Utils::Vector3d &position); // Utils::Vector3d get_stress_at_pos(const Utils::Vector3d& position); diff --git a/src/core/grid_based_algorithms/lbboundaries/LBBoundary.cpp b/src/core/grid_based_algorithms/lbboundaries/LBBoundary.cpp new file mode 100644 index 00000000000..5b6359e501a --- /dev/null +++ b/src/core/grid_based_algorithms/lbboundaries/LBBoundary.cpp @@ -0,0 +1,41 @@ +#include "LBBoundary.hpp" +#include "config.hpp" +#ifdef LB_BOUNDARIES +#include "grid_based_algorithms/lb_interface.hpp" +#include "grid_based_algorithms/lb_walberla_instance.hpp" +namespace LBBoundaries { +Utils::Vector3d LBBoundary::get_force() const { +#if defined(LB_BOUNDARIES) || defined(LB_BOUNDARIES_GPU) + if (lattice_switch == ActiveLB::WALBERLA) { + auto const grid = lb_walberla()->get_grid_dimensions(); + auto const agrid = lb_lbfluid_get_agrid(); + Utils::Vector3d force{0, 0, 0}; + for (auto index_and_pos : lb_walberla()->global_node_indices_positions()) { + // Convert to MD units + auto const index = index_and_pos.first; + auto const pos = index_and_pos.second * agrid; + if (calc_dist(pos) <= 0.) + for (int dx : {-1, 0, 1}) + for (int dy : {-1, 0, 1}) + for (int dz : {-1, 0, 1}) { + Utils::Vector3i shifted_index = + index + + Utils::Vector3i{dx * grid[0], dy * grid[1], dz * grid[2]}; + auto node_force = + lb_walberla()->get_node_boundary_force(shifted_index); + if (node_force) { + force += (*node_force); + } + } // Loop over cells + } // loop over lb cells + boost::mpi::all_reduce(comm_cart, force, std::plus()); + return force; + } else { + return lbboundary_get_force(this); + } +#else + throw std::runtime_error("Needs LB_BOUNDARIES or LB_BOUNDARIES_GPU."); +#endif +} +} // namespace LBBoundaries +#endif diff --git a/src/core/grid_based_algorithms/lbboundaries/LBBoundary.hpp b/src/core/grid_based_algorithms/lbboundaries/LBBoundary.hpp index 73ee5f7193a..76a24145a1d 100644 --- a/src/core/grid_based_algorithms/lbboundaries/LBBoundary.hpp +++ b/src/core/grid_based_algorithms/lbboundaries/LBBoundary.hpp @@ -49,6 +49,12 @@ class LBBoundary { m_shape->calculate_dist(pos, dist, vec); } + double calc_dist(const Utils::Vector3d &pos) const { + double tmp[3], dist; + calc_dist(pos, &dist, tmp); + return dist; + } + void set_shape(std::shared_ptr const &shape) { m_shape = shape; } @@ -64,13 +70,7 @@ class LBBoundary { Shapes::Shape const &shape() const { return *m_shape; } Utils::Vector3d &velocity() { return m_velocity; } Utils::Vector3d &force() { return m_force; } - Utils::Vector3d get_force() const { -#if defined(LB_BOUNDARIES) || defined(LB_BOUNDARIES_GPU) - return lbboundary_get_force(this); -#else - throw std::runtime_error("Needs LB_BOUNDARIES or LB_BOUNDARIES_GPU."); -#endif - } + Utils::Vector3d get_force() const; #ifdef EK_BOUNDARIES // TODO: ugly. Better would be a class EKBoundaries, // deriving from LBBoundaries, but that requires completely diff --git a/src/script_interface/lbboundaries/LBBoundary.hpp b/src/script_interface/lbboundaries/LBBoundary.hpp index 08a6ea67d24..4ee38681a68 100644 --- a/src/script_interface/lbboundaries/LBBoundary.hpp +++ b/src/script_interface/lbboundaries/LBBoundary.hpp @@ -71,7 +71,8 @@ class LBBoundary : public AutoParameters { const auto agrid = lb_lbfluid_get_agrid(); const auto tau = lb_lbfluid_get_tau(); const double unit_conversion = - agrid * agrid * agrid * agrid / rho / tau / tau; + agrid / tau / tau * agrid * agrid * agrid; + // agrid * agrid * agrid * agrid / rho / tau / tau; return m_lbboundary->get_force() * unit_conversion; } return none; diff --git a/testsuite/python/lb_shear.py b/testsuite/python/lb_shear.py index 917818c8046..bb3494d83d4 100644 --- a/testsuite/python/lb_shear.py +++ b/testsuite/python/lb_shear.py @@ -29,10 +29,10 @@ """ -AGRID = 1 +AGRID = .5 VISC = 5.2 -DENS = 2.3 -TIME_STEP = 0.02 +DENS = 2.3 +TIME_STEP = 0.03 # Box size will be H +2 AGRID to make room for walls. # The number of grid cells should be devisible by four and 3 in all directions # for testing on multiple mpi nodes. @@ -139,6 +139,24 @@ def check_profile(self, shear_plane_normal, shear_direction): np.copy(v_measured), np.copy(v_expected[j]) * shear_direction, atol=3E-3) + print(wall1.get_force()) + print(wall2.get_force()) + p_eq = DENS * AGRID**2 / TIME_STEP**2 / 3 + print( + shear_direction * + SHEAR_VELOCITY / + H * + W**2 * + VISC + + p_eq * + W**2 * + shear_plane_normal) + np.testing.assert_allclose( + np.copy(wall1.get_force()), + -np.copy(wall2.get_force()), + atol=1E-4) + np.testing.assert_allclose(np.copy(wall1.get_force()), + shear_direction * SHEAR_VELOCITY / H * W**2 * VISC, atol=2E-4) # Test steady state stress tensor on a node p_eq = DENS * AGRID**2 / TIME_STEP**2 / 3 p_expected = np.diag((p_eq, p_eq, p_eq)) @@ -150,13 +168,6 @@ def check_profile(self, shear_plane_normal, shear_direction): np.testing.assert_allclose(node_stress, p_expected, atol=1E-5, rtol=5E-3) - np.testing.assert_allclose( - np.copy(wall1.get_force()), - -np.copy(wall2.get_force()), - atol=1E-4) - np.testing.assert_allclose(np.copy(wall1.get_force()), - shear_direction * SHEAR_VELOCITY / H * W**2 * VISC, atol=2E-4) - def test(self): x = np.array((1, 0, 0), dtype=float) y = np.array((0, 1, 0), dtype=float) diff --git a/testsuite/python/lb_stokes_sphere.py b/testsuite/python/lb_stokes_sphere.py index 3c2c7a4db99..c1dc8a39b0f 100644 --- a/testsuite/python/lb_stokes_sphere.py +++ b/testsuite/python/lb_stokes_sphere.py @@ -123,5 +123,12 @@ def setUp(self): self.lbf = espressomd.lb.LBFluid(**LB_PARAMS) +@utx.skipIfMissingFeatures(['LB_WALBERLA', 'EXTERNAL_FORCES']) +class LBWalberlaStokes(ut.TestCase, Stokes): + + def setUp(self): + self.lbf = espressomd.lb.LBFluidWalberla(**LB_PARAMS) + + if __name__ == "__main__": ut.main()