From cbe47b992fbda65d3a7e794704750148486d9195 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 2 Nov 2018 18:54:00 +0100 Subject: [PATCH] Merge pull request #2362 from KaiSzuttor/lb_thermostat LB thermostat test --- testsuite/python/CMakeLists.txt | 1 + testsuite/python/dpd.py | 16 ++-- testsuite/python/langevin_thermostat.py | 11 +-- testsuite/python/lb_thermostat.py | 101 ++++++++++++++++++++++++ testsuite/python/tests_common.py | 13 +++ testsuite/python/thermalized_bond.py | 10 +-- 6 files changed, 126 insertions(+), 26 deletions(-) create mode 100644 testsuite/python/lb_thermostat.py diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 968aeb23988..ced330772ce 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -166,6 +166,7 @@ python_test(FILE field_test.py MAX_NUM_PROC 1) python_test(FILE lb_boundary.py MAX_NUM_PROC 2) python_test(FILE lb_streaming.py MAX_NUM_PROC 4) python_test(FILE lb_shear.py MAX_NUM_PROC 2) +python_test(FILE lb_thermostat.py MAX_NUM_PROC 2) if(PY_H5PY) python_test(FILE h5md.py MAX_NUM_PROC 2) diff --git a/testsuite/python/dpd.py b/testsuite/python/dpd.py index 558698b744c..040afa3632d 100644 --- a/testsuite/python/dpd.py +++ b/testsuite/python/dpd.py @@ -18,12 +18,13 @@ # # Tests particle property setters/getters from __future__ import print_function -import unittest as ut -import espressomd import numpy as np -from time import time +import unittest as ut from itertools import product +from time import time +import espressomd +from tests_common import single_component_maxwell @ut.skipIf(not espressomd.has_features("DPD"), "Skipped because feature is disabled") class DPDThermostat(ut.TestCase): @@ -43,11 +44,6 @@ def tearDown(self): s = self.s s.part.clear() - def single_component_maxwell(self, x1, x2, kT): - """Integrate the probability density from x1 to x2 using the trapez rule""" - x = np.linspace(x1, x2, 1000) - return np.trapz(np.exp(-x**2 / (2.*kT)), x)/np.sqrt(2.*np.pi*kT) - def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT): """check the recorded particle distributions in vel againsta histogram with n_bins bins. Drop velocities outside minmax. Check individual histogram bins up to an accuracy of error_tol agaisnt the analytical result for kT.""" for i in range(3): @@ -56,13 +52,13 @@ def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT): bins = hist[1] for j in range(n_bins): found = data[j] - expected = self.single_component_maxwell(bins[j], bins[j+1], kT) + expected = single_component_maxwell(bins[j], bins[j+1], kT) self.assertLessEqual(abs(found - expected), error_tol) def test_aa_verify_single_component_maxwell(self): """Verifies the normalization of the analytical expression.""" self.assertLessEqual( - abs(self.single_component_maxwell(-10, 10, 4.)-1.), 1E-4) + abs(single_component_maxwell(-10, 10, 4.)-1.), 1E-4) def check_total_zero(self): v_total = np.sum(self.s.part[:].v, axis=0) diff --git a/testsuite/python/langevin_thermostat.py b/testsuite/python/langevin_thermostat.py index 0a40046ee0c..933e9f174c9 100644 --- a/testsuite/python/langevin_thermostat.py +++ b/testsuite/python/langevin_thermostat.py @@ -25,6 +25,7 @@ from time import time from espressomd.accumulators import Correlator from espressomd.observables import ParticleVelocities, ParticleBodyAngularVelocities +from tests_common import single_component_maxwell @ut.skipIf(espressomd.has_features("THERMOSTAT_IGNORE_NON_VIRTUAL"), @@ -45,12 +46,6 @@ class LangevinThermostat(ut.TestCase): def setUpClass(cls): np.random.seed(42) - def single_component_maxwell(self, x1, x2, kT): - """Integrate the probability density from x1 to x2 using the trapez rule""" - x = np.linspace(x1, x2, 1000) - return np.trapz(np.exp(-x**2 / (2. * kT)), x) / \ - np.sqrt(2. * np.pi * kT) - def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT): """check the recorded particle distributions in vel againsta histogram with n_bins bins. Drop velocities outside minmax. Check individual histogram bins up to an accuracy of error_tol agaisnt the analytical result for kT.""" for i in range(3): @@ -60,14 +55,14 @@ def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT): bins = hist[1] for j in range(n_bins): found = data[j] - expected = self.single_component_maxwell( + expected = single_component_maxwell( bins[j], bins[j + 1], kT) self.assertLessEqual(abs(found - expected), error_tol) def test_aa_verify_single_component_maxwell(self): """Verifies the normalization of the analytical expression.""" self.assertLessEqual( - abs(self.single_component_maxwell(-10, 10, 4.) - 1.), 1E-4) + abs(single_component_maxwell(-10, 10, 4.) - 1.), 1E-4) def test_global_langevin(self): """Test for global Langevin parameters.""" diff --git a/testsuite/python/lb_thermostat.py b/testsuite/python/lb_thermostat.py new file mode 100644 index 00000000000..ec720c3b5b3 --- /dev/null +++ b/testsuite/python/lb_thermostat.py @@ -0,0 +1,101 @@ +# 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 unittest as ut +import numpy as np + +import espressomd.lb +from tests_common import single_component_maxwell + +""" +Check the Lattice Boltzmann thermostat with respect to the particle velocity distribution. + + +""" + +KT = 2.25 +AGRID = 2.5 +VISC = .7 +DENS = 1.7 +TIME_STEP = 0.01 +LB_PARAMS = {'agrid': AGRID, + 'dens': DENS, + 'visc': VISC, + 'fric': 2.0, + 'tau': TIME_STEP} + + +class LBThermostatCommon(object): + + """Base class of the test that holds the test logic.""" + lbf = None + system = espressomd.System(box_l=[10.0, 10.0, 10.0]) + system.time_step = TIME_STEP + system.cell_system.skin = 0.4 * AGRID + + def prepare(self): + self.system.set_random_state_PRNG() + self.system.actors.clear() + self.system.actors.add(self.lbf) + self.system.part.add( + pos=np.random.random((250, 3)) * self.system.box_l) + self.system.thermostat.set_lb(kT=KT) + + def test_velocity_distribution(self): + self.prepare() + self.system.integrator.run(100) + N = len(self.system.part) + loops = 250 + v_stored = np.zeros((N * loops, 3)) + for i in range(loops): + self.system.integrator.run(10) + v_stored[i * N:(i + 1) * N, :] = self.system.part[:].v + minmax = 5 + n_bins = 5 + error_tol = 0.01 + for i in range(3): + hist = np.histogram(v_stored[:, i], range=( + -minmax, minmax), bins=n_bins, normed=False) + data = hist[0] / float(v_stored.shape[0]) + bins = hist[1] + for j in range(n_bins): + found = data[j] + expected = single_component_maxwell(bins[j], bins[j + 1], KT) + self.assertLessEqual(abs(found - expected), error_tol) + + +@ut.skipIf(not espressomd.has_features( + ['LB']), "Skipping test due to missing features.") +class LBCPUThermostat(ut.TestCase, LBThermostatCommon): + + """Test for the CPU implementation of the LB.""" + + def setUp(self): + self.lbf = espressomd.lb.LBFluid(**LB_PARAMS) + + +@ut.skipIf(not espressomd.has_features( + ['LB_GPU']), "Skipping test due to missing features.") +class LBGPUThermostat(ut.TestCase, LBThermostatCommon): + + """Test for the GPU implementation of the LB.""" + + def setUp(self): + self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS) + + +if __name__ == '__main__': + ut.main() diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py index af3e92d0c97..c87132c3afe 100644 --- a/testsuite/python/tests_common.py +++ b/testsuite/python/tests_common.py @@ -581,3 +581,16 @@ def gaussian_force(r, eps, sig, cutoff): if (r < cutoff): f = eps * r / sig**2 * np.exp(-np.power(r / sig, 2) / 2) return f + + +class DynamicDict(dict): + + def __getitem__(self, key): + value = super(DynamicDict, self).__getitem__(key) + return eval(value, self) if isinstance(value, str) else value + + +def single_component_maxwell(x1, x2, kT): + """Integrate the probability density from x1 to x2 using the trapezoidal rule""" + x = np.linspace(x1, x2, 1000) + return np.trapz(np.exp(-x**2 / (2. * kT)), x) / np.sqrt(2. * np.pi * kT) diff --git a/testsuite/python/thermalized_bond.py b/testsuite/python/thermalized_bond.py index 52daa2c219a..f68ed1d9fc3 100644 --- a/testsuite/python/thermalized_bond.py +++ b/testsuite/python/thermalized_bond.py @@ -23,6 +23,7 @@ import unittest as ut import espressomd +from tests_common import single_component_maxwell @ut.skipIf(not espressomd.has_features(["MASS"]), "Features not available, skipping test!") class ThermalizedBond(ut.TestCase): @@ -39,13 +40,6 @@ class ThermalizedBond(ut.TestCase): def setUpClass(cls): np.random.seed(42) - def single_component_maxwell(self, x1, x2, kT): - """Integrate the probability density from x1 to x2 using the trapez rule""" - - x = np.linspace(x1, x2, 1000) - return np.trapz(np.exp(-x**2 / (2. * kT)), x) / \ - np.sqrt(2. * np.pi * kT) - def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT): """check the recorded particle distributions in vel againsta histogram with n_bins bins. Drop velocities outside minmax. Check individual histogram bins up to an accuracy of error_tol agaisnt the analytical result for kT.""" @@ -57,7 +51,7 @@ def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT): for j in range(n_bins): found = data[j] - expected = self.single_component_maxwell( + expected = single_component_maxwell( bins[j], bins[j + 1], kT) self.assertLessEqual(abs(found - expected), error_tol)