From 4fc322c8af0d0250457647ed588e1cf29cf52a5f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 2 Jul 2020 19:08:16 +0200 Subject: [PATCH] SD: Test CPU/GPU, FT/FTS, checkpointing, sample --- samples/dancing.py | 80 +++++++++++--- testsuite/python/CMakeLists.txt | 10 +- testsuite/python/save_checkpoint.py | 13 +++ testsuite/python/stokesian_dynamics.py | 117 ++++++++++----------- testsuite/python/stokesian_dynamics_cpu.py | 58 ++++++++++ testsuite/python/stokesian_dynamics_gpu.py | 55 ++++++++++ testsuite/python/test_checkpoint.py | 41 ++++++++ testsuite/scripts/samples/CMakeLists.txt | 4 + testsuite/scripts/samples/test_dancing.py | 63 +++++++++++ 9 files changed, 362 insertions(+), 79 deletions(-) create mode 100644 testsuite/python/stokesian_dynamics_cpu.py create mode 100644 testsuite/python/stokesian_dynamics_gpu.py create mode 100644 testsuite/scripts/samples/test_dancing.py diff --git a/samples/dancing.py b/samples/dancing.py index cb8d41afafb..dcbcd780275 100644 --- a/samples/dancing.py +++ b/samples/dancing.py @@ -16,34 +16,86 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . # - +""" +Stokesian Dynamics simulation of particle sedimentation. +Reproduce the trajectory in Figure 5b from :cite:`durlofsky87a`. +""" import espressomd -from espressomd import constraints +import espressomd.constraints +import espressomd.observables +import espressomd.accumulators import numpy as np import matplotlib.pyplot as plt -espressomd.assert_features("STOKESIAN_DYNAMICS") +import argparse + +parser = argparse.ArgumentParser(epilog=__doc__) +group = parser.add_mutually_exclusive_group() +group.add_argument('--cpu', action='store_true', help='Use CPU implementation') +group.add_argument('--gpu', action='store_true', help='Use GPU implementation') +group = parser.add_mutually_exclusive_group() +group.add_argument('--ft', action='store_true', help='Use FT approximation') +group.add_argument('--fts', action='store_true', help='Use FTS approximation') +args = parser.parse_args() + +required_features = ["STOKESIAN_DYNAMICS"] + +if args.gpu: + print("Using GPU implementation") + required_features.append("CUDA") + required_features.append("STOKESIAN_DYNAMICS_GPU") + sd_device = "gpu" +else: + print("Using CPU implementation") + sd_device = "cpu" + +if args.ft: + print("Using FT approximation method") + sd_method = "ft" +else: + print("Using FTS approximation method") + sd_method = "fts" + +espressomd.assert_features(required_features) system = espressomd.System(box_l=[10, 10, 10]) -system.time_step = 1.0 +system.time_step = 1.5 system.cell_system.skin = 0.4 +system.periodicity = [False, False, False] -system.integrator.set_sd(viscosity=1.0, radii={0: 1.0}) +system.integrator.set_sd( + viscosity=1.0, radii={0: 1.0}, device=sd_device, + approximation_method=sd_method) system.part.add(pos=[-5, 0, 0], rotation=[1, 1, 1]) system.part.add(pos=[0, 0, 0], rotation=[1, 1, 1]) system.part.add(pos=[7, 0, 0], rotation=[1, 1, 1]) -gravity = constraints.Gravity(g=[0, -1, 0]) +gravity = espressomd.constraints.Gravity(g=[0, -1, 0]) system.constraints.add(gravity) -intsteps = int(13000 / system.time_step) -pos = np.empty([intsteps, 3 * len(system.part)]) -for i in range(intsteps): - system.integrator.run(1) - for n, p in enumerate(system.part): - pos[i, 3 * n:3 * n + 3] = p.pos +obs = espressomd.observables.ParticlePositions(ids=system.part[:].id) +acc = espressomd.accumulators.TimeSeries(obs=obs, delta_N=1) +system.auto_update_accumulators.add(acc) +acc.update() +intsteps = int(10500 / system.time_step) +system.integrator.run(intsteps) + +positions = acc.time_series() +ref_data = "../testsuite/python/data/dancing.txt" +data = np.loadtxt(ref_data) + +for i in range(3): + plt.plot(positions[:, i, 0], positions[:, i, 1], linestyle='solid') + +plt.gca().set_prop_cycle(None) + +for i in range(0, 6, 2): + plt.plot(data[:, i], data[:, i + 1], linestyle='dashed') -for n, p in enumerate(system.part): - plt.plot(pos[:, 3 * n], pos[:, 3 * n + 1]) +plt.title("Trajectory of sedimenting spheres\nsolid line: simulation " + "({} on {}), dashed line: paper (FTS)" + .format(sd_method.upper(), sd_device.upper())) +plt.xlabel("x") +plt.ylabel("y") plt.show() diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 8bb88dfd4fd..37ff49f305e 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -56,7 +56,7 @@ endfunction(PYTHON_TEST) # Separate features with hyphens, use a period to add an optional flag. foreach( TEST_COMBINATION - lb.cpu-p3m.cpu-lj-therm.lb;lb.gpu-p3m.cpu-lj-therm.lb;ek.gpu;lb.off-therm.npt-int.npt-minimization;lb.off-therm.langevin-int.nvt;lb.off-therm.dpd-int.sd;lb.off-therm.bd-int.bd + lb.cpu-p3m.cpu-lj-therm.lb;lb.gpu-p3m.cpu-lj-therm.lb;ek.gpu;lb.off-therm.npt-int.npt-minimization;lb.off-therm.langevin-int.nvt;lb.off-therm.dpd-int.sd;lb.off-therm.bd-int.bd;lb.off-therm.sdm-int.sdm.cpu;lb.off-therm.sdm-int.sdm.gpu ) if(${TEST_COMBINATION} MATCHES "\\.gpu") set(TEST_LABELS "gpu") @@ -200,10 +200,11 @@ python_test(FILE gpu_availability.py MAX_NUM_PROC 1 LABELS gpu) python_test(FILE features.py MAX_NUM_PROC 1) python_test(FILE galilei.py MAX_NUM_PROC 32) python_test(FILE time_series.py MAX_NUM_PROC 1) -python_test(FILE linear_momentum.py) +python_test(FILE linear_momentum.py MAX_NUM_PROC 4) python_test(FILE linear_momentum_lb.py MAX_NUM_PROC 2 LABELS gpu) python_test(FILE mmm1d.py MAX_NUM_PROC 2 LABELS gpu) -python_test(FILE stokesian_dynamics.py) +python_test(FILE stokesian_dynamics_cpu.py MAX_NUM_PROC 2) +python_test(FILE stokesian_dynamics_gpu.py MAX_NUM_PROC 4 LABELS gpu) python_test(FILE elc.py MAX_NUM_PROC 2) python_test(FILE elc_vs_analytic.py MAX_NUM_PROC 2) python_test(FILE rotation.py MAX_NUM_PROC 1) @@ -225,6 +226,9 @@ add_custom_target( ${CMAKE_CURRENT_BINARY_DIR} COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/ek_common.py ${CMAKE_CURRENT_BINARY_DIR} + COMMAND + ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/stokesian_dynamics.py + ${CMAKE_CURRENT_BINARY_DIR} COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/ek_eof_one_species_base.py diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index 78ffe793041..62576b536e2 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -153,6 +153,9 @@ system.thermostat.set_npt(kT=1.0, gamma0=2.0, gammav=0.1, seed=42) elif 'THERM.DPD' in modes and has_features('DPD'): system.thermostat.set_dpd(kT=1.0, seed=42) + elif 'THERM.SDM' in modes and (has_features('STOKESIAN_DYNAMICS') or has_features('STOKESIAN_DYNAMICS_GPU')): + system.periodicity = [0, 0, 0] + system.thermostat.set_sd(kT=1.0, seed=42) # set integrator if 'INT.NPT' in modes and has_features('NPT'): system.integrator.set_isotropic_npt(ext_pressure=2.0, piston=0.01, @@ -164,6 +167,16 @@ system.integrator.set_nvt() elif 'INT.BD' in modes: system.integrator.set_brownian_dynamics() + elif 'INT.SDM.CPU' in modes and has_features('STOKESIAN_DYNAMICS'): + system.periodicity = [0, 0, 0] + system.integrator.set_sd(approximation_method='ft', viscosity=0.5, + device='cpu', radii={0: 1.5}, + pair_mobility=False, self_mobility=True) + elif 'INT.SDM.GPU' in modes and has_features('STOKESIAN_DYNAMICS_GPU') and espressomd.gpu_available(): + system.periodicity = [0, 0, 0] + system.integrator.set_sd(approximation_method='fts', viscosity=2.0, + device='gpu', radii={0: 1.0}, + pair_mobility=True, self_mobility=False) # set minimization if 'MINIMIZATION' in modes: steepest_descent(system, f_max=1, gamma=10, max_steps=0, diff --git a/testsuite/python/stokesian_dynamics.py b/testsuite/python/stokesian_dynamics.py index 9720eda8989..22f2d60fa9c 100644 --- a/testsuite/python/stokesian_dynamics.py +++ b/testsuite/python/stokesian_dynamics.py @@ -17,26 +17,19 @@ # along with this program. If not, see . # import espressomd -from espressomd import constraints +import espressomd.observables +import espressomd.accumulators +import espressomd.constraints import numpy as np import unittest as ut -import unittest_decorators as utx from tests_common import abspath s = espressomd.System(box_l=[1.0, 1.0, 1.0]) -def skipIfMissingFeatureStokesianDynamics(): - """Specialized unittest skipIf decorator for missing Stokesian Dynamics.""" - if not espressomd.has_features(["STOKESIAN_DYNAMICS"]) and (not espressomd.has_features( - ["STOKESIAN_DYNAMICS_GPU"]) or not espressomd.gpu_available()): - return ut.skip("Skipping test: feature STOKESIAN_DYNAMICS unavailable") - return utx._id - - -@skipIfMissingFeatureStokesianDynamics() class StokesianDynamicsSetupTest(ut.TestCase): system = s + device = 'none' def setUp(self): self.system.box_l = [10] * 3 @@ -44,30 +37,29 @@ def setUp(self): self.system.time_step = 1.0 self.system.cell_system.skin = 0.4 - # unset SD integrator so we can test whether set_sd fails + # unset SD integrator so we can test whether set_sd() fails. # set_nvt() is the only way to ensure that integ_switch is # set to a different value than INTEG_METHOD_SD self.system.integrator.set_nvt() - def test_pbc_checks(self): - + def pbc_checks(self): self.system.periodicity = [0, 0, 1] - with (self.assertRaises(Exception)): + with self.assertRaises(Exception): self.system.integrator.set_sd(viscosity=1.0, - device="cpu", + device=self.device, radii={0: 1.0}) self.system.periodicity = [0, 0, 0] self.system.integrator.set_sd(viscosity=1.0, - device="cpu", + device=self.device, radii={0: 1.0}) - with (self.assertRaises(Exception)): + with self.assertRaises(Exception): self.system.periodicity = [0, 1, 0] -@skipIfMissingFeatureStokesianDynamics() class StokesianDynamicsTest(ut.TestCase): system = s + device = 'none' # Digitized reference data of Figure 5b from # Durlofsky et al., J. Fluid Mech. 180, 21 (1987) @@ -79,61 +71,66 @@ def setUp(self): self.system.periodicity = [0, 0, 0] self.system.cell_system.skin = 0.4 - def falling_spheres(self, time_step, l_factor, t_factor): + def falling_spheres(self, time_step, l_factor, t_factor, + sd_method='fts', sd_short=False): self.system.time_step = time_step self.system.part.add(pos=[-5 * l_factor, 0, 0], rotation=[1, 1, 1]) self.system.part.add(pos=[0 * l_factor, 0, 0], rotation=[1, 1, 1]) self.system.part.add(pos=[7 * l_factor, 0, 0], rotation=[1, 1, 1]) - self.system.integrator.set_sd(viscosity=1.0 / (t_factor * l_factor), - device="cpu", - radii={0: 1.0 * l_factor}) + self.system.integrator.set_sd( + viscosity=1.0 / (t_factor * l_factor), + device=self.device, radii={0: 1.0 * l_factor}, + approximation_method=sd_method) - gravity = constraints.Gravity( + gravity = espressomd.constraints.Gravity( g=[0, -1.0 * l_factor / (t_factor**2), 0]) self.system.constraints.add(gravity) self.system.time_step = 1.0 * t_factor - intsteps = int(8000 * t_factor / self.system.time_step) - pos = np.empty([intsteps + 1, 3 * len(self.system.part)]) - - pos[0, :] = self.system.part[:].pos.flatten() - for i in range(intsteps): - self.system.integrator.run(1) - for n, p in enumerate(self.system.part): - pos[i + 1, 3 * n:3 * n + 3] = p.pos - - for i in range(self.data.shape[0]): - for n in range(3): - x_ref = self.data[i, 2 * n] * l_factor - y_ref = self.data[i, 2 * n + 1] * l_factor - x = pos[:, 3 * n] - y = pos[:, 3 * n + 1] - - if y_ref < -555 * l_factor: - continue - - idx = np.abs(y - y_ref).argmin() - dist = np.sqrt((x_ref - x[idx])**2 + (y_ref - y[idx])**2) - self.assertLess(dist, 0.5 * l_factor) - - def test_default(self): - self.falling_spheres(1.0, 1.0, 1.0) - - def test_rescaled(self): - self.falling_spheres(1.0, 4.5, 2.5) - - def test_different_time_step(self): - self.falling_spheres(0.7, 1.0, 1.0) + obs = espressomd.observables.ParticlePositions(ids=(0, 1, 2)) + acc = espressomd.accumulators.TimeSeries(obs=obs, delta_N=1) + self.system.auto_update_accumulators.add(acc) + acc.update() + + if sd_short: + y_min = -80 + intsteps = 1300 + elif sd_method == 'fts': + y_min = -555 + intsteps = 8000 + else: + y_min = -200 + intsteps = 3000 + intsteps = int(intsteps * t_factor / self.system.time_step) + + self.system.integrator.run(intsteps) + + simul = acc.time_series()[:, :, 0:2] + paper = self.data.reshape([-1, 3, 2]) + + for pid in range(3): + dist = [] + # the simulated trajectory is oversampled by a ratio of + # (90/t_factor):1 compared to the published trajectory + for desired in paper[:, pid] * l_factor: + if desired[1] < y_min * l_factor: + break + # find the closest point in the simulated trajectory + idx = np.abs(simul[:, pid, 1] - desired[1]).argmin() + actual = simul[idx, pid] + dist.append(np.linalg.norm(actual - desired)) + self.assertLess(idx, intsteps, msg='Insufficient sampling') + np.testing.assert_allclose(dist, 0, rtol=0, atol=0.5 * l_factor) def tearDown(self): self.system.constraints.clear() self.system.part.clear() -@skipIfMissingFeatureStokesianDynamics() class StokesianDiffusionTest(ut.TestCase): system = s + device = 'none' kT = 1e-4 R = 1.5 @@ -148,12 +145,12 @@ def setUp(self): self.system.part.add(pos=[0, 0, 0], rotation=[1, 1, 1]) + def check(self): self.system.integrator.set_sd(viscosity=self.eta, - device="cpu", + device=self.device, radii={0: self.R}) self.system.thermostat.set_sd(kT=self.kT, seed=42) - def test(self): intsteps = int(100000 / self.system.time_step) pos = np.empty([intsteps + 1, 3]) orientation = np.empty((intsteps + 1, 3)) @@ -206,7 +203,3 @@ def tearDown(self): self.system.constraints.clear() self.system.part.clear() self.system.thermostat.set_sd(kT=0) - - -if __name__ == '__main__': - ut.main() diff --git a/testsuite/python/stokesian_dynamics_cpu.py b/testsuite/python/stokesian_dynamics_cpu.py new file mode 100644 index 00000000000..4959da98631 --- /dev/null +++ b/testsuite/python/stokesian_dynamics_cpu.py @@ -0,0 +1,58 @@ +# +# Copyright (C) 2019-2020 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 unittest_decorators as utx +import stokesian_dynamics as sd + + +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS"]) +class StokesianDynamicsSetupTest(sd.StokesianDynamicsSetupTest): + device = 'cpu' + + def test_pbc_checks(self): + self.pbc_checks() + + +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS"]) +class StokesianDynamicsTest(sd.StokesianDynamicsTest): + device = 'cpu' + + def test_default(self): + self.falling_spheres(1.0, 1.0, 1.0) + + def test_rescaled(self): + self.falling_spheres(1.0, 4.5, 2.5) + + def test_different_time_step(self): + self.falling_spheres(0.7, 1.0, 1.0) + + def test_default_ft(self): + self.falling_spheres(1.0, 1.0, 1.0, 'ft') + + +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS"]) +class StokesianDiffusionTest(sd.StokesianDiffusionTest): + device = 'cpu' + + def test(self): + self.check() + + +if __name__ == '__main__': + ut.main() diff --git a/testsuite/python/stokesian_dynamics_gpu.py b/testsuite/python/stokesian_dynamics_gpu.py new file mode 100644 index 00000000000..b9ca6d0a73e --- /dev/null +++ b/testsuite/python/stokesian_dynamics_gpu.py @@ -0,0 +1,55 @@ +# +# Copyright (C) 2019-2020 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 unittest_decorators as utx +import stokesian_dynamics as sd + + +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS_GPU"]) +class StokesianDynamicsSetupTest(sd.StokesianDynamicsSetupTest): + device = 'gpu' + + def test_pbc_checks(self): + self.pbc_checks() + + +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS_GPU"]) +class StokesianDynamicsTest(sd.StokesianDynamicsTest): + device = 'gpu' + + def test_default_fts(self): + self.falling_spheres(1.0, 1.0, 1.0, 'fts', sd_short=True) + + def test_default_ft(self): + self.falling_spheres(1.0, 1.0, 1.0, 'ft', sd_short=True) + + +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS_GPU"]) +class StokesianDiffusionTest(sd.StokesianDiffusionTest): + device = 'gpu' + + def test(self): + pass + + +if __name__ == '__main__': + ut.main() diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 4faf7b9d1f9..1ebba9fae27 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -35,6 +35,14 @@ and espressomd.has_features('ELECTROKINETICS')) +def skipIfMissingFeatureStokesianDynamics(): + """Specialized unittest skipIf decorator for missing Stokesian Dynamics.""" + if not espressomd.has_features(["STOKESIAN_DYNAMICS"]) and (not espressomd.has_features( + ["STOKESIAN_DYNAMICS_GPU"]) or not espressomd.gpu_available()): + return ut.skip("Skipping test: feature STOKESIAN_DYNAMICS unavailable") + return utx._id + + class CheckpointTest(ut.TestCase): @classmethod @@ -191,6 +199,14 @@ def test_thermostat_NPT(self): self.assertEqual(thmst['gamma0'], 2.0) self.assertEqual(thmst['gammav'], 0.1) + @skipIfMissingFeatureStokesianDynamics() + @ut.skipIf('THERM.SDM' not in modes, 'SDM thermostat not in modes') + def test_thermostat_SDM(self): + thmst = system.thermostat.get_state()[0] + self.assertEqual(thmst['type'], 'SD') + self.assertEqual(thmst['kT'], 1.0) + self.assertEqual(thmst['seed'], 42) + @utx.skipIfMissingFeatures('NPT') @ut.skipIf('INT.NPT' not in modes, 'NPT integrator not in modes') def test_integrator_NPT(self): @@ -233,6 +249,31 @@ def test_integrator_BD(self): params = integ.get_params() self.assertEqual(params, {}) + @utx.skipIfMissingFeatures('STOKESIAN_DYNAMICS') + @ut.skipIf('INT.SDM.CPU' not in modes, 'SDM CPU integrator not in modes') + def test_integrator_SDM_cpu(self): + integ = system.integrator.get_state() + self.assertIsInstance(integ, espressomd.integrate.StokesianDynamics) + expected_params = { + 'approximation_method': 'ft', 'device': 'cpu', 'radii': {0: 1.5}, + 'viscosity': 0.5, 'lubrication': False, 'pair_mobility': False, + 'self_mobility': True} + params = integ.get_params() + self.assertEqual(params, expected_params) + + @utx.skipIfMissingGPU() + @utx.skipIfMissingFeatures('STOKESIAN_DYNAMICS_GPU') + @ut.skipIf('INT.SDM.GPU' not in modes, 'SDM GPU integrator not in modes') + def test_integrator_SDM_gpu(self): + integ = system.integrator.get_state() + self.assertIsInstance(integ, espressomd.integrate.StokesianDynamics) + expected_params = { + 'approximation_method': 'fts', 'device': 'gpu', 'radii': {0: 1.0}, + 'viscosity': 2.0, 'lubrication': False, 'pair_mobility': True, + 'self_mobility': False} + params = integ.get_params() + self.assertEqual(params, expected_params) + @utx.skipIfMissingFeatures('LENNARD_JONES') @ut.skipIf('LJ' not in modes, "Skipping test due to missing mode.") def test_non_bonded_inter(self): diff --git a/testsuite/scripts/samples/CMakeLists.txt b/testsuite/scripts/samples/CMakeLists.txt index 7fd4d80f588..803d6aa0e83 100644 --- a/testsuite/scripts/samples/CMakeLists.txt +++ b/testsuite/scripts/samples/CMakeLists.txt @@ -24,6 +24,10 @@ add_custom_target( sample_test(FILE test_billiard.py) sample_test(FILE test_chamber_game.py) sample_test(FILE test_constraints.py) +sample_test(FILE test_dancing.py SUFFIX cpu_ft) +sample_test(FILE test_dancing.py SUFFIX cpu_fts) +sample_test(FILE test_dancing.py SUFFIX gpu_ft LABELS "gpu") +sample_test(FILE test_dancing.py SUFFIX gpu_fts LABELS "gpu") sample_test(FILE test_diffusion_coefficient.py) sample_test(FILE test_dpd.py) sample_test(FILE test_drude_bmimpf6.py SUFFIX cpu) diff --git a/testsuite/scripts/samples/test_dancing.py b/testsuite/scripts/samples/test_dancing.py new file mode 100644 index 00000000000..1d322504438 --- /dev/null +++ b/testsuite/scripts/samples/test_dancing.py @@ -0,0 +1,63 @@ +# Copyright (C) 2020 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 importlib_wrapper +import numpy as np + +sd_device = "--gpu" if "gpu" in "@TEST_SUFFIX@" else "--cpu" +sd_method = "--fts" if "fts" in "@TEST_SUFFIX@" else "--ft" +if sd_method == "--fts" and sd_device == "--cpu": + y_min = -555 + intsteps = 5500 +else: + y_min = -200 + intsteps = 2100 + +sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import( + "@SAMPLES_DIR@/dancing.py", script_suffix="@TEST_SUFFIX@", + gpu=(sd_device == "--gpu"), cmd_arguments=[sd_device, sd_method], + intsteps=intsteps, + ref_data="@CMAKE_SOURCE_DIR@/testsuite/python/data/dancing.txt") + + +@skipIfMissingFeatures +class Sample(ut.TestCase): + system = sample.system + + def test_trajectory(self): + # compare the simulated and published trajectories + simul = sample.positions[:, :, 0:2] + paper = sample.data.reshape([-1, 3, 2]) + + for pid in range(3): + dist = [] + # the simulated trajectory is oversampled by a ratio of 60:1 + # compared to the published trajectory (60 +/- 5 to 1) + for desired in paper[:, pid]: + if desired[1] < y_min: + break + # find the closest point in the simulated trajectory + idx = np.abs(simul[:, pid, 1] - desired[1]).argmin() + actual = simul[idx, pid] + dist.append(np.linalg.norm(actual - desired)) + self.assertLess(idx, sample.intsteps, msg='Insufficient sampling') + np.testing.assert_allclose(dist, 0, rtol=0, atol=0.7) + + +if __name__ == "__main__": + ut.main()