Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix 4143 #1

Merged
merged 4 commits into from
Mar 31, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
201 changes: 87 additions & 114 deletions src/python/espressomd/lb.pyx
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
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
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