Skip to content

Commit

Permalink
Merge pull request #1 from jngrad/fix-4143
Browse files Browse the repository at this point in the history 
Fix 4143
  • Loading branch information
@KaiserMartin
KaiserMartin committed Mar 31, 2021
2 parents 15aa38d + 5ceb6f9 commit 13e04c5
Show file tree
Hide file tree
Showing 3 changed files with 123 additions and 132 deletions.
201 changes: 87 additions & 114 deletions src/python/espressomd/lb.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -549,136 +549,109 @@ cdef class LBFluidRoutines:
def __set__(self, value):
raise NotImplementedError

cdef class LBSlice:
cdef np.ndarray x_indices
cdef np.ndarray y_indices
cdef np.ndarray z_indices

class LBSlice:

def __init__(self, key):
shape = lb_lbfluid_get_shape()
self.x_indices, self.y_indices, self.z_indices = self.get_indices(
key, shape[0], shape[1], shape[2])

def get_indices(self, key, shape_x, shape_y, shape_z):
_x_indices = np.atleast_1d(np.arange(shape_x)[key[0]])
_y_indices = np.atleast_1d(np.arange(shape_y)[key[1]])
_z_indices = np.atleast_1d(np.arange(shape_z)[key[2]])
return _x_indices, _y_indices, _z_indices
x_indices = np.atleast_1d(np.arange(shape_x)[key[0]])
y_indices = np.atleast_1d(np.arange(shape_y)[key[1]])
z_indices = np.atleast_1d(np.arange(shape_z)[key[2]])
return x_indices, y_indices, z_indices

