diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 1c02b277167..55f402d2194 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -192,6 +192,8 @@ python_test(FILE lb_boundary.py MAX_NUM_PROC 2 LABELS gpu) python_test(FILE lb_streaming.py MAX_NUM_PROC 4 LABELS gpu) python_test(FILE lb_shear.py MAX_NUM_PROC 2 LABELS gpu) python_test(FILE lb_thermostat.py MAX_NUM_PROC 2 LABELS gpu) +python_test(FILE lb_buoyancy_force.py MAX_NUM_PROC 4 LABELS gpu) +python_test(FILE lb_momentum_conservation.py MAX_NUM_PROC 4 LABELS gpu) python_test(FILE p3m_electrostatic_pressure.py MAX_NUM_PROC 2) python_test(FILE sigint.py DEPENDENCIES sigint_child.py MAX_NUM_PROC 1) python_test(FILE lb_density.py MAX_NUM_PROC 1) diff --git a/testsuite/python/lb_boundary_volume_force.py b/testsuite/python/lb_boundary_volume_force.py index d23ca756d79..1d154172dea 100644 --- a/testsuite/python/lb_boundary_volume_force.py +++ b/testsuite/python/lb_boundary_volume_force.py @@ -22,6 +22,8 @@ import espressomd.lbboundaries import espressomd.shapes +from tests_common import count_fluid_nodes + """ Checks force on lb boundaries for a fluid with a uniform volume force @@ -48,14 +50,6 @@ class LBBoundaryForceCommon: system.time_step = TIME_STEP system.cell_system.skin = 0.4 * AGRID - def count_fluid_nodes(self): - # Count non-boundary nodes: - fluid_nodes = 0 - for n in self.lbf.nodes(): - if not n.boundary: - fluid_nodes += 1 - return fluid_nodes - def test(self): """ Integrate the LB fluid until steady state is reached within a certain @@ -74,7 +68,7 @@ def test(self): self.system.lbboundaries.add(wall1) self.system.lbboundaries.add(wall2) - fluid_nodes = self.count_fluid_nodes() + fluid_nodes = count_fluid_nodes(self.lbf) self.system.integrator.run(20) diff = float("inf") diff --git a/testsuite/python/lb_buoyancy_force.py b/testsuite/python/lb_buoyancy_force.py new file mode 100644 index 00000000000..4c4945223ef --- /dev/null +++ b/testsuite/python/lb_buoyancy_force.py @@ -0,0 +1,128 @@ +# Copyright (C) 2010-2018 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import espressomd +from espressomd import lb, lbboundaries, shapes, has_features +import unittest as ut +import numpy as np +import unittest_decorators as utx + +from tests_common import count_fluid_nodes + +# Define the LB Parameters +TIME_STEP = 0.01 +AGRID = 0.5 +KVISC = 6 +DENS = 2 +G = 0.08 +BOX_SIZE = 24 * AGRID + +LB_PARAMS = {'agrid': AGRID, + 'dens': DENS, + 'visc': KVISC, + 'tau': TIME_STEP, + 'ext_force_density': [0, DENS * G, 0]} +# System setup +RADIUS = 8 * AGRID + + +class Buoyancy(object): + """ + Tests buoyancy force on a sphere in a closed box of lb fluid and + the overall force balance + + """ + lbf = None + system = espressomd.System(box_l=[BOX_SIZE] * 3) + system.time_step = TIME_STEP + system.cell_system.skin = 0.01 + + def test(self): + self.system.actors.clear() + self.system.lbboundaries.clear() + self.system.actors.add(self.lbf) + + # Setup walls + for i in range(3): + n = np.zeros(3) + n[i] = 1 + self.system.lbboundaries.add( + lbboundaries.LBBoundary(shape=shapes.Wall( + normal=-n, dist=-(self.system.box_l[i] - AGRID)))) + + self.system.lbboundaries.add(lbboundaries.LBBoundary( + shape=shapes.Wall( + normal=n, dist=AGRID))) + + # setup sphere without slip in the middle + sphere = lbboundaries.LBBoundary(shape=shapes.Sphere( + radius=RADIUS, center=self.system.box_l / 2, direction=1)) + + self.system.lbboundaries.add(sphere) + + sphere_volume = 4. / 3. * np.pi * RADIUS**3 + + # Equilibration + last_force = -999999 + self.system.integrator.run(10) + while True: + self.system.integrator.run(10) + force = np.linalg.norm(sphere.get_force()) + + if np.linalg.norm(force - last_force) < 0.01: + break + last_force = force + + # Check force balance + boundary_force = np.zeros(3) + for b in self.system.lbboundaries: + boundary_force += b.get_force() + + fluid_nodes = count_fluid_nodes(self.lbf) + fluid_volume = fluid_nodes * AGRID**3 + applied_force = fluid_volume * np.array(LB_PARAMS['ext_force_density']) + + np.testing.assert_allclose( + boundary_force, + applied_force, + atol=0.08 * np.linalg.norm(applied_force)) + + # Check buoyancy force on the sphere + expected_force = np.array( + [0, -sphere_volume * DENS * G, 0]) + np.testing.assert_allclose( + np.copy(sphere.get_force()), expected_force, + atol=np.linalg.norm(expected_force) * 0.02) + + +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures(["LB_BOUNDARIES_GPU", "EXTERNAL_FORCES"]) +class LBGPUBuoyancy(ut.TestCase, Buoyancy): + + def setUp(self): + self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS) + + +@utx.skipIfMissingFeatures(["LB_BOUNDARIES", "EXTERNAL_FORCES"]) +class LBCPUBuoyancy(ut.TestCase, Buoyancy): + + def setUp(self): + self.lbf = espressomd.lb.LBFluid(**LB_PARAMS) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/lb_momentum_conservation.py b/testsuite/python/lb_momentum_conservation.py new file mode 100644 index 00000000000..1dd7c85bf6d --- /dev/null +++ b/testsuite/python/lb_momentum_conservation.py @@ -0,0 +1,103 @@ +# Copyright (C) 2010-2018 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import espressomd +from espressomd import lb, lbboundaries, shapes, has_features +import unittest as ut +import numpy as np +import sys + +# Define the LB Parameters +TIME_STEP = 0.1 +AGRID = 1.0 +KVISC = 5 +DENS = 1 +BOX_SIZE = 6 * AGRID +F = 1. / BOX_SIZE**3 +GAMMA = 15 + +LB_PARAMS = {'agrid': AGRID, + 'dens': DENS, + 'visc': KVISC, + 'tau': TIME_STEP, + 'ext_force_density': [0, F, 0]} + + +class Momentum(object): + """ + Tests momentum conservatoin for an LB coupled to a particle, where opposing + forces are applied to LB and particle. The test should uncover issues + with boundary and ghost layer handling. + + """ + lbf = None + system = espressomd.System(box_l=[BOX_SIZE] * 3) + system.time_step = TIME_STEP + system.cell_system.skin = 0.01 + + def test(self): + self.system.actors.clear() + self.system.part.clear() + self.system.actors.add(self.lbf) + self.system.thermostat.set_lb(LB_fluid=self.lbf, gamma=GAMMA, seed=1) + + applied_force = self.system.volume() * np.array( + LB_PARAMS['ext_force_density']) + p = self.system.part.add( + pos=(0, 0, 0), ext_force=-applied_force, v=[.1, .2, .3]) + + # Reach steady state + self.system.integrator.run(500) + v_final = np.copy(p.v) + momentum = self.system.analysis.linear_momentum() + + for i in range(10): + self.system.integrator.run(50) + # check that momentum stays constant + np.testing.assert_allclose( + self.system.analysis.linear_momentum(), momentum, atol=2E-4) + + # Check that particle velocity is stationary + # up to the acceleration of 1/2 time step + np.testing.assert_allclose(np.copy(p.v), v_final, atol=2.2E-3) + + # Make sure, the particle has crossed the periodic boundaries + self.assertGreater( + np.amax( + np.abs(v_final) * + self.system.time), + BOX_SIZE) + + +@ut.skipIf(not espressomd.gpu_available() or not espressomd.has_features( + ['EXTERNAL_FORCES']), "Skipping test due to missing features.") +class LBGPUMomentum(ut.TestCase, Momentum): + + def setUp(self): + self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS) + + +@ut.skipIf(not espressomd.has_features( + ['EXTERNAL_FORCES']), "Skipping test due to missing features.") +class LBCPUMomentum(ut.TestCase, Momentum): + + def setUp(self): + self.lbf = espressomd.lb.LBFluid(**LB_PARAMS) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py index 8955c76f732..f04dca7c4ee 100644 --- a/testsuite/python/tests_common.py +++ b/testsuite/python/tests_common.py @@ -614,3 +614,14 @@ def single_component_maxwell(x1, x2, kT): def lists_contain_same_elements(list1, list2): return len(list1) == len(list2) and sorted(list1) == sorted(list2) + + +def count_fluid_nodes(lbf): + """Counts the non-boundary nodes in the passed lb fluid instance.""" + + fluid_nodes = 0 + for n in lbf.nodes(): + if not n.boundary: + fluid_nodes += 1 + + return fluid_nodes