From 91bff8bad5d9b8e3b139c16861052454eb1f8ec7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 8 Jul 2024 18:07:20 +0200 Subject: [PATCH] Refactor Lees-Edwards code Capture Lees-Edwards parameters by reference, remove mention of Python in the core, improve code coverage. --- src/core/lb/particle_coupling.cpp | 19 +++++---- .../InterpolateAndShiftAtBoundary.hpp | 39 +++++++------------ .../lb_lees_edwards_particle_coupling.py | 31 +++++---------- 3 files changed, 34 insertions(+), 55 deletions(-) diff --git a/src/core/lb/particle_coupling.cpp b/src/core/lb/particle_coupling.cpp index cfe3ce622b..9bfa71e7ed 100644 --- a/src/core/lb/particle_coupling.cpp +++ b/src/core/lb/particle_coupling.cpp @@ -123,12 +123,12 @@ static void positions_in_halo_impl(Utils::Vector3d const &pos_folded, pos_folded + Utils::hadamard_product(box_geo.length(), shift); if (box_geo.type() == BoxType::LEES_EDWARDS) { - // note: Python-style modulo: 0<=offset positions_in_halo(Utils::Vector3d const &pos, return res; } -Utils::Vector3d -lees_edwards_vel_shift(const Utils::Vector3d &pos_shifted_by_box_l, - const Utils::Vector3d &orig_pos, - const BoxGeometry &box_geo) { +static auto lees_edwards_vel_shift(Utils::Vector3d const &pos_shifted_by_box_l, + Utils::Vector3d const &orig_pos, + BoxGeometry const &box_geo) { Utils::Vector3d vel_shift{{0., 0., 0.}}; if (box_geo.type() == BoxType::LEES_EDWARDS) { auto le = box_geo.lees_edwards_bc(); @@ -291,7 +290,7 @@ void ParticleCoupling::kernel(std::vector const &particles) { auto const random_force = get_noise_term(p); force_on_particle = drag_force + random_force; ++it_interpolated_velocities; - ++it_positions_velocity_coupling++; + ++it_positions_velocity_coupling; } auto force_on_fluid = -force_on_particle; diff --git a/src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp b/src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp index d24dd54b46..489e81be9d 100644 --- a/src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp +++ b/src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp @@ -31,16 +31,6 @@ #include #include -namespace { -template T python_modulo(T a, T b) { - T res = std::fmod(a, b); - if (res < 0) - return res + b; - else - return res; -} -} // namespace - namespace walberla { /** @@ -118,26 +108,22 @@ class InterpolateAndShiftAtBoundary { // the target auto const prefactor = ((slab_dir == m_slab_max) ? FloatType{-1} : FloatType{1}); - auto const offset = get_pos_offset() * prefactor; - auto const folded_offset = python_modulo(offset, length); + auto const offset = static_cast(get_pos_offset()) * prefactor; + auto const folded_offset = modulo(offset, length); // 0<=folded_offset(cell[dir]) + folded_offset; + auto const folded_source_pos = modulo(source_pos, length); + // 0 <= folded_source_pos < length source1[dir] = cell_idx_c(std::floor(folded_source_pos)); - // 0<=source1[dir]get(cell, q) = field->get(source1, q) * weight1 + field->get(source2, q) * weight2; } @@ -152,6 +138,11 @@ class InterpolateAndShiftAtBoundary { } } + FloatType modulo(FloatType a, FloatType b) const { + auto const res = std::fmod(a, b); + return (res < FloatType{0}) ? res + b : res; + } + private: std::shared_ptr m_blocks; BlockDataID m_field_id; diff --git a/testsuite/python/lb_lees_edwards_particle_coupling.py b/testsuite/python/lb_lees_edwards_particle_coupling.py index bb212cdc93..713a43e5fe 100644 --- a/testsuite/python/lb_lees_edwards_particle_coupling.py +++ b/testsuite/python/lb_lees_edwards_particle_coupling.py @@ -54,28 +54,14 @@ def coupling_weight(pos_lb_units, node_idx, lb_shape): # minimum image convention, node and coupling position can be # at different sides of a periodic boundary dx = np.abs(min_image_dist(pos_lb_units, node_idx, lb_shape)) - # If the coupling point is >=1 lattice constant away from the node, no coupling. - if np.any(dx >= 1): - weight = 0 - else: - # distance pos to node via module with lattice constant 1 - weight = np.product(1 - dx) + # If the coupling point is >=1 lattice constant away from the node, + # no coupling. Otherwise, distance pos to node via module with lattice + # constant 1 + weight = 0. if np.any(dx >= 1.) else np.product(1. - dx) return weight -class MockLBF: - shape = np.array((10, 10, 10)) - agrid = 1 - - def __getitem__(self, idx): - assert within_grid(idx, self.shape) - return f"node {idx}" - - -mock_lbf = MockLBF() - - def le_aware_lb_nodes_around_pos( folded_pos, lbf, le_pos_offset, shear_direction, shear_plane_normal): """Returns LB node(s) relevant for interpolation around the given position""" @@ -156,13 +142,16 @@ def lb_node(idx): return lbf[fold_index(idx, lbf.shape)] @utx.skipIfMissingFeatures("WALBERLA") +@ut.skipIf(np.prod(system.cell_system.node_grid) != 1, "Requires 1 MPI rank") class LBLeesEdwardsParticleCoupling(ut.TestCase): + """Test LB Lees-Edwards corner cases with a random RNG seed (smoke test)""" + def test_viscous_coupling_with_offset(self): system.lb = None system.time_step = 1 system.cell_system.skin = 0.1 system.cell_system.set_n_square() - offset = 2 * (np.random.random() - 1) * 3 * system.box_l[1] + offset = (np.random.random() - 1.) * 6. * system.box_l[1] protocol = lees_edwards.LinearShear( shear_velocity=0, initial_pos_offset=offset, time_0=0.) system.lees_edwards.set_boundary_conditions( @@ -244,10 +233,10 @@ def v_x(x): return np.interp( if abs(y <= 0.5): pref = -1. dist_to_unshifted_lb_nodes = 0.5 - y - elif y >= system.box_l[1] - 0.5: + else: + assert y >= system.box_l[1] - 0.5 pref = 1. dist_to_unshifted_lb_nodes = y - (system.box_l[2] - 0.5) - else: raise Exception() vel_shift = pref * shear_vel xs = 0.5 + np.arange(lbf.shape[0]) ys = [v_x(x - pref * pos_offset) for x in xs]