diff --git a/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp b/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp index c9fbe803b03..13b430e612c 100644 --- a/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp +++ b/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp @@ -63,9 +63,9 @@ #include <walberla_bridge/lattice_boltzmann/LeesEdwardsPack.hpp> #include <utils/Vector.hpp> +#include <utils/index.hpp> #include <utils/interpolation/bspline_3d.hpp> #include <utils/math/make_lin_space.hpp> -#include <utils/index.hpp> #include <array> #include <bitset> @@ -819,9 +819,10 @@ class LBWalberlaImpl : public LBWalberlaBase { auto const &lattice = get_lattice(); auto const n_ghost_layers = lattice.get_ghost_layers(); auto const blocks = lattice.get_blocks(); - if ((shear_direction == 0u and blocks->getXSize() != 1u) or (shear_direction == 2u and blocks->getZSize() != 1u)) { - throw std::domain_error( - "Lees-Edwards LB doesn't support domain decomposition along the shear direction."); + if ((shear_direction == 0u and blocks->getXSize() != 1u) or + (shear_direction == 2u and blocks->getZSize() != 1u)) { + throw std::domain_error("Lees-Edwards LB doesn't support domain " + "decomposition along the shear direction."); } auto const agrid = FloatType_c(lattice.get_grid_dimensions()[shear_plane_normal]); @@ -917,16 +918,17 @@ class LBWalberlaImpl : public LBWalberlaBase { template <typename F> void mapping_block_to_local(std::optional<CellInterval> const &bci, - std::optional<CellInterval> const &ci, - Utils::Vector3i const &block_offset, - Utils::Vector3i const &lower_corner, - F&& func) const { - auto const local_grid = Utils::Vector3i{{ci->max().x() - ci->min().x() + 1, - ci->max().y() - ci->min().y() + 1, - ci->max().z() - ci->min().z() + 1}}; - auto const block_grid = Utils::Vector3i{{bci->max().x() - bci->min().x() + 1, - bci->max().y() - bci->min().y() + 1, - bci->max().z() - bci->min().z() + 1}}; + std::optional<CellInterval> const &ci, + Utils::Vector3i const &block_offset, + Utils::Vector3i const &lower_corner, + F &&func) const { + auto const local_grid = Utils::Vector3i{ + {ci->max().x() - ci->min().x() + 1, ci->max().y() - ci->min().y() + 1, + ci->max().z() - ci->min().z() + 1}}; + auto const block_grid = + Utils::Vector3i{{bci->max().x() - bci->min().x() + 1, + bci->max().y() - bci->min().y() + 1, + bci->max().z() - bci->min().z() + 1}}; auto const lower_cell = bci->min(); auto const upper_cell = bci->max(); // In the loop, x,y,z are in block coordinates @@ -935,20 +937,17 @@ class LBWalberlaImpl : public LBWalberlaBase { // to block coordinates for (auto x = lower_cell.x(); x <= upper_cell.x(); ++x) { for (auto y = lower_cell.y(); y <= upper_cell.y(); ++y) { - for (auto z = lower_cell.z(); z <= upper_cell.z(); ++z) { - auto const node = block_offset + Utils::Vector3i{{x, y, z}}; - auto const local_index = Utils::get_linear_index(node[0] - lower_corner[0], - node[1] - lower_corner[1], - node[2] - lower_corner[2], - local_grid, - Utils::MemoryOrder::ROW_MAJOR); - auto const block_index = Utils::get_linear_index(x - lower_cell.x(), - y - lower_cell.y(), - z - lower_cell.z(), - block_grid, - Utils::MemoryOrder::ROW_MAJOR); - func(block_index, local_index, node); - } + for (auto z = lower_cell.z(); z <= upper_cell.z(); ++z) { + auto const node = block_offset + Utils::Vector3i{{x, y, z}}; + auto const local_index = Utils::get_linear_index( + node[0] - lower_corner[0], node[1] - lower_corner[1], + node[2] - lower_corner[2], local_grid, + Utils::MemoryOrder::ROW_MAJOR); + auto const block_index = Utils::get_linear_index( + x - lower_cell.x(), y - lower_cell.y(), z - lower_cell.z(), + block_grid, Utils::MemoryOrder::ROW_MAJOR); + func(block_index, local_index, node); + } } } } @@ -973,26 +972,24 @@ class LBWalberlaImpl : public LBWalberlaBase { assert(values.size() == 3u * bci->numCells()); values_size += 3u * bci->numCells(); - auto func = [&values, &out, this] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - if (m_boundary->node_is_boundary(node)) { - auto const &vec = - m_boundary->get_node_value_at_boundary(node); - for (uint_t f = 0u; f < 3u; ++f) { - out[static_cast<unsigned int>(3u * local_index + f)] = - double_c(vec[f]); - } - } else { - for (uint_t f = 0u; f < 3u; ++f) { - out[static_cast<unsigned int>(3u * local_index + f)] = - double_c(values[static_cast<unsigned int>( - 3u * block_index + f)]); - } - } - }; - - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + auto func = [&values, &out, this](uint_t block_index, + uint_t local_index, + Utils::Vector3i node) { + if (m_boundary->node_is_boundary(node)) { + auto const &vec = m_boundary->get_node_value_at_boundary(node); + for (uint_t f = 0u; f < 3u; ++f) { + out[static_cast<unsigned int>(3u * local_index + f)] = + double_c(vec[f]); + } + } else { + for (uint_t f = 0u; f < 3u; ++f) { + out[static_cast<unsigned int>(3u * local_index + f)] = double_c( + values[static_cast<unsigned int>(3u * block_index + f)]); + } + } + }; + + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); } } assert(values_size == 3u * ci->numCells()); @@ -1022,17 +1019,17 @@ class LBWalberlaImpl : public LBWalberlaBase { std::vector<FloatType> values = std::vector<FloatType>( static_cast<unsigned int>(3u * bci->numCells())); - auto func = [&values, &velocity] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - for (uint_t f = 0u; f < 3u; ++f) { - values[static_cast<unsigned int>(3u * block_index + f)] = - numeric_cast<FloatType>( - velocity[static_cast<unsigned int>(3u * local_index + f)]); - } - }; - - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + auto func = [&values, &velocity](uint_t block_index, + uint_t local_index, + Utils::Vector3i node) { + for (uint_t f = 0u; f < 3u; ++f) { + values[static_cast<unsigned int>(3u * block_index + f)] = + numeric_cast<FloatType>(velocity[static_cast<unsigned int>( + 3u * local_index + f)]); + } + }; + + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); lbm::accessor::Velocity::set(pdf_field, vel_field, force_field, values, *bci); } @@ -1267,16 +1264,16 @@ class LBWalberlaImpl : public LBWalberlaBase { auto const values = lbm::accessor::Vector::get(field, *bci); assert(values.size() == 3u * bci->numCells()); - auto func = [&values, &out, this] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - for (uint_t f = 0u; f < 3u; ++f) { - out[static_cast<unsigned int>(3u * local_index + f)] = - values[static_cast<unsigned int>(3u * block_index + f)]; - } - }; + auto func = [&values, &out, this](uint_t block_index, + uint_t local_index, + Utils::Vector3i node) { + for (uint_t f = 0u; f < 3u; ++f) { + out[static_cast<unsigned int>(3u * local_index + f)] = + values[static_cast<unsigned int>(3u * block_index + f)]; + } + }; - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); } } } @@ -1305,17 +1302,16 @@ class LBWalberlaImpl : public LBWalberlaBase { std::vector<FloatType> values = std::vector<FloatType>( static_cast<unsigned int>(3u * bci->numCells())); - auto func = [&values, &force] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - for (uint_t f = 0u; f < 3u; ++f) { - values[static_cast<unsigned int>(3u * block_index + f)] = - numeric_cast<FloatType>( - force[static_cast<unsigned int>(3u * local_index + f)]); - } - }; - - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + auto func = [&values, &force](uint_t block_index, uint_t local_index, + Utils::Vector3i node) { + for (uint_t f = 0u; f < 3u; ++f) { + values[static_cast<unsigned int>(3u * block_index + f)] = + numeric_cast<FloatType>( + force[static_cast<unsigned int>(3u * local_index + f)]); + } + }; + + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); lbm::accessor::Force::set(pdf_field, vel_field, force_field, values, *bci); } @@ -1384,17 +1380,17 @@ class LBWalberlaImpl : public LBWalberlaBase { auto const values = lbm::accessor::Population::get(pdf_field, *bci); assert(values.size() == stencil_size() * bci->numCells()); - auto func = [&values, &out, this] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - for (uint_t f = 0u; f < stencil_size(); ++f) { - out[static_cast<unsigned int>(stencil_size() * local_index + f)] = - values[static_cast<unsigned int>( - stencil_size() * block_index + f)]; - } - }; - - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + auto func = [&values, &out, this](uint_t block_index, + uint_t local_index, + Utils::Vector3i node) { + for (uint_t f = 0u; f < stencil_size(); ++f) { + out[static_cast<unsigned int>(stencil_size() * local_index + f)] = + values[static_cast<unsigned int>( + stencil_size() * block_index + f)]; + } + }; + + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); } } } @@ -1421,19 +1417,18 @@ class LBWalberlaImpl : public LBWalberlaBase { std::vector<FloatType> values = std::vector<FloatType>( static_cast<unsigned int>(stencil_size() * bci->numCells())); - auto func = [&values, &population, this] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - for (uint_t f = 0u; f < stencil_size(); ++f) { - values[static_cast<unsigned int>( - stencil_size() * block_index + f)] = - numeric_cast<FloatType>( - population[static_cast<unsigned int>( - stencil_size() * local_index + f)]); - } - }; - - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + auto func = [&values, &population, this](uint_t block_index, + uint_t local_index, + Utils::Vector3i node) { + for (uint_t f = 0u; f < stencil_size(); ++f) { + values[static_cast<unsigned int>(stencil_size() * block_index + + f)] = + numeric_cast<FloatType>(population[static_cast<unsigned int>( + stencil_size() * local_index + f)]); + } + }; + + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); lbm::accessor::Population::set(pdf_field, vel_field, force_field, values, *bci); } @@ -1486,13 +1481,12 @@ class LBWalberlaImpl : public LBWalberlaBase { auto const values = lbm::accessor::Density::get(pdf_field, *bci); assert(values.size() == bci->numCells()); - auto func = [&values, &out] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - out[local_index] = values[block_index]; - }; + auto func = [&values, &out](uint_t block_index, uint_t local_index, + Utils::Vector3i node) { + out[local_index] = values[block_index]; + }; - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); } } } @@ -1516,13 +1510,13 @@ class LBWalberlaImpl : public LBWalberlaBase { std::vector<FloatType> values = std::vector<FloatType>(bci->numCells()); - auto func = [&values, &density] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - values[block_index] = numeric_cast<FloatType>(density[local_index]); - }; + auto func = [&values, &density](uint_t block_index, + uint_t local_index, + Utils::Vector3i node) { + values[block_index] = numeric_cast<FloatType>(density[local_index]); + }; - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); lbm::accessor::Density::set(pdf_field, values, *bci); } } @@ -1569,18 +1563,17 @@ class LBWalberlaImpl : public LBWalberlaBase { if (auto const bci = get_block_interval(lower_corner, upper_corner, block_offset, block)) { - auto func = [&out, this] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - if (m_boundary->node_is_boundary(node)) { - out[local_index] = - to_vector3d(m_boundary->get_node_value_at_boundary(node)); - } else { - out[local_index] = std::nullopt; - } - }; - - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + auto func = [&out, this](uint_t block_index, uint_t local_index, + Utils::Vector3i node) { + if (m_boundary->node_is_boundary(node)) { + out[local_index] = + to_vector3d(m_boundary->get_node_value_at_boundary(node)); + } else { + out[local_index] = std::nullopt; + } + }; + + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); } } assert(out.size() == ci->numCells()); @@ -1603,21 +1596,21 @@ class LBWalberlaImpl : public LBWalberlaBase { if (auto const bci = get_block_interval(lower_corner, upper_corner, block_offset, block)) { - auto func = [&lattice, &block, &velocity, this] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - auto const bc = get_block_and_cell(lattice, node, false); - assert(bc->block->getAABB() == block.getAABB()); - auto const &opt = velocity[local_index]; - if (opt) { - m_boundary->set_node_value_at_boundary( - node, to_vector3<FloatType>(*opt), *bc); - } else { - m_boundary->remove_node_from_boundary(node, *bc); - } - }; - - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + auto func = [&lattice, &block, &velocity, + this](uint_t block_index, uint_t local_index, + Utils::Vector3i node) { + auto const bc = get_block_and_cell(lattice, node, false); + assert(bc->block->getAABB() == block.getAABB()); + auto const &opt = velocity[local_index]; + if (opt) { + m_boundary->set_node_value_at_boundary( + node, to_vector3<FloatType>(*opt), *bc); + } else { + m_boundary->remove_node_from_boundary(node, *bc); + } + }; + + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); } } } @@ -1668,13 +1661,12 @@ class LBWalberlaImpl : public LBWalberlaBase { if (auto const bci = get_block_interval(lower_corner, upper_corner, block_offset, block)) { - auto func = [&out, this] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - out[local_index] = m_boundary->node_is_boundary(node); - }; + auto func = [&out, this](uint_t block_index, uint_t local_index, + Utils::Vector3i node) { + out[local_index] = m_boundary->node_is_boundary(node); + }; - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); } } assert(out.size() == ci->numCells()); @@ -1744,18 +1736,18 @@ class LBWalberlaImpl : public LBWalberlaBase { auto values = lbm::accessor::PressureTensor::get(pdf_field, *bci); assert(values.size() == 9u * bci->numCells()); - auto func = [&values, &out, this] (uint_t block_index, - uint_t local_index, - Utils::Vector3i node) { - pressure_tensor_correction(std::span<FloatType, 9ul>( - &values[static_cast<unsigned int>(9u * block_index)], 9ul)); - for (uint_t f = 0u; f < 9u; ++f) { - out[static_cast<unsigned int>(9u * local_index + f)] = - values[static_cast<unsigned int>(9u * block_index + f)]; - } - }; - - mapping_block_to_local(bci, ci, block_offset, lower_corner, func); + auto func = [&values, &out, this](uint_t block_index, + uint_t local_index, + Utils::Vector3i node) { + pressure_tensor_correction(std::span<FloatType, 9ul>( + &values[static_cast<unsigned int>(9u * block_index)], 9ul)); + for (uint_t f = 0u; f < 9u; ++f) { + out[static_cast<unsigned int>(9u * local_index + f)] = + values[static_cast<unsigned int>(9u * block_index + f)]; + } + }; + + mapping_block_to_local(bci, ci, block_offset, lower_corner, func); } } } diff --git a/src/walberla_bridge/tests/CMakeLists.txt b/src/walberla_bridge/tests/CMakeLists.txt index fa3ddc09946..7b3a85ab1b9 100644 --- a/src/walberla_bridge/tests/CMakeLists.txt +++ b/src/walberla_bridge/tests/CMakeLists.txt @@ -30,8 +30,7 @@ function(ESPRESSO_ADD_TEST) espresso::walberla_codegen_cuda) endif() if(${TEST_SRC} MATCHES ".*\.cu$") - target_link_libraries(${TEST_NAME} PRIVATE espresso::walberla::cuda_flags - espresso::walberla_cuda) + target_link_libraries(${TEST_NAME} PRIVATE espresso::walberla::cuda_flags) else() target_link_libraries(${TEST_NAME} PRIVATE espresso::walberla::cpp_flags) endif()