From 4284116825b052c931d02a31a80b47492cb6e0e2 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Mon, 8 Jul 2024 18:06:07 +0200 Subject: [PATCH 1/2] Fix and test LB Lees-Edwards particle coupling --- src/core/lb/particle_coupling.cpp | 45 ++- .../InterpolateAndShiftAtBoundary.hpp | 41 ++- .../src/lattice_boltzmann/LBWalberlaImpl.hpp | 2 +- .../lb_lees_edwards_particle_coupling.py | 340 +++++++++++++++--- 4 files changed, 361 insertions(+), 67 deletions(-) diff --git a/src/core/lb/particle_coupling.cpp b/src/core/lb/particle_coupling.cpp index ae6e62c849..cfe3ce622b 100644 --- a/src/core/lb/particle_coupling.cpp +++ b/src/core/lb/particle_coupling.cpp @@ -123,12 +123,15 @@ 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 std::numeric_limits::epsilon()) - pos_shifted[le.shear_direction] += le.pos_offset; - if (normal_shift < -std::numeric_limits::epsilon()) - pos_shifted[le.shear_direction] -= le.pos_offset; + auto folded_offset = + std::fmod(le.pos_offset, box_geo.length()[le.shear_direction]); + if (folded_offset < 0.) { + folded_offset += box_geo.length()[le.shear_direction]; + } + pos_shifted[le.shear_direction] += + shift[le.shear_plane_normal] * folded_offset; } if (in_box(pos_shifted, halo_lower_corner, halo_upper_corner)) { @@ -163,6 +166,26 @@ std::vector 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) { + Utils::Vector3d vel_shift{{0., 0., 0.}}; + if (box_geo.type() == BoxType::LEES_EDWARDS) { + auto le = box_geo.lees_edwards_bc(); + auto normal_shift = + (pos_shifted_by_box_l - orig_pos)[le.shear_plane_normal]; + // normal_shift is +,- box_l or 0 up to floating point errors + if (normal_shift > std::numeric_limits::epsilon()) { + vel_shift[le.shear_direction] -= le.shear_velocity; + } + if (normal_shift < -std::numeric_limits::epsilon()) { + vel_shift[le.shear_direction] += le.shear_velocity; + } + } + return vel_shift; +} + namespace LB { Utils::Vector3d ParticleCoupling::get_noise_term(Particle const &p) const { @@ -242,6 +265,7 @@ void ParticleCoupling::kernel(std::vector const &particles) { auto const &domain_upper_corner = m_local_box.my_right(); auto it_interpolated_velocities = interpolated_velocities.begin(); auto it_positions_force_coupling = positions_force_coupling.begin(); + auto it_positions_velocity_coupling = positions_velocity_coupling.begin(); auto it_positions_force_coupling_counter = positions_force_coupling_counter.begin(); for (auto ptr : coupled_particles) { @@ -254,11 +278,20 @@ void ParticleCoupling::kernel(std::vector const &particles) { #endif Utils::Vector3d force_on_particle = {}; if (coupling_mode == particle_force) { - auto const &v_fluid = *it_interpolated_velocities; + auto &v_fluid = *it_interpolated_velocities; + if (m_box_geo.type() == BoxType::LEES_EDWARDS) { + // Account for the case where the interpolated velocity has been read + // from a ghost of the particle across the LE boundary (or vice verssa) + // Then the particle velocity is shifted by +,- the LE shear velocity + auto const vel_correction = lees_edwards_vel_shift( + *it_positions_velocity_coupling, p.pos(), m_box_geo); + v_fluid += vel_correction; + } auto const drag_force = lb_drag_force(p, m_thermostat.gamma, v_fluid); auto const random_force = get_noise_term(p); force_on_particle = drag_force + random_force; ++it_interpolated_velocities; + ++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 72b5c32697..d24dd54b46 100644 --- a/src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp +++ b/src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp @@ -31,6 +31,16 @@ #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 { /** @@ -94,8 +104,6 @@ class InterpolateAndShiftAtBoundary { auto const dir = m_shear_direction; auto const dim = cell_idx_c(m_blocks->getNumberOfCells(*block, dir)); auto const length = numeric_cast(dim); - auto const weight = - std::abs(std::fmod(get_pos_offset() + length, FloatType{1})); // setup slab auto field = block->template getData(m_field_id); @@ -111,24 +119,27 @@ class InterpolateAndShiftAtBoundary { 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); + // 0<=folded_offset(source1[dir]) + offset)) % - dim; - source1[dir] = cell_idx_c(static_cast(source1[dir]) + length); - source1[dir] = cell_idx_c(source1[dir] % dim); - + int target_idx = cell[dir]; + FloatType const source_pos = FloatType(target_idx) + folded_offset; + auto folded_source_pos = python_modulo(source_pos, length); + // 0<=folded_source_pos (source2[dir]) + offset)) % - dim; - source2[dir] = cell_idx_c(static_cast(source2[dir]) + length); - source2[dir] = cell_idx_c(source2[dir] % dim); + cell_idx_c(python_modulo(FloatType(source1[dir] + 1), length)); + // ineger values between 0 and length -1 inclusive - for (uint_t f = 0; f < FieldType::F_SIZE; ++f) { - tmp_field->get(cell, f) = field->get(source1, f) * (1 - weight) + - field->get(source2, f) * weight; + for (uint_t q = 0; q < FieldType::F_SIZE; ++q) { + tmp_field->get(cell, q) = + field->get(source1, q) * weight1 + field->get(source2, q) * weight2; } tmp_field->get(cell, m_shear_direction) -= prefactor * shift; } diff --git a/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp b/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp index 98b0679f4b..1f1a8bb24e 100644 --- a/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp +++ b/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp @@ -535,13 +535,13 @@ class LBWalberlaImpl : public LBWalberlaBase { void integrate_pull_scheme() { auto const &blocks = get_lattice().get_blocks(); - integrate_reset_force(blocks); // Handle boundaries integrate_boundaries(blocks); // LB stream integrate_stream(blocks); // LB collide integrate_collide(blocks); + integrate_reset_force(blocks); // Refresh ghost layers ghost_communication_pdfs(); } diff --git a/testsuite/python/lb_lees_edwards_particle_coupling.py b/testsuite/python/lb_lees_edwards_particle_coupling.py index 103197c721..bb212cdc93 100644 --- a/testsuite/python/lb_lees_edwards_particle_coupling.py +++ b/testsuite/python/lb_lees_edwards_particle_coupling.py @@ -21,22 +21,279 @@ import espressomd import espressomd.lb import numpy as np +import itertools +import copy import unittest_decorators as utx +from tests_common import fold_index + +system = espressomd.System(box_l=[10, 10, 10]) + + +def unit_vec(k): + res = np.zeros(3) + res[k] = 1 + return res + + +def within_grid(idx, shape): + return np.all(idx >= 0) and np.all(idx < shape) + + +def min_image_dist(a, b, l): + res = b - a + for i in range(3): + while res[i] < -l[i] / 2: res += l[i] + while res[i] >= l[i] / 2: res -= l[i] + return res + + +def coupling_weight(pos_lb_units, node_idx, lb_shape): + # 1. For the coupling weights it does not matter on which side of the + # lb_node the position is + # 2. To determine the lb node to position distance, we need + # 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) + + 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""" + + # helper to get lb node from index with folding + def lb_node(idx): return lbf[fold_index(idx, lbf.shape)] + + # center of first lb lattice point is at 0.5 in each Cartesian direction. + + # determine position in lb units for 3 cases; + # unshifted, lees-edwards shifted to the left and to the right + pos_unshifted_lb_units = folded_pos / lbf.agrid - 0.5 # relative to node centers + shear_vec = unit_vec(shear_direction) + pos_shifted_left_lb_units = ( + folded_pos - shear_vec * le_pos_offset) / lbf.agrid - 0.5 + pos_shifted_right_lb_units = ( + folded_pos + shear_vec * le_pos_offset) / lbf.agrid - 0.5 + + # Particle couples to its 8 neighboring lattice sites + # when a particle is at a lattice site in any coordinate, + # the right hand side neighbor is included + # but the coupling weight for that neighbor is 0 + + # find the lower left lb node to which the particle couples + lower_idx_unshifted = np.array(np.floor(pos_unshifted_lb_units), dtype=int) + lower_idx_shifted_left = np.array( + np.floor(pos_shifted_left_lb_units), dtype=int) + lower_idx_shifted_right = np.array( + np.floor(pos_shifted_right_lb_units), dtype=int) + + ijks = np.array(list(itertools.product([0, 1], repeat=3))) + indices_unshifted = [lower_idx_unshifted + ijk for ijk in ijks] + + # Nodes with an index within the primary box in shear_plane_normal direction + # do not need Lees-Edwards handling + dont_need_shift = [ + idx for idx in indices_unshifted if within_grid( + idx[shear_plane_normal], + lbf.shape[shear_plane_normal])] + unshifted_nodes = [lb_node(idx) for idx in dont_need_shift] + unshifted_weights = [ + coupling_weight(pos_unshifted_lb_units, idx, lbf.shape) + for idx in dont_need_shift] + + # Handle indices which are not in the primary box in the sheare plane + # normal + to_be_shifted_left = [ + (idx, ijk) for idx, ijk in zip(indices_unshifted, ijks) + if idx[shear_plane_normal] >= lbf.shape[shear_plane_normal]] + to_be_shifted_right = [(idx, ijk) for idx, ijk in zip( + indices_unshifted, ijks) if idx[shear_plane_normal] < 0] + + # replace the index in shear direction + shifted_left = copy.deepcopy(to_be_shifted_left) + for idx, ijk in shifted_left: + idx[shear_direction] = lower_idx_shifted_left[shear_direction] + \ + ijk[shear_direction] + shifted_right = copy.deepcopy(to_be_shifted_right) + for idx, ijk in shifted_right: + idx[shear_direction] = lower_idx_shifted_right[shear_direction] + \ + ijk[shear_direction] + + weights_shifted_left = [ + coupling_weight(pos_shifted_left_lb_units, idx, lbf.shape) + for idx, _ in shifted_left] + weights_shifted_right = [ + coupling_weight(pos_shifted_right_lb_units, idx, lbf.shape) + for idx, _ in shifted_right] + + shifted_nodes = [lb_node(idx) for idx, _ in (shifted_left + shifted_right)] + shifted_weights = weights_shifted_left + weights_shifted_right + + np.testing.assert_allclose(sum(unshifted_weights + shifted_weights), 1) + assert len(shifted_nodes + unshifted_nodes) == 8 + assert len(set(shifted_nodes + unshifted_nodes)) == 8 # no duplicates + + return unshifted_nodes, shifted_nodes, unshifted_weights, shifted_weights @utx.skipIfMissingFeatures("WALBERLA") class LBLeesEdwardsParticleCoupling(ut.TestCase): - def test(self): - system = espressomd.System(box_l=[10, 10, 10]) + 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] + protocol = lees_edwards.LinearShear( + shear_velocity=0, initial_pos_offset=offset, time_0=0.) + system.lees_edwards.set_boundary_conditions( + shear_direction="x", shear_plane_normal="y", protocol=protocol) + lbf = espressomd.lb.LBFluidWalberla( + agrid=1., density=1., kinematic_viscosity=1., tau=system.time_step) + system.lb = lbf + system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1) + for _ in range(10): + system.part.clear() + lbf[:, :, :].velocity = np.zeros(3) + + x = np.random.random() * system.box_l[0] + z = np.random.random() * system.box_l[2] + # within 0.5 of the lees-edwards boundary + y = (np.random.random() - 0.5) % system.box_l[1] + pos = np.array((x, y, z)) + p = system.part.add(pos=pos) + v0 = np.random.random(3) - 1 / 2 + + nodes_unshifted, nodes_shifted, weights_unshifted, weights_shifted = \ + le_aware_lb_nodes_around_pos(pos, lbf, offset, 0, 1) + all_nodes = nodes_unshifted + nodes_shifted + all_weights = weights_unshifted + weights_shifted + + for n in all_nodes: + n.velocity = v0 + + system.integrator.run(1) + # Gather forces applied to the LB by the particle coupling + lb_force = np.sum( + np.array([n.last_applied_force for n in all_nodes]), axis=0) + + # total force on lb = - force on particle? + np.testing.assert_allclose(lb_force, -np.copy(p.f)) + + # validate our assumptions about which lb nodes get a force + # from the coupling. Exactly the nodes listed in `nodes` + # should have received a force during coupling. + lb_nodes_with_force_idx = sorted( + [n.index for n in lbf[:, :, :] if np.any(n.last_applied_force != 0)]) + expected_nodes_idx = sorted( + [n.index for n, w in zip(all_nodes, all_weights) if w > 0]) + np.testing.assert_array_equal( + lb_nodes_with_force_idx, expected_nodes_idx) + + # force on individual nodes + for n, w in zip(all_nodes, all_weights): + np.testing.assert_allclose( + np.copy(n.last_applied_force), -w * np.copy(p.f)) + + def check_velocity_interpolation(self, pos_offset, shear_vel, test_positions): + system.lb = None + system.part.clear() system.time_step = 1 system.cell_system.skin = 0.1 system.cell_system.set_n_square() + system.time = 0 + protocol = lees_edwards.LinearShear( + shear_velocity=shear_vel, initial_pos_offset=pos_offset, time_0=0.) + system.lees_edwards.set_boundary_conditions( + shear_direction="x", shear_plane_normal="y", protocol=protocol) + lbf = espressomd.lb.LBFluidWalberla( + agrid=1., density=1., kinematic_viscosity=1., tau=system.time_step) + system.lb = lbf + system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1) + system.part.clear() + + def v_x(x): return np.interp( + x, [0.5, lbf.shape[0] - .5], [0, lbf.shape[0] - 1], period=lbf.shape[0]) + nodes_at_y_boundary = list( + lbf[:, 0, :]) + list(lbf[:, lbf.shape[1] - 1, :]) + for n in nodes_at_y_boundary: + node_x = 0.5 + n.index[0] + n.velocity = [v_x(node_x), 0, 0] + for pos in test_positions: + y = pos[1] + if abs(y <= 0.5): + pref = -1. + dist_to_unshifted_lb_nodes = 0.5 - y + elif 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] + def v_x_shifted(x): return np.interp( + x, xs, ys, period=system.box_l[0]) + unshifted_vel = v_x(pos[0]) + shifted_vel = v_x_shifted(pos[0]) + vel_shift + weight_unshifted = 1 - dist_to_unshifted_lb_nodes + weight_shifted = 1 - weight_unshifted + expected_vel = np.array( + [weight_unshifted * unshifted_vel + weight_shifted * shifted_vel, 0, 0]) + observed_vel = np.copy(lbf.get_interpolated_velocity(pos=pos)) + np.testing.assert_allclose(observed_vel, expected_vel) + + def test_vel_interpol_all(self): + n = 25 + xs = np.linspace(0, system.box_l[0], n) + y_ls = [0.2] * n + y_us = [system.box_l[1] - .2] * n + zs = np.random.random(n) * system.box_l[2] + pos_lower = np.vstack((xs, y_ls, zs)).T + pos_upper = np.vstack((xs, y_us, zs)).T + pos_all = np.vstack((pos_lower, pos_upper)) + # non-integer offset + pos_offsets = 100 * system.box_l[0] * (np.random.random(10) - .5) + for pos_offset in pos_offsets: + self.check_velocity_interpolation( + pos_offset, 2 * np.random.random() - 1, pos_all) - offset = 1 - idx = int(offset) + def test_viscous_coupling_with_shear_vel(self): + # Place a co-moving particle close to the LE boundary in shear flow. + # Check that it remains force-free. This is only the case, + # if the periodic images in the halo region calculate + # the drag force including the LE shear velocity. + box_l = system.box_l + system.lb = None + system.part.clear() + system.time_step = 0.1 + system.time = 0 + system.cell_system.skin = 0.1 + system.cell_system.set_n_square() + v_shear = 2. * (np.random.random() - 0.5) protocol = lees_edwards.LinearShear( - shear_velocity=0., initial_pos_offset=offset, time_0=0.) + shear_velocity=v_shear, + initial_pos_offset=(np.random.random() - 0.5) * 5. * box_l[0], + time_0=np.random.random()) system.lees_edwards.set_boundary_conditions( shear_direction="x", shear_plane_normal="y", protocol=protocol) @@ -44,50 +301,43 @@ def test(self): agrid=1., density=1., kinematic_viscosity=1., tau=system.time_step) system.lb = lbf system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1) + system.integrator.run(5000) + for n in lbf[:, :, :]: + np.testing.assert_allclose(n.velocity[1:], [0, 0], atol=1E-8) + pos = np.random.random(3) * box_l + p = system.part.add(pos=pos, v=lbf.get_interpolated_velocity(pos=pos)) + np.testing.assert_allclose(p.v[1:], [0, 0], atol=1E-8) + for _ in range(1000): + system.integrator.run(1, reuse_forces=True) + np.testing.assert_allclose(np.copy(p.f), np.zeros(3), atol=2E-6) - pos = [system.box_l[0] / 2., 0., system.box_l[0] / 2.] - p = system.part.add(pos=pos) - v0 = np.array([1, 2, 3]) - mid_x = lbf.shape[0] // 2 - mid_z = lbf.shape[2] // 2 - - upper_y = lbf.shape[1] - 1 - nodes = [lbf[mid_x - 1, 0, mid_z], - lbf[mid_x, 0, mid_z - 1], - lbf[mid_x - 1, 0, mid_z], - lbf[mid_x, 0, mid_z], - lbf[mid_x - 1 + idx, upper_y, mid_z], - lbf[mid_x + idx, upper_y, mid_z - 1], - lbf[mid_x - 1 + idx, upper_y, mid_z], - lbf[mid_x + idx, upper_y, mid_z]] - for n in nodes: - n.velocity = v0 + def test_momentum_conservation(self): + system.lb = None + system.part.clear() + system.time_step = 0.01 + system.cell_system.skin = 0.1 + system.cell_system.set_n_square() + v_shear = np.random.random() - 0.5 + protocol = lees_edwards.LinearShear( + shear_velocity=v_shear, initial_pos_offset=13.7, time_0=0.) + system.lees_edwards.set_boundary_conditions( + shear_direction="x", shear_plane_normal="y", protocol=protocol) + lbf = espressomd.lb.LBFluidWalberla( + agrid=1., density=1., kinematic_viscosity=1., tau=system.time_step) + system.lb = lbf + system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1) + pos = (0, 0, 0) + p = system.part.add(pos=pos, v=(0, 0, 0)) system.integrator.run(1) - lb_forces = np.array([n.last_applied_force for n in nodes]) - lb_force = np.sum(lb_forces, axis=0) - np.testing.assert_allclose(lb_force, -np.copy(p.f)) - for f in lb_forces: - np.testing.assert_allclose(f, lb_forces[0]) - - lbf[:, :, :].velocity = [0, 0, 0] - - lower_nodes = nodes[:4] - upper_nodes = nodes[4:] - for n in lower_nodes: - n.velocity = v0 - for n in upper_nodes: - n.velocity = - v0 - p.update(dict(pos=pos, v=np.zeros(3))) - np.testing.assert_allclose( - np.copy(lbf.get_interpolated_velocity(pos=pos)), - np.zeros(3)) - system.integrator.run(1) - np.testing.assert_allclose(np.copy(p.pos), pos) - np.testing.assert_allclose(np.copy(p.f), np.zeros(3)) - for n in nodes: + initial_mom = np.copy(system.analysis.linear_momentum()) + for _ in range(100): + system.integrator.run(1) + np.testing.assert_allclose(-np.copy(p.f), np.copy( + np.sum(lbf[:, :, :].last_applied_force, axis=(0, 1, 2))), atol=1E-9) + current_mom = np.copy(system.analysis.linear_momentum()) np.testing.assert_allclose( - np.copy(n.last_applied_force), np.zeros(3)) + initial_mom[1:], current_mom[1:], atol=2E-7) if __name__ == '__main__': 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 2/2] 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]