From f67262937585a84ae7171471025768b5fbd88401 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Fri, 9 Aug 2019 16:35:56 +0200 Subject: [PATCH] WIP: Boundary support for Walberla --- src/core/grid_based_algorithms/LbWalberla.cpp | 86 +++++++++++++++---- src/core/grid_based_algorithms/LbWalberla.hpp | 24 ++++-- .../grid_based_algorithms/lb_boundaries.cpp | 46 +++++++++- src/core/unit_tests/LbWalberla.cpp | 52 +++++++++++ testsuite/python/lb_shear.py | 27 ++++-- 5 files changed, 202 insertions(+), 33 deletions(-) diff --git a/src/core/grid_based_algorithms/LbWalberla.cpp b/src/core/grid_based_algorithms/LbWalberla.cpp index 630c7aca6c6..63931f8fede 100644 --- a/src/core/grid_based_algorithms/LbWalberla.cpp +++ b/src/core/grid_based_algorithms/LbWalberla.cpp @@ -4,6 +4,7 @@ #include "communication.hpp" #include "lb_walberla_instance.hpp" #include "utils/Vector.hpp" +#include "utils/math/make_lin_space.hpp" #include "utils/mpi/gatherv.hpp" #include "blockforest/Initialization.h" @@ -69,7 +70,9 @@ LbWalberla::LbWalberla(double viscosity, double density, double agrid, } } m_grid_dimensions = grid_dimensions; - printf("grid: %d %d %d, node: %d %d %d\n",grid_dimensions[0],grid_dimensions[1],grid_dimensions[2],node_grid[0],node_grid[1],node_grid[2]); + printf("grid: %d %d %d, node: %d %d %d\n", grid_dimensions[0], + grid_dimensions[1], grid_dimensions[2], node_grid[0], node_grid[1], + node_grid[2]); m_blocks = blockforest::createUniformBlockGrid( uint_c(node_grid[0]), // blocks in x direction @@ -108,7 +111,7 @@ LbWalberla::LbWalberla(double viscosity, double density, double agrid, LB_boundary_handling(m_flag_field_id, m_pdf_field_id), "boundary handling"); - empty_flags_for_boundary(m_blocks, m_boundary_handling_id); + clear_boundaries(); // 1 timestep each time the integrate function is called m_time_loop = @@ -125,20 +128,19 @@ LbWalberla::LbWalberla(double viscosity, double density, double agrid, ResetForce>( m_pdf_field_id, m_force_field_id, m_force_field_from_md_id, m_boundary_handling_id); - m_time_loop->add() - << // timeloop::BeforeFunction(communication, "communication") - timeloop::Sweep( - Boundary_handling_t::getBlockSweep(m_boundary_handling_id), - "boundary handling"); + m_time_loop->add() << timeloop::BeforeFunction(communication, "communication") + << timeloop::Sweep(Boundary_handling_t::getBlockSweep( + m_boundary_handling_id), + "boundary handling"); m_time_loop->add() << timeloop::Sweep(makeSharedSweep(m_reset_force), "Reset force fields"); - m_time_loop->add() - << timeloop::AfterFunction(communication, "communication") - << timeloop::Sweep( - domain_decomposition::makeSharedSweep( - lbm::makeCellwiseSweep( - m_pdf_field_id, m_flag_field_id, Fluid_flag)), - "LB stream & collide"); + // m_time_loop->add() + // << timeloop::AfterFunction(communication, "communication") + m_time_loop->add() << timeloop::Sweep( + domain_decomposition::makeSharedSweep( + lbm::makeCellwiseSweep( + m_pdf_field_id, m_flag_field_id, Fluid_flag)), + "LB stream & collide"); m_force_distributor_id = field::addDistributor( @@ -181,7 +183,7 @@ LbWalberla::get_block_and_cell(const Utils::Vector3i &node, // Try to find a block which has the cell as ghost layer for (auto b = m_blocks->begin(); b != m_blocks->end(); ++b) { if (b->getAABB() - .getExtended(m_skin / m_agrid) + .getExtended(1) .contains(real_c(node[0]), real_c(node[1]), real_c(node[2]))) { block = &(*b); break; @@ -222,7 +224,7 @@ IBlock *LbWalberla::get_block(const Utils::Vector3d &pos, bool LbWalberla::set_node_velocity_at_boundary(const Utils::Vector3i node, const Utils::Vector3d v) { - auto bc = get_block_and_cell(node); + auto bc = get_block_and_cell(node,true); // Return if we don't have the cell. if (!bc) return false; @@ -423,4 +425,56 @@ bool LbWalberla::pos_in_local_domain(const Utils::Vector3d &pos) const { return (block != nullptr); } +std::vector> +LbWalberla::node_indices_positions() { + + std::vector> res; + for (auto block = m_blocks->begin(); block != m_blocks->end(); ++block) { + auto left = block->getAABB().min(); + // Lattice constant is 1, node centers are offset by .5 + Utils::Vector3d pos_offset = + to_vector3d(left) + Utils::Vector3d::broadcast(.5); + + // Lattice constant is 1, so cast left corner position to ints + Utils::Vector3i index_offset = + Utils::Vector3i{int(left[0]), int(left[1]), int(left[2])}; + + // Get field data which knows about the indices + // In the loop, x,y,z are in block-local coordinates + auto pdf_field = block->getData(m_pdf_field_id); + for (int x : Utils::make_lin_space(int(0), int(pdf_field->xSize()), + int(pdf_field->xSize()), false)) { + for (int y : Utils::make_lin_space(int(0), int(pdf_field->ySize()), + int(pdf_field->ySize()), false)) { + for (int z : Utils::make_lin_space(int(0), int(pdf_field->zSize()), + int(pdf_field->zSize()), false)) { + res.push_back( + {index_offset + Utils::Vector3i{x, y, z}, + pos_offset + Utils::Vector3d{double(x), double(y), double(z)}}); + } + } + } + } + return res; +} + +std::vector> +LbWalberla::global_node_indices_positions() { + + std::vector> res; + for (int x : Utils::make_lin_space(int(0), int(m_grid_dimensions[0]), + int(m_grid_dimensions[0]), false)) { + for (int y : Utils::make_lin_space(int(0), int(m_grid_dimensions[1]), + int(m_grid_dimensions[1]), false)) { + for (int z : Utils::make_lin_space(int(0), int(m_grid_dimensions[2]), + int(m_grid_dimensions[2]), false)) { + res.push_back( + {Utils::Vector3i{x, y, z}, + Utils::Vector3d::broadcast(.5) +Utils::Vector3d{double(x), double(y), double(z)}}); + } + } + } + return res; +} + #endif diff --git a/src/core/grid_based_algorithms/LbWalberla.hpp b/src/core/grid_based_algorithms/LbWalberla.hpp index 109c75bb9b9..5b43e7443e3 100644 --- a/src/core/grid_based_algorithms/LbWalberla.hpp +++ b/src/core/grid_based_algorithms/LbWalberla.hpp @@ -130,7 +130,10 @@ class LbWalberla { Flag_field_t *flag_field = block->getData(m_flag_field_id); Pdf_field_t *pdf_field = block->getData(m_pdf_field_id); - const uint8_t fluid = flag_field->registerFlag(Fluid_flag); + // const uint8_t fluid = flag_field->registerFlag(Fluid_flag); + const auto fluid = flag_field->flagExists(Fluid_flag) + ? flag_field->getFlag(Fluid_flag) + : flag_field->registerFlag(Fluid_flag); return new Boundary_handling_t( "boundary handling", flag_field, fluid, @@ -155,7 +158,7 @@ class LbWalberla { void integrate(); std::pair get_local_domain() { // We only have one block per mpi rank - assert(m_blocks->begin()++ == m_blocks->end()); + assert(++(m_blocks->begin()) == m_blocks->end()); auto const ab = m_blocks->begin()->getAABB(); return {to_vector3d(ab.min()), to_vector3d(ab.max())}; @@ -261,19 +264,22 @@ class LbWalberla { walberla::field::TrilinearFieldInterpolator; - void empty_flags_for_boundary( - std::shared_ptr &blocks, - const walberla::BlockDataID &boundary_handling_id) { +public: + std::vector> + node_indices_positions(); + std::vector> + global_node_indices_positions(); + void clear_boundaries() { using namespace walberla; const CellInterval &domain_bb_in_global_cell_coordinates = - blocks->getDomainCellBB(); - for (auto block = blocks->begin(); block != blocks->end(); ++block) { + m_blocks->getDomainCellBB(); + for (auto block = m_blocks->begin(); block != m_blocks->end(); ++block) { Boundary_handling_t *boundary_handling = - block->getData(boundary_handling_id); + block->getData(m_boundary_handling_id); CellInterval domain_bb(domain_bb_in_global_cell_coordinates); - blocks->transformGlobalToBlockLocalCellInterval(domain_bb, *block); + m_blocks->transformGlobalToBlockLocalCellInterval(domain_bb, *block); boundary_handling->fillWithDomain(domain_bb); } diff --git a/src/core/grid_based_algorithms/lb_boundaries.cpp b/src/core/grid_based_algorithms/lb_boundaries.cpp index c1462c67ba1..6e607ec4d63 100644 --- a/src/core/grid_based_algorithms/lb_boundaries.cpp +++ b/src/core/grid_based_algorithms/lb_boundaries.cpp @@ -34,6 +34,7 @@ #include "grid_based_algorithms/lattice.hpp" #include "grid_based_algorithms/lb.hpp" #include "grid_based_algorithms/lb_interface.hpp" +#include "grid_based_algorithms/lb_walberla_instance.hpp" #include "grid_based_algorithms/lbgpu.hpp" #include "lbboundaries/LBBoundary.hpp" @@ -271,9 +272,52 @@ void lb_init_boundaries() { } } #endif - } + } else if (lattice_switch == ActiveLB::WALBERLA) { +#ifdef LB_WALBERLA +#if defined(LB_BOUNDARIES) + Utils::Vector3i offset; + int the_boundary = -1; + + lb_walberla()->clear_boundaries(); + + auto const agrid = lb_lbfluid_get_agrid(); + + 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; + + int n = 0; + for (auto it = lbboundaries.begin(); it != lbboundaries.end(); + ++it, ++n) { + double dist; + Utils::Vector3d tmp; + (**it).calc_dist(pos, &dist, tmp.data()); + + if (dist <= 0) { + + // Set boundaries on the ghost layers + auto const grid = lb_walberla()->get_grid_dimensions(); + 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]}; + if (lb_walberla()->set_node_velocity_at_boundary( + shifted_index, (**it).velocity() / lb_lbfluid_get_lattice_speed())) { + printf("Boundary: %d, %d %d %d, %g %g %g\n", n, shifted_index[0], shifted_index[1], + shifted_index[2], pos[0], pos[1], pos[2]); + } + } + break; + } // if dist <=0 + } // loop over boundaries + } // Loop over cells + } // lattice switch is WALBERLA +#endif +#endif } + Utils::Vector3d lbboundary_get_force(LBBoundary const *lbb) { Utils::Vector3d force{}; #if defined(LB_BOUNDARIES) || defined(LB_BOUNDARIES_GPU) diff --git a/src/core/unit_tests/LbWalberla.cpp b/src/core/unit_tests/LbWalberla.cpp index c0665eb86e2..692ef67f4fd 100644 --- a/src/core/unit_tests/LbWalberla.cpp +++ b/src/core/unit_tests/LbWalberla.cpp @@ -75,6 +75,58 @@ BOOST_AUTO_TEST_CASE(boundary) { } } +// BOOST_AUTO_TEST_CASE(boundary_flow_single_node) { +// Vector3d vel = {0.2, 3.8, 0}; +// LbWalberla lb = LbWalberla(viscosity, density, agrid, tau, box_dimensions, +// node_grid, skin); +// Vector3i node{5,5,5}; +// if (lb.node_in_local_domain(node)) { +// BOOST_CHECK(lb.set_node_velocity_at_boundary(node, vel)); +// } +// for( int i=0;i<10;i++) { +// lb.integrate(); +// } +// for (int j=0;j<9;j++) { +// auto v = lb.get_node_velocity(Vector3i{j,node[1],node[2]}); +// if (v) { +// if (j!=node[0]) { +// BOOST_CHECK((*v)[0]>1E-4); +// BOOST_CHECK((*v)[1]>1E-4); +// BOOST_CHECK(fabs((*v)[2])<1E-12); +// printf("%d,%g %g %g\n",j,(*v)[0],(*v)[1],(*v)[2]); +// } +// } +// } +//} + +BOOST_AUTO_TEST_CASE(boundary_flow_shear) { + Vector3d vel = {0.2, -0.3, 0}; + LbWalberla lb = LbWalberla(viscosity, density, agrid, tau, box_dimensions, + node_grid, skin); + for (int x = 1; x < grid_dimensions[0] - 1; x++) { + for (int y = 1; y < grid_dimensions[1] - 1; y++) { + Vector3i node{x, y, 1}; + if (lb.node_in_local_domain(node)) { + BOOST_CHECK(lb.set_node_velocity_at_boundary(node, vel)); + } + node[2] = grid_dimensions[2] - 4; + if (lb.node_in_local_domain(node)) { + BOOST_CHECK(lb.set_node_velocity_at_boundary(node, -vel)); + } + } + } + + for (int i = 0; i < 150; i++) { + lb.integrate(); + } + for (int j = 0; j < grid_dimensions[2]; j++) { + auto v = lb.get_node_velocity(Vector3i{0, 0, j}); + if (v) { + printf("%d,%g %g %g\n", j, (*v)[0], (*v)[1], (*v)[2]); + } + } +} + BOOST_AUTO_TEST_CASE(velocity) { LbWalberla lb = LbWalberla(viscosity, density, agrid, tau, box_dimensions, node_grid, skin); diff --git a/testsuite/python/lb_shear.py b/testsuite/python/lb_shear.py index 05c1d313548..ac65f8099b9 100644 --- a/testsuite/python/lb_shear.py +++ b/testsuite/python/lb_shear.py @@ -29,8 +29,8 @@ """ -AGRID = 0.6 -VISC = 3.2 +AGRID = 1 +VISC = 5.2 DENS = 2.3 TIME_STEP = 0.02 # Box size will be H +2 AGRID to make room for walls. @@ -79,7 +79,7 @@ def shear_flow(x, t, nu, v, h, k_max): return v * u -class LBShearCommon(object): +class LBShearCommon: """Base class of the test that holds the test logic.""" lbf = None @@ -115,22 +115,26 @@ def check_profile(self, shear_plane_normal, shear_direction): self.system.lbboundaries.add(wall2) t0 = self.system.time - sample_points = int(H / AGRID - 1) + sample_points = np.arange(H/AGRID+2) for i in range(9): self.system.integrator.run(50) v_expected = shear_flow( - x=(np.arange(0, sample_points) + .5) * AGRID, + x=(sample_points -.5) * AGRID, t=self.system.time - t0, nu=VISC, v=SHEAR_VELOCITY, h=H, k_max=100) - for j in range(1, sample_points): - ind = np.max(((1, 1, 1), shear_plane_normal * j + 1), 0) + + # We omit the boundary nodes themselves, which have undefined velocities + for j in np.array(sample_points,dtype=int)[1:-1]: + ind = (1, 1, 1) ind = np.array(ind, dtype=int) + ind[np.argmax(shear_plane_normal)]=j v_measured = self.lbf[ind[0], ind[1], ind[2]].velocity + print(j,v_measured,v_expected[j]) np.testing.assert_allclose( np.copy(v_measured), np.copy(v_expected[j]) * shear_direction, atol=3E-3) @@ -183,6 +187,15 @@ class LBGPUShear(ut.TestCase, LBShearCommon): def setUp(self): self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS) +@utx.skipIfMissingFeatures(['LB_WALBERLA', 'LB_BOUNDARIES']) +class LBWalberlaShear(ut.TestCase, LBShearCommon): + + """Test for the Walberla implementation of the LB.""" + + def setUp(self): + self.lbf = espressomd.lb.LBFluidWalberla(**LB_PARAMS) + + if __name__ == '__main__': ut.main()