diff --git a/src/core/grid_based_algorithms/lb_interface.cpp b/src/core/grid_based_algorithms/lb_interface.cpp index 6508e501016..860403b0a4e 100644 --- a/src/core/grid_based_algorithms/lb_interface.cpp +++ b/src/core/grid_based_algorithms/lb_interface.cpp @@ -1107,22 +1107,23 @@ void lb_lbfluid_load_checkpoint(const std::string &filename, int binary) { } } -bool lb_lbnode_is_index_valid(const Utils::Vector3i &ind) { - auto within_bounds = [](const Utils::Vector3i &ind, - const Utils::Vector3i &limits) { - return ind < limits && ind >= Utils::Vector3i{}; - }; +Utils::Vector3i lb_lbfluid_get_shape() { if (lattice_switch == ActiveLB::GPU) { #ifdef CUDA - return within_bounds(ind, {static_cast(lbpar_gpu.dim_x), - static_cast(lbpar_gpu.dim_y), - static_cast(lbpar_gpu.dim_z)}); + return {static_cast(lbpar_gpu.dim_x), + static_cast(lbpar_gpu.dim_y), + static_cast(lbpar_gpu.dim_z)}; #endif } if (lattice_switch == ActiveLB::CPU) { - return within_bounds(ind, lblattice.global_grid); + return lblattice.global_grid; } - return false; + throw std::runtime_error("No LB active"); +} + +bool lb_lbnode_is_index_valid(const Utils::Vector3i &ind) { + auto const limit = lb_lbfluid_get_shape(); + return ind < limit && ind >= Utils::Vector3i{}; } double lb_lbnode_get_density(const Utils::Vector3i &ind) { diff --git a/src/core/grid_based_algorithms/lb_interface.hpp b/src/core/grid_based_algorithms/lb_interface.hpp index 4539735c4fb..aec7647b427 100644 --- a/src/core/grid_based_algorithms/lb_interface.hpp +++ b/src/core/grid_based_algorithms/lb_interface.hpp @@ -216,6 +216,11 @@ void lb_lbfluid_load_checkpoint(const std::string &filename, int binary); */ bool lb_lbnode_is_index_valid(const Utils::Vector3i &ind); +/** + * @brief returns the shape of the LB fluid lattice + */ +Utils::Vector3i lb_lbfluid_get_shape(); + void lb_lbfluid_on_lb_params_change(LBParam field); Utils::Vector3d lb_lbfluid_calc_fluid_momentum(); diff --git a/src/python/espressomd/lb.pxd b/src/python/espressomd/lb.pxd index 711dd2f2778..c770baac51a 100644 --- a/src/python/espressomd/lb.pxd +++ b/src/python/espressomd/lb.pxd @@ -76,6 +76,7 @@ cdef extern from "grid_based_algorithms/lb_interface.hpp": ActiveLB lb_lbfluid_get_lattice_switch() except + Vector6d lb_lbfluid_get_stress() except + bool lb_lbnode_is_index_valid(const Vector3i & ind) except + + Vector3i lb_lbfluid_get_shape() except + const Vector3d lb_lbnode_get_velocity(const Vector3i & ind) except + void lb_lbnode_set_velocity(const Vector3i & ind, const Vector3d & u) except + double lb_lbnode_get_density(const Vector3i & ind) except + diff --git a/src/python/espressomd/lb.pyx b/src/python/espressomd/lb.pyx index b06febec7a1..ac785c821d7 100644 --- a/src/python/espressomd/lb.pyx +++ b/src/python/espressomd/lb.pyx @@ -20,6 +20,7 @@ from __future__ import print_function, absolute_import, division include "myconfig.pxi" import os import cython +import itertools import numpy as np cimport numpy as np from libc cimport stdint @@ -244,6 +245,28 @@ cdef class HydrodynamicInteraction(Actor): def _deactivate_method(self): lb_lbfluid_set_lattice_switch(NONE) + + property shape: + def __get__(self): + cdef Vector3i shape = lb_lbfluid_get_shape() + return (shape[0], shape[1], shape[2]) + + property stress: + def __get__(self): + cdef Vector6d res + res = lb_lbfluid_get_stress() + return array_locked(( + res[0], res[1], res[2], res[3], res[4], res[5])) + + def __set__(self, value): + raise NotImplementedError + + def nodes(self): + """Provides a generator for iterating over all lb nodes""" + + shape = self.shape + for i, j, k in itertools.product(range(shape[0]), range(shape[1]), range(shape[2])): + yield self[i, j, k] # LBFluid main class @@ -261,17 +284,7 @@ cdef class LBFluid(HydrodynamicInteraction): self.validate_params() self._set_lattice_switch() self._set_params_in_es_core() - - property stress: - def __get__(self): - cdef Vector6d res - res = lb_lbfluid_get_stress() - return array_locked(( - res[0], res[1], res[2], res[3], res[4], res[5])) - - def __set__(self, value): - raise NotImplementedError - + IF CUDA: cdef class LBFluidGPU(HydrodynamicInteraction): """ @@ -330,6 +343,10 @@ cdef class LBFluidRoutines(object): self.node[2] = key[2] if not lb_lbnode_is_index_valid(self.node): raise ValueError("LB node index out of bounds") + + property index: + def __get__(self): + return (self.node[0], self.node[1], self.node[2]) property velocity: def __get__(self): diff --git a/testsuite/python/lb.py b/testsuite/python/lb.py index 6d047145483..90c3dbdada1 100644 --- a/testsuite/python/lb.py +++ b/testsuite/python/lb.py @@ -108,14 +108,6 @@ def test_mass_momentum_thermostat(self): all_temp_particle = [] all_temp_fluid = [] - # Cache the lb nodes - lb_nodes = [] - n_nodes = int(self.params['box_l'] / self.params['agrid']) - for i in range(n_nodes): - for j in range(n_nodes): - for k in range(n_nodes): - lb_nodes.append(self.lbf[i, j, k]) - # Integration for i in range(self.params['int_times']): self.system.integrator.run(self.params['int_steps']) @@ -125,14 +117,15 @@ def test_mass_momentum_thermostat(self): fluid_temp = 0.0 # Go over lb lattice - for lb_node in lb_nodes: + for lb_node in self.lbf.nodes(): dens = lb_node.density fluid_mass += dens fluid_temp += np.sum(lb_node.velocity**2) * dens # Normalize - fluid_mass /= len(lb_nodes) - fluid_temp *= self.system.volume() / (3. * len(lb_nodes)**2) + fluid_mass /= np.product(self.lbf.shape) + fluid_temp *= self.system.volume() / ( + 3. * np.product(self.lbf.shape)**2) # check mass conversation self.assertAlmostEqual(fluid_mass, self.params["dens"], @@ -190,10 +183,8 @@ def test_stress_tensor(self): system.integrator.run(10) stress = np.zeros((3, 3)) agrid = self.params["agrid"] - for i in range(int(system.box_l[0] / agrid)): - for j in range(int(system.box_l[1] / agrid)): - for k in range(int(system.box_l[2] / agrid)): - stress += self.lbf[i, j, k].stress + for n in self.lbf.nodes(): + stress += n.stress stress /= system.volume() / agrid**3 @@ -202,7 +193,6 @@ def test_stress_tensor(self): obs_stress = np.array([[obs_stress[0], obs_stress[1], obs_stress[3]], [obs_stress[1], obs_stress[2], obs_stress[4]], [obs_stress[3], obs_stress[4], obs_stress[5]]]) - print(stress / obs_stress) np.testing.assert_allclose(stress, obs_stress, atol=1E-10) def test_lb_node_set_get(self): @@ -215,6 +205,13 @@ def test_lb_node_set_get(self): tau=self.system.time_step, ext_force_density=[0, 0, 0]) self.system.actors.add(self.lbf) + + self.assertEqual(self.lbf.shape, + ( + int(self.system.box_l[0] / self.params["agrid"]), + int(self.system.box_l[1] / self.params["agrid"]), + int(self.system.box_l[2] / self.params["agrid"]))) + v_fluid = np.array([1.2, 4.3, 0.2]) self.lbf[0, 0, 0].velocity = v_fluid np.testing.assert_allclose( @@ -223,6 +220,8 @@ def test_lb_node_set_get(self): self.lbf[0, 0, 0].density = density self.assertAlmostEqual(self.lbf[0, 0, 0].density, density, delta=1e-4) + self.assertEqual(self.lbf[3, 2, 1].index, (3, 2, 1)) + def test_parameter_change_without_seed(self): self.system.actors.clear() self.lbf = self.lb_class( @@ -324,9 +323,9 @@ def test_a_ext_force_density(self): # (force is applied only to the second half of the first integration step) fluid_velocity = np.array(ext_force_density) * self.system.time_step * ( n_time_steps - 0.5) / self.params['dens'] - for n in list(itertools.combinations(range(int(self.system.box_l[0] / self.params['agrid'])), 3)): + for n in self.lbf.nodes(): np.testing.assert_allclose( - np.copy(self.lbf[n].velocity), fluid_velocity, atol=1E-6) + np.copy(n.velocity), fluid_velocity, atol=1E-6) class TestLBCPU(TestLB, ut.TestCase):