diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index a9de5d3039..b7b3eacd47 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.cpp @@ -18,6 +18,7 @@ */ #include "LocalBox.hpp" #include "Particle.hpp" +#include "algorithm/periodic_fold.hpp" #include "cells.hpp" #include "communication.hpp" #include "config/config.hpp" @@ -162,27 +163,34 @@ bool in_local_halo(Utils::Vector3d const &pos) { * @brief Return a vector of positions shifted by +,- box length in each * coordinate */ -std::vector positions_in_halo(Utils::Vector3d pos, - const BoxGeometry &box) { - std::vector res; +std::vector> +positions_in_halo(Utils::Vector3d pos, const BoxGeometry &box) { + std::vector> res; for (int i : {-1, 0, 1}) { for (int j : {-1, 0, 1}) { for (int k : {-1, 0, 1}) { Utils::Vector3d shift{{double(i), double(j), double(k)}}; Utils::Vector3d pos_shifted = pos + Utils::hadamard_product(box.length(), shift); - + 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 - pos)[le.shear_plane_normal]; - if (normal_shift > 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 = Algorithm::periodic_fold( + le.pos_offset, box_geo.length()[le.shear_direction]); + if (normal_shift > std::numeric_limits::epsilon()) { + pos_shifted[le.shear_direction] += folded_offset; + vel_shift[le.shear_direction] -= le.shear_velocity; + } + if (normal_shift < -std::numeric_limits::epsilon()) { + pos_shifted[le.shear_direction] -= folded_offset; + vel_shift[le.shear_direction] += le.shear_velocity; + } } if (in_local_halo(pos_shifted)) { - res.push_back(pos_shifted); + res.push_back({pos_shifted, vel_shift}); + // std::cout << "coupling at "< 0.0, p.id(), rng_counter); @@ -274,7 +286,8 @@ void couple_particle(Particle &p, bool couple_virtual, double noise_amplitude, // couple positions including shifts by one box length to add // forces to ghost layers - for (auto pos : positions_in_halo(p.pos(), box_geo)) { + for (auto shift : positions_in_halo(p.pos(), box_geo)) { + auto const pos = shift.first; if (in_local_domain(pos)) { /* Particle is in our LB volume, so this node * is responsible to adding its force */ diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.hpp b/src/core/grid_based_algorithms/lb_particle_coupling.hpp index c93a98ca00..2357cdba6f 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.hpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.hpp @@ -96,8 +96,8 @@ Utils::Vector3d lb_particle_coupling_noise(bool enabled, int part_id, const OptionalCounter &rng_counter); // internal function exposed for unit testing -std::vector positions_in_halo(Utils::Vector3d pos, - const BoxGeometry &box); +std::vector> +positions_in_halo(Utils::Vector3d pos, const BoxGeometry &box); // internal function exposed for unit testing void couple_particle(Particle &p, bool couple_virtual, double noise_amplitude, diff --git a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp b/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp index 6c06d9d665..af0f5f4105 100644 --- a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp +++ b/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp @@ -59,7 +59,8 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) { return; } if (should_be_coupled(p, coupled_ghost_particles)) { - for (auto pos : positions_in_halo(p.pos(), box_geo)) { + for (auto shift : positions_in_halo(p.pos(), box_geo)) { + auto const pos = shift.first; add_md_force(pos * to_lb_units, -p.force(), time_step); } } @@ -71,7 +72,8 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) { return; } if (should_be_coupled(p, coupled_ghost_particles)) { - for (auto pos : positions_in_halo(p.pos(), box_geo)) { + for (auto shift : positions_in_halo(p.pos(), box_geo)) { + auto const pos = shift.first; add_md_force(pos * to_lb_units, -p.force(), time_step); } } diff --git a/src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp b/src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp index 72b5c32697..4466dcc724 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 { +double python_modulo(double a, double b) { + double 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,35 @@ 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); + auto const weight1 = 1 - std::fmod(folded_offset, FloatType{1}); + auto const weight2 = std::fmod(folded_offset, FloatType{1}); for (auto const &&cell : ci) { Cell source1 = cell; Cell source2 = cell; - source1[dir] = cell_idx_c(std::floor( - static_cast(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]; + auto ghost_fold = [dim](FloatType x) { + // std::cout << "ghost fold " << x<< " "; + while (x > dim) { + x -= dim; + }; + while (x < -1) { + x += dim; + }; + // std::cout <(source2[dir]) + offset)) % - dim; - source2[dir] = cell_idx_c(static_cast(source2[dir]) + length); - source2[dir] = cell_idx_c(source2[dir] % dim); - - 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; + cell_idx_c(std::floor(ghost_fold(target_idx + folded_offset + 1))); + // std::cout << offset<<"->"<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 6c230b1585..4ed9d71e54 100644 --- a/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp +++ b/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp @@ -33,6 +33,7 @@ #include #include #include +#include #include #include diff --git a/testsuite/python/lb_lees_edwards_particle_coupling.py b/testsuite/python/lb_lees_edwards_particle_coupling.py index 373d2eca25..91c48f0684 100644 --- a/testsuite/python/lb_lees_edwards_particle_coupling.py +++ b/testsuite/python/lb_lees_edwards_particle_coupling.py @@ -23,18 +23,18 @@ import numpy as np import unittest_decorators as utx +system = espressomd.System(box_l=[10, 10, 10]) + @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.actors.clear() system.time_step = 1 system.cell_system.skin = 0.1 system.cell_system.set_n_square() - offset = 1 - idx = int(offset) + offset = 0.9 # (np.random.random()-1/2) * 5*system.box_l[0] protocol = lees_edwards.LinearShear( shear_velocity=0., initial_pos_offset=offset, time_0=0.) system.lees_edwards.set_boundary_conditions( @@ -45,49 +45,119 @@ def test(self): system.actors.add(lbf) system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1) + idx = int(offset % system.box_l[0]) # lb x index incl offset + # the particle is placed in the center of 8 lb points, i.e., + # 0.5 away from the planes with the grid points 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_x = lbf.shape[0] // 2 # lb index for the origin particle 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]] + upper_y = lbf.shape[1] - 1 # y is shear plane noremal + # LB Nodes surrounding the particle. + # On the particle's side + unshifted = [lbf[mid_x - 1, 0, mid_z - 1], + lbf[mid_x, 0, mid_z - 1], + lbf[mid_x - 1, 0, mid_z], + lbf[mid_x, 0, mid_z]] + # across the Lees Edwads boundary conditoin + # if offset ^1<.5, it couples to the left neighbor + # otherwise to the right + if (offset % 1) >= .5: extra = 1 + else: extra = 0 + shifted_left = [lbf[(mid_x - 1 + idx + extra) % lbf.shape[0], upper_y, mid_z - 1], + lbf[(mid_x - 1 + idx + extra) % lbf.shape[0], upper_y, mid_z]] + shifted_right = [lbf[(mid_x + idx + extra) % lbf.shape[0], upper_y, mid_z - 1], + lbf[(mid_x + idx + extra) % lbf.shape[0], upper_y, mid_z]] + nodes = shifted_left + shifted_right + unshifted + + v0 = np.array([1, 2, 3]) for n in nodes: n.velocity = v0 system.integrator.run(1) + + # Gather forces applied to the LB by the particle coupling lb_forces = np.array([n.last_applied_force for n in nodes]) lb_force = np.sum(lb_forces, axis=0) + + # total force on lb = - force on particle? np.testing.assert_allclose(lb_force, -np.copy(p.f)) - for f in lb_forces: - np.testing.assert_allclose(f, lb_forces[0]) + # Particle couples to 8 nodes. On the sie of the particle + # each lb node should have 1/8 of the force + for n in unshifted: + np.testing.assert_allclose( + np.copy(n.last_applied_force), -np.copy(p.f) / 8) + + # Across the lees edwards boundary, forces ahead and behind + # the particle in shear directiion are differently distributed + # For offset 0 (in the middle between two nodes) + # left and right weighs are equal (0.5) + weight_r = (offset + .5 - idx) % 1 + weight_l = 1 - weight_r + np.testing.assert_allclose(weight_l + weight_r, 1) + for w in (weight_l, weight_r): + assert w >= 0 and w <= 1 + # The total weight is the product of the weights + # in all 3 Cartesian directions. These are + # weight_l/r in the shear direction and 1/2 in the other two + for n in shifted_left: + np.testing.assert_allclose( + np.copy(n.last_applied_force), -np.copy(p.f) * 1 / 2 * 1 / 2 * weight_l) + for n in shifted_right: + np.testing.assert_allclose( + np.copy(n.last_applied_force), -np.copy(p.f) * 1 / 2 * 1 / 2 * weight_r) + # Check the LB velocity interpolation + # at a positoin, where the ghost layer of the + # lees edwards shear plane contributes + lbf[:, :, :].velocity = np.zeros(3) lbf[:, :, :].velocity = [0, 0, 0] - lower_nodes = nodes[:4] - upper_nodes = nodes[4:] + lower_nodes = unshifted + upper_nodes = shifted_left + shifted_right for n in lower_nodes: n.velocity = v0 for n in upper_nodes: n.velocity = - v0 + # When the offset modulo 1 is not zero + # in the ghost layer, neighboring cell swith zero velocity contribute + expected_v = 1 / 2 * v0 +\ + 1 / 4 * -v0 + 1 / 4 * (1 - abs(1 - (offset % 1))) * -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: - np.testing.assert_allclose( - np.copy(n.last_applied_force), np.zeros(3)) + expected_v) + + def test_viscous_coupling_with_shear_vel(self): + # Places a co-moving particle close to the LE boundary + # in shear flow. chesk that it remains force free + # this is only the case, if the periodic imagesin the + # halo regoin calculate the drag force including the LE + # shear velocity. + system.actors.clear() + system.time_step = 0.1 + system.cell_system.skin = 0.1 + system.cell_system.set_n_square() + v_shear = 0.5 + protocol = lees_edwards.LinearShear( + shear_velocity=v_shear, initial_pos_offset=0, 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.actors.add(lbf) + system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1) + system.integrator.run(200) + p = system.part.add( + pos=( + 0, 0, 0), v=lbf.get_interpolated_velocity( + pos=( + 0, 0, 0))) + for _ in range(100): + system.integrator.run(1) + np.testing.assert_allclose(p.f, np.zeros(3), atol=1E-7) if __name__ == '__main__':