def get_values(self, _x_indices, _y_indices,
_z_indices, prop_name, shape_res):
def get_values(self, x_indices, y_indices,
z_indices, prop_name, shape_res):
res = np.zeros(
(_x_indices.size,
_y_indices.size,
_z_indices.size,
(x_indices.size,
y_indices.size,
z_indices.size,
*shape_res))
for i, x in enumerate(_x_indices):
for j, y in enumerate(_y_indices):
for k, z in enumerate(_z_indices):
for i, x in enumerate(x_indices):
for j, y in enumerate(y_indices):
for k, z in enumerate(z_indices):
res[i, j, k] = getattr(LBFluidRoutines(
np.array([x, y, z])), prop_name)
return res

def set_values(self, _x_indices, _y_indices, _z_indices, prop_name, value):
for i, x in enumerate(_x_indices):
for j, y in enumerate(_y_indices):
for k, z in enumerate(_z_indices):
def set_values(self, x_indices, y_indices, z_indices, prop_name, value):
for i, x in enumerate(x_indices):
for j, y in enumerate(y_indices):
for k, z in enumerate(z_indices):
setattr(LBFluidRoutines(
np.array([x, y, z])), prop_name, value[i, j, k])

property density:
def __get__(self):
prop_name = "density"
shape_res = (1,)
return np.squeeze(self.get_values(
self.x_indices, self.y_indices, self.z_indices, prop_name, shape_res), axis=3)

def __set__(self, value):
prop_name = "density"
self.set_values(
self.x_indices,
self.y_indices,
self.z_indices,
prop_name,
value)

property index:
def __get__(self):
prop_name = "index"
shape_res = (3,)
return self.get_values(
self.x_indices, self.y_indices, self.z_indices, prop_name, shape_res)

property velocity:
def __get__(self):
prop_name = "velocity"
shape_res = (3,)
return self.get_values(
self.x_indices, self.y_indices, self.z_indices, prop_name, shape_res)

def __set__(self, value):
prop_name = "velocity"
if value.shape == (len(self.x_indices), len(
self.y_indices), len(self.z_indices), 3):
self.set_values(
self.x_indices,
self.y_indices,
self.z_indices,
prop_name,
value)
else:
raise ValueError(
"Input-dimensions of velocity array", value.shape, "does not match slice dimensions",
(len(self.x_indices), len(self.y_indices), len(self.z_indices), 3), ".")

property pressure_tensor:
def __get__(self):
prop_name = "pressure_tensor"
shape_res = (3, 3)
return self.get_values(
self.x_indices, self.y_indices, self.z_indices, prop_name, shape_res)

def __set__(self, value):
raise NotImplementedError

property pressure_tensor_neq:
def __get__(self):
prop_name = "pressure_tensor_neq"
shape_res = (3, 3)
return self.get_values(
self.x_indices, self.y_indices, self.z_indices, prop_name, shape_res)

def __set__(self, value):
raise NotImplementedError

property population:
def __get__(self):
prop_name = "population"
shape_res = (19,)
return self.get_values(
self.x_indices, self.y_indices, self.z_indices, prop_name, shape_res)

def __set__(self, value):
prop_name = "population"
if value.shape == (len(self.x_indices), len(
self.y_indices), len(self.z_indices), 19):
self.set_values(
self.x_indices,
self.y_indices,
self.z_indices,
prop_name,
value)
else:
raise ValueError(
"Input-dimensions of population array", value.shape, "does not match slice dimensions",
(len(self.x_indices), len(self.y_indices), len(self.z_indices), 19), ".")
@property
def density(self):
return np.squeeze(self.get_values(
self.x_indices, self.y_indices, self.z_indices, "density", (1,)), axis=3)

@density.setter
def density(self, value):
self.set_values(self.x_indices, self.y_indices, self.z_indices,
"density", value)

@property
def index(self):
return self.get_values(
self.x_indices, self.y_indices, self.z_indices, "index", (3,))

@property
def velocity(self):
return self.get_values(
self.x_indices, self.y_indices, self.z_indices, "velocity", (3,))

@velocity.setter
def velocity(self, value):
s = (len(self.x_indices), len(self.y_indices), len(self.z_indices), 3)
if value.shape == s:
self.set_values(self.x_indices, self.y_indices, self.z_indices,
"velocity", value)
else:
raise ValueError(
f"Input-dimensions of velocity array {value.shape} does not match slice dimensions {s}.")

@property
def pressure_tensor(self):
return self.get_values(
self.x_indices, self.y_indices, self.z_indices, "pressure_tensor", (3, 3))

@pressure_tensor.setter
def pressure_tensor(self, value):
raise NotImplementedError

@property
def pressure_tensor_neq(self):
return self.get_values(
self.x_indices, self.y_indices, self.z_indices, "pressure_tensor_neq", (3, 3))

@pressure_tensor_neq.setter
def pressure_tensor_neq(self, value):
raise NotImplementedError

@property
def population(self):
return self.get_values(
self.x_indices, self.y_indices, self.z_indices, "population", (19,))

@population.setter
def population(self, value):
s = (len(self.x_indices), len(self.y_indices), len(self.z_indices), 19)
if value.shape == s:
self.set_values(self.x_indices, self.y_indices, self.z_indices,
"population", value)
else:
raise ValueError(
f"Input-dimensions of population array {value.shape} does not match slice dimensions {s}.")

property boundary:
def __get__(self):
shape_res = (1,)
prop_name = "boundary"
return np.squeeze(self.get_values(
self.x_indices, self.y_indices, self.z_indices, prop_name, shape_res), axis=3)
@property
def boundary(self):
return np.squeeze(self.get_values(
self.x_indices, self.y_indices, self.z_indices, "boundary", (1,)), axis=3)

def __set__(self, value):
raise NotImplementedError
@boundary.setter
def boundary(self, value):
raise NotImplementedError
2 changes: 1 addition & 1 deletion testsuite/python/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,7 +187,7 @@ python_test(FILE actor.py MAX_NUM_PROC 1)
python_test(FILE drude.py MAX_NUM_PROC 2)
python_test(FILE thermalized_bond.py MAX_NUM_PROC 4)
python_test(FILE thole.py MAX_NUM_PROC 4)
python_test(FILE test_lb_slice.py MAX_NUM_PROC 1)
python_test(FILE lb_slice.py MAX_NUM_PROC 1)
python_test(FILE lb_switch.py MAX_NUM_PROC 1 LABELS gpu)
python_test(FILE lb_boundary_velocity.py MAX_NUM_PROC 1)
python_test(FILE lb_boundary_volume_force.py MAX_NUM_PROC 4)
Expand Down
52 changes: 35 additions & 17 deletions testsuite/python/test_lb_slice.py → testsuite/python/lb_slice.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,11 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
import espressomd.lb
import espressomd.lbboundaries
import espressomd.shapes
import unittest as ut
import unittest_decorators as utx
import numpy as np


@utx.skipIfMissingFeatures(["LB_BOUNDARIES"])
class TestLBSliceTest(ut.TestCase):
class LBSliceTest(ut.TestCase):

"""This simple test first writes random numbers and then reads them
to same slices of LB nodes and compares if the results are the same,
Expand All @@ -33,56 +29,78 @@ class TestLBSliceTest(ut.TestCase):
system = espressomd.System(box_l=[10.0, 10.0, 10.0])
system.time_step = .01
system.cell_system.skin = 0.1
np.random.seed(seed=42)

def test_LBSlice(self):
def test_slicing(self):
system = self.system

lb_fluid = espressomd.lb.LBFluid(
agrid=1.0, dens=1., visc=1., tau=0.01)
system.actors.add(lb_fluid)

"velocity on test slice [:-1, :-1, -1]"
# velocity on test slice [:-1, :-1, -1]
input_vel = np.random.rand(9, 9, 9, 3)
lb_fluid[:-1, :-1, :-1].velocity = input_vel
output_vel = lb_fluid[:-1, :-1, :-1].velocity
np.testing.assert_array_almost_equal(input_vel, output_vel)

"density on test slice [1:-1:2, 5, 3:6:2]"
with self.assertRaisesRegex(ValueError, r"Input-dimensions of velocity array \(9, 9, 9, 2\) does not match slice dimensions \(9, 9, 9, 3\)"):
lb_fluid[:-1, :-1, :-1].velocity = input_vel[:, :, :, :2]

# density on test slice [1:-1:2, 5, 3:6:2]
input_dens = np.random.rand(4, 1, 2)
lb_fluid[1:-1:2, 5, 3:6:2].density = input_dens
output_dens = lb_fluid[1:-1:2, 5, 3:6:2].density
np.testing.assert_array_almost_equal(input_dens, output_dens)

"population on test slice [:, :, :]"
# population on test slice [:, :, :]
input_pop = np.random.rand(10, 10, 10, 19)
lb_fluid[:, :, :].population = input_pop
output_pop = lb_fluid[:, :, :].population
np.testing.assert_array_almost_equal(input_pop, output_pop)

"pressure tensor on test slice [3, 6, 2:5], should be of shape (1, 1, 3, 3, 3)"
with self.assertRaisesRegex(ValueError, r"Input-dimensions of population array \(10, 10, 10, 5\) does not match slice dimensions \(10, 10, 10, 19\)"):
lb_fluid[:, :, :].population = input_pop[:, :, :, :5]

# pressure tensor on test slice [3, 6, 2:5], should be of shape
# (1, 1, 3, 3, 3)
output_pressure_shape = lb_fluid[3, 6, 2:5].pressure_tensor.shape
should_pressure_shape = (1, 1, 3, 3, 3)
np.testing.assert_array_almost_equal(
output_pressure_shape, should_pressure_shape)

"pressure tensor neq on test slice [3, 6, 2:10], should be of shape (2, 1, 8, 3, 3)"
with self.assertRaises(NotImplementedError):
lb_fluid[3, 6, 2:5].pressure_tensor = [1, 2]

# pressure tensor neq on test slice [3, 6, 2:10], should be of shape
# (2, 1, 8, 3, 3)
output_pressure_neq_shape = lb_fluid[3:5,
6:7, 2:10].pressure_tensor_neq.shape
should_pressure_neq_shape = (2, 1, 8, 3, 3)
np.testing.assert_array_almost_equal(
output_pressure_neq_shape, should_pressure_neq_shape)

"index on test slice [1, 1:5, 6:], should be of shape (1, 4, 4, 3)"
with self.assertRaises(NotImplementedError):
lb_fluid[3, 6, 2:5].pressure_tensor_neq = [1, 2]

# index on test slice [1, 1:5, 6:], should be of shape (1, 4, 4, 3)
output_index_shape = lb_fluid[1, 1:5, 6:].index.shape
should_index_shape = (1, 4, 4, 3)
np.testing.assert_array_almost_equal(
output_index_shape, should_index_shape)

"boundary on test slice [1:, 1:, 1:], should be of shape (9, 9, 9)"
output_boundary_shape = lb_fluid[1:, 1:, 1:].boundary.shape
should_boundary_shape = (9, 9, 9)
np.testing.assert_array_almost_equal(
output_boundary_shape, should_boundary_shape)
with self.assertRaisesRegex(AttributeError, "can't set attribute"):
lb_fluid[1, 1:5, 6:].index = [1, 2]

# boundary on test slice [1:, 1:, 1:], should be of shape (9, 9, 9)
if espressomd.has_features('LB_BOUNDARIES'):
output_boundary_shape = lb_fluid[1:, 1:, 1:].boundary.shape
should_boundary_shape = (9, 9, 9)
np.testing.assert_array_almost_equal(
output_boundary_shape, should_boundary_shape)

with self.assertRaises(NotImplementedError):
lb_fluid[3, 6, 2:5].pressure_tensor_neq = [1, 2]


if __name__ == "__main__":
Expand Down

0 comments on commit 13e04c5

Please sign in to comment.