Skip to content
Permalink
Browse files

Fix LBGPU velocity interpolation close to boundary (#3593)

Fixes #

Description of changes:
LB GPu velocity interpolation used incorrect velocities at boundary nodes. This lead to incorrect resutls for interpolations closer than one lattice constant to the boundary.
This adds a missing unit conversion 1/lattice_speed for boundary nodes, which fixes the issue.

Tests for interpolation at the boundary position and at the mid-point between the boundary and the next lattice site.

This should go into 4.1.3
  • Loading branch information
kodiakhq committed Mar 26, 2020
2 parents bf30a6c + ecaedad commit 0479f39ebab008a068b0738c8a5df504ae3040b1
@@ -1372,8 +1372,10 @@ __device__ __inline__ float3 node_velocity(float rho_eq, LB_nodes_gpu n_a,
auto const boundary_index = n_a.boundary[index];

if (boundary_index) {
auto const inv_lattice_speed = para->tau / para->agrid;
auto const &u = n_a.boundary_velocity[index];
return make_float3(u[0], u[1], u[2]);
return make_float3(inv_lattice_speed * u[0], inv_lattice_speed * u[1],
inv_lattice_speed * u[2]);
}

auto const rho = rho_eq + calc_mode_x_from_n(n_a, index, 0);
@@ -384,7 +384,8 @@ def __init__(self, system, **kwargs):
self.highlighted_particle = {}
self.particle_attributes = []
for d in dir(ParticleHandle):
if type(getattr(ParticleHandle, d)) == type(ParticleHandle.pos):
if isinstance(getattr(ParticleHandle, d),
type(ParticleHandle.pos)):
if d not in ["pos_folded"]:
self.particle_attributes.append(d)
self.max_len_attr = max([len(a) for a in self.particle_attributes])
@@ -1371,14 +1372,14 @@ def display():

# pylint: disable=unused-argument
def keyboard_up(button, x, y):
if type(button) is bytes:
if isinstance(button, bytes):
button = button.decode("utf-8")
self.keyboard_manager.keyboard_up(button)
return

# pylint: disable=unused-argument
def keyboard_down(button, x, y):
if type(button) is bytes:
if isinstance(button, bytes):
button = button.decode("utf-8")
self.keyboard_manager.keyboard_down(button)
return
@@ -18,6 +18,7 @@
import unittest_decorators as utx
import numpy as np
import itertools
import sys

import espressomd
import espressomd.shapes
@@ -57,9 +58,9 @@ class LBInterpolation:
def set_boundaries(self, velocity):
"""Place boundaries *not* exactly on a LB node."""
wall_shape1 = espressomd.shapes.Wall(
normal=[1, 0, 0], dist=0.6 * AGRID)
normal=[1, 0, 0], dist=AGRID)
wall_shape2 = espressomd.shapes.Wall(
normal=[-1, 0, 0], dist=-(BOX_L - 0.6 * AGRID))
normal=[-1, 0, 0], dist=-(BOX_L - AGRID))
self.system.lbboundaries.add(
espressomd.lbboundaries.LBBoundary(shape=wall_shape1))
self.system.lbboundaries.add(
@@ -72,9 +73,25 @@ def test_interpolated_velocity(self):
"""
self.set_boundaries([0.0, 0.0, V_BOUNDARY])
self.system.integrator.run(250)
# Shear plane for boundary 1
# for pos in itertools.product((AGRID,), np.arange(0.5 * AGRID, BOX_L, AGRID), np.arange(0.5 * AGRID, BOX_L, AGRID)):
# np.testing.assert_almost_equal(self.lbf.get_interpolated_velocity(pos)[2], 0.0)
# Check interpolated vel at upper boundary. The node position is at
# box_l[0]-agrid/2.
np.testing.assert_allclose(
np.copy(self.lbf.get_interpolated_velocity(
[self.system.box_l[0] - AGRID / 2, 0, 0])),
np.array([0, 0, V_BOUNDARY]))

# Check interpolated velocity involving boundary and neighboring node
# The boundary node index is lbf.shape[0]-1, so -2 refers to the
# node in front of the boundary.
node_next_to_boundary = self.lbf[
self.lbf.shape[0] - 2, 0, 0]
# The midpoint between the boundary and
# that node is box_l - agrid.
np.testing.assert_allclose(
np.copy(self.lbf.get_interpolated_velocity(
[self.system.box_l[0] - AGRID, 0, 0])),
0.5 * (np.array([0, 0, V_BOUNDARY]) + node_next_to_boundary.velocity))

# Bulk
for pos in itertools.product(
np.arange(1.5 * AGRID, BOX_L - 1.5 * AGRID, 0.5 * AGRID),
@@ -93,9 +110,15 @@ def test_mach_limit_check(self):
"""
max_vel = 0.31 * AGRID / TAU
print("Begin: Test error generation")
sys.stdout.flush()
sys.stderr.flush()
with self.assertRaises(Exception):
self.set_boundaries([0.0, 0.0, max_vel])
self.system.integrator.run(1)
sys.stdout.flush()
sys.stderr.flush()
print("End: Test error generation")


@utx.skipIfMissingFeatures(['LB_BOUNDARIES'])

0 comments on commit 0479f39

Please sign in to comment.
You can’t perform that action at this time.