Skip to content

Commit

Permalink
Formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
@RudolfWeeber
RudolfWeeber committed Jul 4, 2024
1 parent b95fd23 commit 156d195
Show file tree
Hide file tree
Showing 3 changed files with 61 additions and 49 deletions.
7 changes: 5 additions & 2 deletions src/core/lb/particle_coupling.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -277,8 +277,11 @@ void ParticleCoupling::kernel(std::vector<Particle *> const &particles) {
#ifndef LEES_EDWARDS
auto const drag_force = lb_drag_force(p, m_thermostat.gamma, v_fluid);
#else
auto vel_correction = lees_edvards_vel_shift(m_box_geo.folded_position(p.pos()),*it_positions_force_coupling,m_box_geo);
auto const drag_force = lb_drag_force(p, m_thermostat.gamma, v_fluid-vel_correction);
auto vel_correction =
lees_edvards_vel_shift(m_box_geo.folded_position(p.pos()),
*it_positions_force_coupling, m_box_geo);
auto const drag_force =
lb_drag_force(p, m_thermostat.gamma, v_fluid - vel_correction);
#endif
auto const random_force = get_noise_term(p);
force_on_particle = drag_force + random_force;
Expand Down
5 changes: 3 additions & 2 deletions src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,8 +131,9 @@ class InterpolateAndShiftAtBoundary {
auto folded_source_pos = python_modulo(source_pos, length);
// 0<=folded_source_pos <length
source1[dir] = cell_idx_c(std::floor(folded_source_pos));
// 0<=source1[dir]<length, i.e. integer value sbetw. 0 and length-1 inclusive
source2[dir] = python_modulo(source1[dir]+1, length);
// 0<=source1[dir]<length, i.e. integer value sbetw. 0 and length-1
// inclusive
source2[dir] = python_modulo(source1[dir] + 1, length);
// ineger values between 0 and length -1 inclusive

for (uint_t q = 0; q < FieldType::F_SIZE; ++q) {
Expand Down
98 changes: 53 additions & 45 deletions testsuite/python/lb_lees_edwards_particle_coupling.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,25 +38,27 @@ def unit_vec(k):
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 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 posiiton is
# 2. To determine the lb node to position distance, we need
# minimum image convetion, node and coupling position can be at different sides
# of a periodic boudnary
dx = np.abs(min_image_dist(pos_lb_units,node_idx, lb_shape))
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
if np.any(dx >= 1):
weight = 0
else:
# distance pos to node via module with lattice constant 1
weight =np.product(1-dx)
# distance pos to node via module with lattice constant 1
weight = np.product(1 - dx)

return weight

Expand Down Expand Up @@ -164,7 +166,6 @@ def lb_node(idx): return lbf[fold_index(idx, lbf.shape)]
return unshifted_nodes, shifted_nodes, unshifted_weights, shifted_weights



@utx.skipIfMissingFeatures("WALBERLA")
class LBLeesEdwardsParticleCoupling(ut.TestCase):
def test_viscous_coupling_with_offset(self):
Expand Down Expand Up @@ -239,46 +240,51 @@ def check_velocity_interpolation(self, pos_offset, shear_vel, test_positions):
system.lb = lbf
system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1)
system.part.clear()
v_x = lambda x: 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,:])

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]
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
pref = -1
dist_to_unshifted_lb_nodes = 0.5 - y
elif y >= system.box_l[1] - .5:
pref = 1
dist_to_unshifted_lb_nodes = y-(system.box_l[2]-.5)
pref = 1
dist_to_unshifted_lb_nodes = y - (system.box_l[2] - .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]
v_x_shifted = lambda x: np.interp(x,xs,ys,period=system.box_l[0])
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])
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, 0, pos_all)
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, 0, pos_all)

def test_viscous_coupling_with_shear_vel(self):
# Places a co-moving particle close to the LE boundary
Expand All @@ -302,12 +308,12 @@ def test_viscous_coupling_with_shear_vel(self):
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)
for n in lbf[:, :, :]:
np.testing.assert_allclose(n.velocity[1:], [0, 0], atol=1E-8)
pos = np.random.random(3) * system.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)
np.testing.assert_allclose(p.v[1:], [0, 0], atol=1E-8)
for _ in range(1000):
system.integrator.run(1)
np.testing.assert_allclose(np.copy(p.f), np.zeros(3), atol=2E-6)
Expand Down Expand Up @@ -335,9 +341,11 @@ def test_momentum_conservation(self):
initial_mom = np.copy(system.analysis.linear_momentum())
for i 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)
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(initial_mom[1:], current_mom[1:], atol=2E-7)
np.testing.assert_allclose(
initial_mom[1:], current_mom[1:], atol=2E-7)


if __name__ == '__main__':
Expand Down

0 comments on commit 156d195

Please sign in to comment.