From 8bd7b898919d8fa0c966634bcf314589f1cfa932 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 8 May 2023 20:15:00 +0200 Subject: [PATCH] walberla: Restore all LB samples --- samples/object_in_fluid/motivation.py | 36 +++++++++---------- testsuite/python/ek_diffusion.py | 5 ++- testsuite/python/tests_common.py | 27 ++++++-------- testsuite/scripts/samples/CMakeLists.txt | 2 +- .../test_object_in_fluid__motivation.py | 19 +++++++--- 5 files changed, 45 insertions(+), 44 deletions(-) diff --git a/samples/object_in_fluid/motivation.py b/samples/object_in_fluid/motivation.py index 44d554596dd..52f5316ae82 100644 --- a/samples/object_in_fluid/motivation.py +++ b/samples/object_in_fluid/motivation.py @@ -26,6 +26,7 @@ espressomd.assert_features(required_features) import os +import tqdm import argparse import warnings @@ -46,8 +47,8 @@ boxX = 22.0 boxY = 14.0 -boxZ = 15.0 -time_step = 0.1 +boxZ = 6.0 +time_step = 0.05 system = espressomd.System(box_l=(boxX, boxY, boxZ)) system.time_step = time_step @@ -63,18 +64,15 @@ # creating the RBCs cell0 = oif.OifCell(cell_type=cell_type, particle_type=0, origin=[5.0, 5.0, 3.0]) -cell1 = oif.OifCell(cell_type=cell_type, - particle_type=1, origin=[5.0, 5.0, 7.0]) # cell-wall interactions -system.non_bonded_inter[0, 10].soft_sphere.set_params( - a=0.0001, n=1.2, cutoff=0.1, offset=0.0) -system.non_bonded_inter[1, 10].soft_sphere.set_params( +system.non_bonded_inter[cell0.particle_type, 10].soft_sphere.set_params( a=0.0001, n=1.2, cutoff=0.1, offset=0.0) # fluid -lbf = espressomd.lb.LBFluidWalberla(agrid=1., density=1., kinematic_viscosity=1.5, - tau=0.1, ext_force_density=[0.002, 0., 0.]) +lbf = espressomd.lb.LBFluidWalberla( + agrid=1., density=1., kinematic_viscosity=1.5, tau=system.time_step, + ext_force_density=[0.025, 0., 0.], single_precision=True) system.actors.add(lbf) system.thermostat.set_lb(LB_fluid=lbf, gamma=1.5) @@ -140,15 +138,15 @@ system.constraints.add(shape=shape, particle_type=10) -maxCycle = 50 +def write_cells_vtk(i): + filepath = os.path.join(output_path, "cell{cell_id}_{index}.vtk") + cell0.output_vtk_pos_folded(file_name=filepath.format(cell_id=0, index=i)) + + +maxCycle = 100 # main integration loop -cell0.output_vtk_pos_folded(file_name=os.path.join(output_path, "cell0_0.vtk")) -cell1.output_vtk_pos_folded(file_name=os.path.join(output_path, "cell1_0.vtk")) -for i in range(1, maxCycle): - system.integrator.run(steps=500) - cell0.output_vtk_pos_folded( - file_name=os.path.join(output_path, f"cell0_{i}.vtk")) - cell1.output_vtk_pos_folded( - file_name=os.path.join(output_path, f"cell1_{i}.vtk")) - print(f"time: {i * time_step:.1f}") +for i in tqdm.tqdm(range(maxCycle)): + write_cells_vtk(i) + system.integrator.run(steps=100) +write_cells_vtk(maxCycle) print("Simulation completed.") diff --git a/testsuite/python/ek_diffusion.py b/testsuite/python/ek_diffusion.py index 9f2f4861d75..5ea77a6ea44 100644 --- a/testsuite/python/ek_diffusion.py +++ b/testsuite/python/ek_diffusion.py @@ -95,8 +95,7 @@ def detail_test_diffusion(self, single_precision: bool): # check that the density is conserved np.testing.assert_almost_equal( np.sum(simulated_density), self.DENSITY, decimal_precision) - if np.any(simulated_density < 0.): - self.fail("EK density array contains negative densities!") + assert np.all(simulated_density >= 0.), "EK has negative densities" # check that the maximum is in the right place peak = np.unravel_index( @@ -124,7 +123,7 @@ def detail_test_diffusion(self, single_precision: bool): calc_density, simulated_density, atol=1e-5, - rtol=0) + rtol=0.) if __name__ == "__main__": diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py index 23c7eb17850..86f1dff7ec6 100644 --- a/testsuite/python/tests_common.py +++ b/testsuite/python/tests_common.py @@ -18,6 +18,7 @@ # import pathlib +import itertools import numpy as np @@ -356,24 +357,16 @@ def fold_index(idx, shape): def get_lb_nodes_around_pos(pos, lbf): - """Returns lb node(s) relevant for interpolation around the given position""" - - pos_lb_units = pos / lbf.agrid - .5 # Rel to node centers + """Returns LB node(s) relevant for interpolation around the given position""" + pos_lb_units = pos / lbf.agrid - 0.5 # relative to node centers lower_left_index = np.array(np.floor(pos_lb_units), dtype=int) - # Position exactly on a node? Return just that node - if np.all(lower_left_index == pos_lb_units): - return [lbf[fold_index(lower_left_index, lbf.shape)]] - - # Otherwise return 8 surrounding nodes - res = [] - for i in (0, 1): - for j in (0, 1): - for k in (0, 1): - res.append(lbf[fold_index(lower_left_index + - np.array((i, j, k), dtype=int), lbf.shape)]) - return res + nodes = [] + for i, j, k in itertools.product([0, 1], repeat=3): + index = lower_left_index + np.array((i, j, k), dtype=int) + nodes.append(lbf[fold_index(index, lbf.shape)]) + return nodes def random_dipoles(n_particles): @@ -416,11 +409,11 @@ def check_non_bonded_loop_trace(ut_obj, system, cutoff=None): msg=f"Extra pair from core {p}") if (p[0], p[1]) in py_distances: np.testing.assert_allclose( - np.copy(p[4]), -np.copy(py_distances[p[0], p[1]])) + np.copy(p[4]), -py_distances[p[0], p[1]]) del py_distances[p[0], p[1]] elif (p[1], p[0]) in py_distances: np.testing.assert_allclose( - np.copy(p[4]), np.copy(py_distances[p[1], p[0]])) + np.copy(p[4]), py_distances[p[1], p[0]]) del py_distances[p[1], p[0]] for ids, dist in py_distances.items(): diff --git a/testsuite/scripts/samples/CMakeLists.txt b/testsuite/scripts/samples/CMakeLists.txt index 93bc6a462c1..78d4a9c2675 100644 --- a/testsuite/scripts/samples/CMakeLists.txt +++ b/testsuite/scripts/samples/CMakeLists.txt @@ -99,7 +99,7 @@ sample_test(FILE test_visualization_elc.py) sample_test(FILE test_visualization_npt.py) sample_test(FILE test_visualization_poiseuille.py) sample_test(FILE test_widom_insertion.py) -# sample_test(FILE test_object_in_fluid__motivation.py) # TODO WALBERLA +sample_test(FILE test_object_in_fluid__motivation.py) sample_test(FILE test_immersed_boundary.py) add_custom_target( diff --git a/testsuite/scripts/samples/test_object_in_fluid__motivation.py b/testsuite/scripts/samples/test_object_in_fluid__motivation.py index 8836577acff..a5b3777bfec 100644 --- a/testsuite/scripts/samples/test_object_in_fluid__motivation.py +++ b/testsuite/scripts/samples/test_object_in_fluid__motivation.py @@ -19,6 +19,7 @@ import importlib_wrapper import os import pathlib +import numpy as np os.chdir("@SAMPLES_DIR@/object_in_fluid") sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import( @@ -38,15 +39,25 @@ def test_file_generation(self): "wallBottom.vtk", "wallBack.vtk", "wallFront.vtk"] - for i in [0, 1]: - for j in range(sample.maxCycle): - basenames.append(f"cell{i}_{j}.vtk") + for j in range(sample.maxCycle + 1): + basenames.append(f"cell0_{j}.vtk") # test .vtk files exist path_vtk_root = pathlib.Path("output") for name in basenames: filepath = path_vtk_root / "sim0" / name - self.assertTrue(filepath.is_file(), f"File {filepath} not created") + self.assertTrue( + filepath.is_file(), + f"File '{filepath}' not created") + + # make sure we are still in the LB linear regime + lb_density = np.copy(sample.lbf[:, :, :].density) + np.testing.assert_allclose(lb_density, 1., atol=2e-3) + + # verify cell momentum + cell_vel = np.mean(self.system.part.all().v, axis=0) + np.testing.assert_allclose( + cell_vel, [1.54e-2, 1.8e-3, 0.], rtol=1e-2, atol=1e-6) if __name__ == "__main__":