From e0c1f27d6fbdbaa7928ce00a50ee520d0403e77b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 15 Aug 2023 17:18:27 +0200 Subject: [PATCH] Rewrite LB/EK interface * remove the actors list * enforce SOLID principles * add fine-grained event hooks --- doc/sphinx/advanced_methods.rst | 2 +- doc/sphinx/ek.rst | 61 +++-- doc/sphinx/integration.rst | 2 +- doc/sphinx/lb.rst | 21 +- .../active_matter/active_matter.ipynb | 2 +- .../electrokinetics/electrokinetics.ipynb | 30 +- .../lattice_boltzmann_poiseuille_flow.ipynb | 4 +- .../lattice_boltzmann_sedimentation.ipynb | 4 +- .../lattice_boltzmann_theory.ipynb | 2 +- doc/tutorials/polymers/polymers.ipynb | 4 +- .../raspberry_electrophoresis.ipynb | 2 +- maintainer/benchmarks/lb.py | 2 +- .../sampleImmersedBoundary.py | 3 +- samples/lb_circular_couette.py | 2 +- samples/lb_four_roller_mill.py | 2 +- samples/lb_planar_couette.py | 2 +- samples/lb_profile.py | 2 +- samples/lbf.py | 2 +- samples/object_in_fluid/motivation.py | 3 +- samples/visualization_lbboundaries.py | 2 +- samples/visualization_poiseuille.py | 2 +- src/core/CMakeLists.txt | 8 +- src/core/analysis/statistics.cpp | 7 +- .../CMakeLists.txt | 15 +- src/core/ek/EKNone.hpp | 39 +++ .../EKReactions.hpp | 28 +- src/core/ek/EKWalberla.cpp | 125 +++++++++ src/core/ek/EKWalberla.hpp | 77 ++++++ .../Implementation.hpp} | 36 +-- src/core/ek/Solver.cpp | 139 ++++++++++ src/core/ek/Solver.hpp | 101 +++++++ .../ek_reactions.cpp | 2 +- .../ek_reactions.hpp => ek/utils.hpp} | 21 +- src/core/event.cpp | 58 ++-- src/core/forces.cpp | 5 +- .../grid_based_algorithms/ek_container.cpp | 110 -------- .../grid_based_algorithms/lb_interface.cpp | 170 ------------ .../grid_based_algorithms/lb_interface.hpp | 131 --------- .../lb_interpolation.cpp | 60 ---- .../lb_walberla_instance.cpp | 104 ------- .../lb_walberla_instance.hpp | 58 ---- src/core/integrate.cpp | 99 +++++-- src/core/integrate.hpp | 18 ++ src/core/lb/CMakeLists.txt | 24 ++ .../Implementation.hpp} | 44 +-- src/core/lb/LBNone.hpp | 57 ++++ src/core/lb/LBWalberla.cpp | 89 ++++++ src/core/lb/LBWalberla.hpp | 89 ++++++ src/core/lb/Solver.cpp | 222 +++++++++++++++ src/core/lb/Solver.hpp | 160 +++++++++++ .../particle_coupling.cpp} | 94 +++---- .../particle_coupling.hpp} | 60 ++-- src/core/lb/utils.hpp | 30 ++ ...BFluxDensityProfileAtParticlePositions.cpp | 9 +- .../CylindricalLBVelocityProfile.cpp | 7 +- ...alLBVelocityProfileAtParticlePositions.cpp | 7 +- .../observables/LBFluidPressureTensor.hpp | 14 +- src/core/observables/LBVelocityProfile.cpp | 7 +- src/core/system/System.hpp | 6 +- src/core/system/System.impl.hpp | 3 + src/core/unit_tests/ek_interface_test.cpp | 192 +++++++++---- .../unit_tests/lb_particle_coupling_test.cpp | 153 +++++++---- .../VirtualSitesInertialessTracers.cpp | 41 +-- src/python/espressomd/actors.py | 99 ------- src/python/espressomd/electrokinetics.py | 41 +-- src/python/espressomd/system.py | 64 ++++- src/python/espressomd/thermostat.pxd | 5 +- src/python/espressomd/thermostat.pyx | 12 +- src/python/espressomd/visualization.py | 25 +- .../lees_edwards/LeesEdwards.hpp | 8 +- src/script_interface/walberla/EKContainer.hpp | 174 ++++++------ src/script_interface/walberla/EKFFT.hpp | 4 +- src/script_interface/walberla/EKReactions.hpp | 22 +- src/script_interface/walberla/LBFluid.cpp | 16 +- src/script_interface/walberla/LBFluid.hpp | 4 +- .../walberla/LatticeWalberla.hpp | 6 +- .../electrokinetics/EKContainer.hpp | 59 ++-- src/walberla_bridge/tests/CMakeLists.txt | 6 +- .../tests/ek_kernels_unit_tests.cpp | 167 ++++++++++++ ...it_tests.cpp => lb_kernels_unit_tests.cpp} | 18 +- testsuite/python/CMakeLists.txt | 1 - testsuite/python/actor.py | 258 ------------------ testsuite/python/array_properties.py | 5 +- testsuite/python/ek_boundary.py | 6 +- testsuite/python/ek_bulk_reactions.py | 13 +- testsuite/python/ek_diffusion.py | 8 +- testsuite/python/ek_eof.py | 10 +- testsuite/python/ek_fixeddensity.py | 6 +- testsuite/python/ek_fixedflux.py | 10 +- testsuite/python/ek_indexed_reactions.py | 14 +- testsuite/python/ek_interface.py | 73 +++-- testsuite/python/ek_noflux.py | 8 +- testsuite/python/engine_lb.py | 4 +- testsuite/python/ibm.py | 1 - testsuite/python/lattice_vtk.py | 10 +- testsuite/python/lb.py | 88 +++--- testsuite/python/lb_boundary.py | 4 +- testsuite/python/lb_boundary_velocity.py | 4 +- testsuite/python/lb_boundary_volume_force.py | 4 +- testsuite/python/lb_buoyancy_force.py | 4 +- testsuite/python/lb_circular_couette.py | 4 +- testsuite/python/lb_electrohydrodynamics.py | 4 +- testsuite/python/lb_interpolation.py | 4 +- testsuite/python/lb_lees_edwards.py | 47 ++-- .../lb_lees_edwards_particle_coupling.py | 2 +- testsuite/python/lb_mass_conservation.py | 4 +- testsuite/python/lb_momentum_conservation.py | 4 +- testsuite/python/lb_planar_couette.py | 4 +- testsuite/python/lb_poiseuille.py | 4 +- testsuite/python/lb_poiseuille_cylinder.py | 4 +- testsuite/python/lb_pressure_tensor.py | 11 +- testsuite/python/lb_shear.py | 4 +- testsuite/python/lb_slice.py | 4 +- testsuite/python/lb_stats.py | 4 +- testsuite/python/lb_stokes_sphere.py | 4 +- testsuite/python/lb_streaming.py | 4 +- testsuite/python/lb_thermo_virtual.py | 26 +- testsuite/python/lb_thermostat.py | 4 +- testsuite/python/linear_momentum_lb.py | 4 +- testsuite/python/long_range_actors.py | 2 +- testsuite/python/observable_cylindricalLB.py | 4 +- testsuite/python/observable_profileLB.py | 6 +- testsuite/python/save_checkpoint.py | 11 +- testsuite/python/test_checkpoint.py | 20 +- .../python/virtual_sites_tracers_common.py | 10 +- 125 files changed, 2498 insertions(+), 1837 deletions(-) rename src/core/{grid_based_algorithms => ek}/CMakeLists.txt (53%) create mode 100644 src/core/ek/EKNone.hpp rename src/core/{grid_based_algorithms => ek}/EKReactions.hpp (72%) create mode 100644 src/core/ek/EKWalberla.cpp create mode 100644 src/core/ek/EKWalberla.hpp rename src/core/{grid_based_algorithms/ek_container.hpp => ek/Implementation.hpp} (57%) create mode 100644 src/core/ek/Solver.cpp create mode 100644 src/core/ek/Solver.hpp rename src/core/{grid_based_algorithms => ek}/ek_reactions.cpp (97%) rename src/core/{grid_based_algorithms/ek_reactions.hpp => ek/utils.hpp} (66%) delete mode 100644 src/core/grid_based_algorithms/ek_container.cpp delete mode 100644 src/core/grid_based_algorithms/lb_interface.cpp delete mode 100644 src/core/grid_based_algorithms/lb_interface.hpp delete mode 100644 src/core/grid_based_algorithms/lb_interpolation.cpp delete mode 100644 src/core/grid_based_algorithms/lb_walberla_instance.cpp delete mode 100644 src/core/grid_based_algorithms/lb_walberla_instance.hpp create mode 100644 src/core/lb/CMakeLists.txt rename src/core/{grid_based_algorithms/lb_interpolation.hpp => lb/Implementation.hpp} (51%) create mode 100644 src/core/lb/LBNone.hpp create mode 100644 src/core/lb/LBWalberla.cpp create mode 100644 src/core/lb/LBWalberla.hpp create mode 100644 src/core/lb/Solver.cpp create mode 100644 src/core/lb/Solver.hpp rename src/core/{grid_based_algorithms/lb_particle_coupling.cpp => lb/particle_coupling.cpp} (79%) rename src/core/{grid_based_algorithms/lb_particle_coupling.hpp => lb/particle_coupling.hpp} (77%) create mode 100644 src/core/lb/utils.hpp delete mode 100644 src/python/espressomd/actors.py create mode 100644 src/walberla_bridge/tests/ek_kernels_unit_tests.cpp rename src/walberla_bridge/tests/{kernels_unit_tests.cpp => lb_kernels_unit_tests.cpp} (92%) delete mode 100644 testsuite/python/actor.py diff --git a/doc/sphinx/advanced_methods.rst b/doc/sphinx/advanced_methods.rst index 7676d831980..a2ccc248511 100644 --- a/doc/sphinx/advanced_methods.rst +++ b/doc/sphinx/advanced_methods.rst @@ -460,7 +460,7 @@ Specification of fluid and movement lbf = espressomd.lb.LBFluidWalberla(agrid=1, density=1.0, kinematic_viscosity=1.5, tau=time_step, ext_force_density=[0.002, 0.0, 0.0]) - system.actors.add(lbf) + self.system.lb = lbf This part of the script specifies the fluid that will get the system moving. Here ``agrid`` :math:`=\Delta x` is the spatial discretisation diff --git a/doc/sphinx/ek.rst b/doc/sphinx/ek.rst index a2e52502eda..b5fba56ed88 100644 --- a/doc/sphinx/ek.rst +++ b/doc/sphinx/ek.rst @@ -161,19 +161,23 @@ Here is a minimal working example:: system.time_step = 0.01 system.cell_system.skin = 1.0 - ek_lattice = espressomd.electrokinetics.LatticeWalberla(agrid=0.5, n_ghost_layers=1) - ek_solver = espressomd.electrokinetics.EKNone(lattice=ek_lattice) - system.ekcontainer.solver = ek_solver - system.ekcontainer.tau = system.time_step + lattice = espressomd.electrokinetics.LatticeWalberla(agrid=0.5, n_ghost_layers=1) + ek_solver = espressomd.electrokinetics.EKNone(lattice=lattice) + system.ekcontainer = espressomd.electrokinetics.EKContainer( + solver=ek_solver, tau=system.time_step) where ``system.ekcontainer`` is the EK system, ``ek_solver`` is the Poisson solver (here ``EKNone`` doesn't actually solve the electrostatic field, but -instead imposes a zero field), and ``ek_lattice`` contains the grid parameters. +instead imposes a zero field), and ``lattice`` contains the grid parameters. In this setup, the EK system doesn't contain any species. The following sections will show how to add species that can diffuse, advect, react and/or electrostatically interact. An EK system can be set up at the same time as a LB system. +To detach an EK system, use the following syntax:: + + system.ekcontainer = None + .. _Diffusive species: Diffusive species @@ -181,9 +185,10 @@ Diffusive species :: ek_species = espressomd.electrokinetics.EKSpecies( - lattice=ek_lattice, + lattice=lattice, single_precision=False, kT=1.0, + tau=system.time_step, density=0.85, valency=0.0, diffusion=0.1, @@ -234,8 +239,8 @@ Checkpointing :: - ek.save_checkpoint(path, binary) - ek.load_checkpoint(path, binary) + ek_species.save_checkpoint(path, binary) + ek_species.ekcontainer.load_checkpoint(path, binary) The first command saves all of the EK nodes' properties to an ASCII (``binary=False``) or binary (``binary=True``) format respectively. @@ -260,7 +265,7 @@ one or multiple fluid field data into a single file using # create a VTK callback that automatically writes every 10 EK steps ek_vtk = espressomd.electrokinetics.VTKOutput( identifier="ek_vtk_automatic", observables=vtk_obs, delta_N=10) - ek.add_vtk_writer(vtk=ek_vtk) + ek_species.add_vtk_writer(vtk=ek_vtk) system.integrator.run(100) # can be deactivated ek_vtk.disable() @@ -269,7 +274,7 @@ one or multiple fluid field data into a single file using # create a VTK callback that writes only when explicitly called ek_vtk_on_demand = espressomd.electrokinetics.VTKOutput( identifier="ek_vtk_now", observables=vtk_obs) - ek.add_vtk_writer(vtk=ek_vtk_on_demand) + ek_species.add_vtk_writer(vtk=ek_vtk_on_demand) ek_vtk_on_demand.write() Currently only supports the species density. @@ -306,12 +311,14 @@ is available through :class:`~espressomd.io.vtk.VTKReader`:: system.cell_system.skin = 0.4 system.time_step = 0.1 - lattice = espressomd.electrokinetics.LatticeWalberla(agrid=1.) - species = espressomd.electrokinetics.EKSpecies( - lattice=lattice, density=1., kT=1., diffusion=0.1, valency=0., - advection=False, friction_coupling=False, tau=system.time_step) - system.ekcontainer.tau = species.tau - system.ekcontainer.add(species) + lattice = espressomd.electrokinetics.LatticeWalberla(agrid=1., n_ghost_layers=1) + ek_solver = espressomd.electrokinetics.EKNone(lattice=lattice) + ek_species = espressomd.electrokinetics.EKSpecies( + lattice=lattice, density=1., kT=1., diffusion=0.1, valency=0., + advection=False, friction_coupling=False, tau=system.time_step) + system.ekcontainer = espressomd.electrokinetics.EKContainer( + solver=ek_solver, tau=ek_species.tau) + system.ekcontainer.add(ek_species) system.integrator.run(10) vtk_reader = espressomd.io.vtk.VTKReader() @@ -327,7 +334,7 @@ is available through :class:`~espressomd.io.vtk.VTKReader`:: identifier=label_vtk, delta_N=0, observables=["density"], base_folder=str(path_vtk_root)) - species.add_vtk_writer(vtk=ek_vtk) + ek_species.add_vtk_writer(vtk=ek_vtk) ek_vtk.write() # read VTK file @@ -335,7 +342,7 @@ is available through :class:`~espressomd.io.vtk.VTKReader`:: vtk_density = vtk_grids[label_density] # check VTK values match node values - ek_density = np.copy(lbf[:, :, :].density) + ek_density = np.copy(ek_species[:, :, :].density) np.testing.assert_allclose(vtk_density, ek_density, rtol=1e-10, atol=0.) .. _Setting up EK boundary conditions: @@ -361,17 +368,19 @@ One can set (or update) the boundary conditions of individual nodes:: system.cell_system.skin = 0.1 system.time_step = 0.01 lattice = espressomd.electrokinetics.LatticeWalberla(agrid=0.5, n_ghost_layers=1) + ek_solver = espressomd.electrokinetics.EKNone(lattice=lattice) ek_species = espressomd.electrokinetics.EKSpecies( - kT=1.5, lattice=self.lattice, density=0.85, valency=0., diffusion=0.1, + kT=1.5, lattice=lattice, density=0.85, valency=0., diffusion=0.1, advection=False, friction_coupling=False, tau=system.time_step) - system.ekcontainer.tau = species.tau + system.ekcontainer = espressomd.electrokinetics.EKContainer( + solver=ek_solver, tau=ek_species.tau) system.ekcontainer.add(ek_species) # set node fixed density boundary conditions - lbf[0, 0, 0].boundary = espressomd.electrokinetics.DensityBoundary(1.) + ek_species[0, 0, 0].boundary = espressomd.electrokinetics.DensityBoundary(1.) # update node fixed density boundary conditions - lbf[0, 0, 0].boundary = espressomd.electrokinetics.DensityBoundary(2.) + ek_species[0, 0, 0].boundary = espressomd.electrokinetics.DensityBoundary(2.) # remove node boundary conditions - lbf[0, 0, 0].boundary = None + ek_species[0, 0, 0].boundary = None .. _Shape-based EK boundary conditions: @@ -387,10 +396,12 @@ Adding a shape-based boundary is straightforward:: system.cell_system.skin = 0.1 system.time_step = 0.01 lattice = espressomd.electrokinetics.LatticeWalberla(agrid=0.5, n_ghost_layers=1) + ek_solver = espressomd.electrokinetics.EKNone(lattice=lattice) ek_species = espressomd.electrokinetics.EKSpecies( - kT=1.5, lattice=self.lattice, density=0.85, valency=0.0, diffusion=0.1, + kT=1.5, lattice=lattice, density=0.85, valency=0.0, diffusion=0.1, advection=False, friction_coupling=False, tau=system.time_step) - system.ekcontainer.tau = species.tau + system.ekcontainer = espressomd.electrokinetics.EKContainer( + solver=ek_solver, tau=ek_species.tau) system.ekcontainer.add(ek_species) # set fixed density boundary conditions wall = espressomd.shapes.Wall(normal=[1., 0., 0.], dist=2.5) diff --git a/doc/sphinx/integration.rst b/doc/sphinx/integration.rst index b7f0bd34eaf..cf521a64d6e 100644 --- a/doc/sphinx/integration.rst +++ b/doc/sphinx/integration.rst @@ -659,7 +659,7 @@ parameter ``gamma``. To enable the LB thermostat, use:: import espressomd.lb system = espressomd.System(box_l=[1, 1, 1]) lbf = espressomd.lb.LBFluidWalberla(agrid=1, density=1, kinematic_viscosity=1, tau=0.01) - system.actors.add(lbf) + self.system.lb = lbf system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1.5) No other thermostatting mechanism is necessary diff --git a/doc/sphinx/lb.rst b/doc/sphinx/lb.rst index d23e928e037..30828cd6055 100644 --- a/doc/sphinx/lb.rst +++ b/doc/sphinx/lb.rst @@ -41,13 +41,13 @@ The following minimal example illustrates how to use the LBM in |es|:: system = espressomd.System(box_l=[10, 20, 30]) system.time_step = 0.01 system.cell_system.skin = 0.4 - lb = espressomd.lb.LBFluidWalberla(agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.01) - system.actors.add(lb) + lbf = espressomd.lb.LBFluidWalberla(agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.01) + system.lb = lbf system.integrator.run(100) To use the GPU-accelerated variant, replace line 6 in the example above by:: - lb = espressomd.lb.LBFluidWalberlaGPU(agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.01) + lbf = espressomd.lb.LBFluidWalberlaGPU(agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.01) .. note:: Feature ``CUDA`` required for the GPU-accelerated variant @@ -99,6 +99,11 @@ units to be applied to the fluid. Before running a simulation at least the following parameters must be set up: ``agrid``, ``tau``, ``kinematic_viscosity``, ``density``. +To detach the LBM solver, use this syntax:: + + system.lb = None + + Performance considerations ^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -321,7 +326,7 @@ is available through :class:`~espressomd.io.vtk.VTKReader`:: lbf = espressomd.lb.LBFluidWalberla( agrid=1., tau=0.1, density=1., kinematic_viscosity=1.) - system.actors.add(lbf) + system.lb = lbf system.integrator.run(10) vtk_reader = espressomd.io.vtk.VTKReader() @@ -378,8 +383,8 @@ of the LBM in analogy to the example for the CPU given in section system = espressomd.System(box_l=[10, 20, 30]) system.time_step = 0.01 system.cell_system.skin = 0.4 - lb = espressomd.lb.LBFluidWalberlaGPU(agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.01) - system.actors.add(lb) + lbf = espressomd.lb.LBFluidWalberlaGPU(agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.01) + system.lb = lbf system.integrator.run(100) .. _Electrohydrodynamics: @@ -429,7 +434,7 @@ One can set (or update) the slip velocity of individual nodes:: system.cell_system.skin = 0.1 system.time_step = 0.01 lbf = espressomd.lb.LBFluidWalberla(agrid=0.5, density=1.0, kinematic_viscosity=1.0, tau=0.01) - system.actors.add(lbf) + system.lb = lbf # make one node a boundary node with a slip velocity lbf[0, 0, 0].boundary = espressomd.lb.VelocityBounceBack([0, 0, 1]) # update node for no-slip boundary conditions @@ -450,7 +455,7 @@ Adding a shape-based boundary is straightforward:: system.cell_system.skin = 0.1 system.time_step = 0.01 lbf = espressomd.lb.LBFluidWalberla(agrid=0.5, density=1.0, kinematic_viscosity=1.0, tau=0.01) - system.actors.add(lbf) + system.lb = lbf # set up shear flow between two sliding walls wall1 = espressomd.shapes.Wall(normal=[+1., 0., 0.], dist=2.5) lbf.add_boundary_from_shape(shape=wall1, velocity=[0., +0.05, 0.]) diff --git a/doc/tutorials/active_matter/active_matter.ipynb b/doc/tutorials/active_matter/active_matter.ipynb index 08046e92401..3bdab0431dd 100644 --- a/doc/tutorials/active_matter/active_matter.ipynb +++ b/doc/tutorials/active_matter/active_matter.ipynb @@ -903,7 +903,7 @@ " density=HYDRO_PARAMS['dens'],\n", " kinematic_viscosity=HYDRO_PARAMS['visc'],\n", " tau=HYDRO_PARAMS['time_step'])\n", - "system.actors.add(lbf)\n", + "system.lb = lbf\n", "system.thermostat.set_lb(LB_fluid=lbf, gamma=HYDRO_PARAMS['gamma'], seed=42)\n", "```" ] diff --git a/doc/tutorials/electrokinetics/electrokinetics.ipynb b/doc/tutorials/electrokinetics/electrokinetics.ipynb index a20ec75b40e..cde7cc5e494 100644 --- a/doc/tutorials/electrokinetics/electrokinetics.ipynb +++ b/doc/tutorials/electrokinetics/electrokinetics.ipynb @@ -354,6 +354,7 @@ "source": [ "# Initializing espresso modules and the numpy package\n", "import espressomd\n", + "import espressomd.lb\n", "import espressomd.electrokinetics\n", "import espressomd.shapes\n", "\n", @@ -451,8 +452,10 @@ "outputs": [], "source": [ "viscosity_kinematic = viscosity_dynamic / density_water\n", - "lbf = espressomd.lb.LBFluidWalberla(lattice=lattice, density=density_water, kinematic_viscosity=viscosity_kinematic, tau=dt, single_precision=single_precision)\n", - "system.actors.add(lbf)" + "lbf = espressomd.lb.LBFluidWalberla(lattice=lattice, density=density_water,\n", + " kinematic_viscosity=viscosity_kinematic,\n", + " tau=dt, single_precision=single_precision)\n", + "system.lb = lbf" ] }, { @@ -461,10 +464,11 @@ "metadata": {}, "outputs": [], "source": [ - "eksolver = espressomd.electrokinetics.EKFFT(lattice=lattice, permittivity=permittivity, single_precision=single_precision)\n", + "eksolver = espressomd.electrokinetics.EKFFT(lattice=lattice,\n", + " permittivity=permittivity,\n", + " single_precision=single_precision)\n", "\n", - "system.ekcontainer.solver = eksolver\n", - "system.ekcontainer.tau = dt" + "system.ekcontainer = espressomd.electrokinetics.EKContainer(solver=eksolver, tau=dt)" ] }, { @@ -474,7 +478,10 @@ "outputs": [], "source": [ "density_counterions = -2.0 * sigma / width\n", - "ekspecies = espressomd.electrokinetics.EKSpecies(lattice=lattice, density=0.0, kT=kT, diffusion=D, valency=valency, advection=True, friction_coupling=True, ext_efield=ext_force_density, single_precision=single_precision, tau=dt)\n", + "ekspecies = espressomd.electrokinetics.EKSpecies(\n", + " lattice=lattice, density=0.0, kT=kT, diffusion=D, valency=valency,\n", + " advection=True, friction_coupling=True, ext_efield=ext_force_density,\n", + " single_precision=single_precision, tau=dt)\n", "system.ekcontainer.add(ekspecies)" ] }, @@ -484,7 +491,10 @@ "metadata": {}, "outputs": [], "source": [ - "ekwallcharge = espressomd.electrokinetics.EKSpecies(lattice=lattice, density=0.0, kT=kT, diffusion=0., valency=-valency, advection=False, friction_coupling=False, ext_efield=[0, 0, 0], single_precision=single_precision, tau=dt)\n", + "ekwallcharge = espressomd.electrokinetics.EKSpecies(\n", + " lattice=lattice, density=0.0, kT=kT, diffusion=0., valency=-valency,\n", + " advection=False, friction_coupling=False, ext_efield=[0, 0, 0],\n", + " single_precision=single_precision, tau=dt)\n", "system.ekcontainer.add(ekwallcharge)" ] }, @@ -527,8 +537,10 @@ "outputs": [], "source": [ "for shape_obj in (wall_left, wall_right):\n", - " ekspecies.add_boundary_from_shape(shape=shape_obj, value=[0., 0., 0.], boundary_type=espressomd.electrokinetics.FluxBoundary)\n", - " ekspecies.add_boundary_from_shape(shape=shape_obj, value=0.0, boundary_type=espressomd.electrokinetics.DensityBoundary)\n", + " ekspecies.add_boundary_from_shape(shape=shape_obj, value=[0., 0., 0.],\n", + " boundary_type=espressomd.electrokinetics.FluxBoundary)\n", + " ekspecies.add_boundary_from_shape(shape=shape_obj, value=0.0,\n", + " boundary_type=espressomd.electrokinetics.DensityBoundary)\n", " lbf.add_boundary_from_shape(shape=shape_obj, velocity=[0., 0., 0.])" ] }, diff --git a/doc/tutorials/lattice_boltzmann/lattice_boltzmann_poiseuille_flow.ipynb b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_poiseuille_flow.ipynb index b0d51f19add..90eca9bf034 100644 --- a/doc/tutorials/lattice_boltzmann/lattice_boltzmann_poiseuille_flow.ipynb +++ b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_poiseuille_flow.ipynb @@ -105,7 +105,7 @@ "solution2_first": true }, "source": [ - "Create a lattice-Boltzmann actor and append it to the list of system actors. Use the GPU implementation of LB.\n", + "Create a lattice-Boltzmann actor and attach it to the system. Use the GPU implementation of LB.\n", "\n", "You can refer to section [setting up a LB fluid](https://espressomd.github.io/doc/lb.html#setting-up-a-lb-fluid)\n", "in the user guide." @@ -123,7 +123,7 @@ " kinematic_viscosity=VISCOSITY,\n", " tau=TIME_STEP,\n", " ext_force_density=FORCE_DENSITY)\n", - "system.actors.add(lbf)\n", + "system.lb = lbf\n", "```" ] }, diff --git a/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb index a8100e9cbdb..bbc442c04d9 100644 --- a/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb +++ b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb @@ -329,7 +329,7 @@ "\n", "**Exercise:**\n", "* create an unthermalized lattice-Boltzmann object with viscosity 1, density 1, grid size equal to\n", - " ``spacing`` and LB time step equal to the MD time step and add this object to the system list of actors\n", + " ``spacing`` and LB time step equal to the MD time step and add this object to the system\n", " ([user guide](https://espressomd.github.io/doc/lb.html#setting-up-a-lb-fluid))\n", "* activate particle coupling to the fluid by setting the LB thermostat with $\\gamma=15$\n", " ([user guide](https://espressomd.github.io/doc/lb.html#coupling-lb-to-a-md-simulation))" @@ -346,7 +346,7 @@ " density=1.,\n", " kinematic_viscosity=1.,\n", " tau=system.time_step, kT=0.)\n", - "system.actors.add(lbf)\n", + "system.lb = lbf\n", "system.thermostat.set_lb(LB_fluid=lbf, gamma=15., seed=0)\n", "```" ] diff --git a/doc/tutorials/lattice_boltzmann/lattice_boltzmann_theory.ipynb b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_theory.ipynb index c4ffcef9d42..bb177147687 100644 --- a/doc/tutorials/lattice_boltzmann/lattice_boltzmann_theory.ipynb +++ b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_theory.ipynb @@ -273,7 +273,7 @@ "# Initialize an LBFluidWalberla with the minimum set of valid parameters.\n", "lbf = espressomd.lb.LBFluidWalberla(agrid=1, density=10, kinematic_viscosity=.1, tau=0.01)\n", "# Activate the LB by adding it to the System's actor list.\n", - "system.actors.add(lbf)\n", + "system.lb = lbf\n", "```" ] }, diff --git a/doc/tutorials/polymers/polymers.ipynb b/doc/tutorials/polymers/polymers.ipynb index 40d4d3bbaa2..738a8fff77c 100644 --- a/doc/tutorials/polymers/polymers.ipynb +++ b/doc/tutorials/polymers/polymers.ipynb @@ -300,7 +300,7 @@ " lbf = espressomd.lb.LBFluidWalberla(kT=kT, seed=42, agrid=1, density=1,\n", " kinematic_viscosity=5, tau=system.time_step,\n", " single_precision=True)\n", - " system.actors.add(lbf)\n", + " system.lb = lbf\n", " system.thermostat.set_lb(LB_fluid=lbf, gamma=gamma, seed=42)" ] }, @@ -440,7 +440,7 @@ "\n", " # reset system\n", " system.part.clear()\n", - " system.actors.clear()\n", + " system.lb = None\n", " system.thermostat.turn_off()\n", " system.auto_update_accumulators.clear()\n", "\n", diff --git a/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb b/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb index 03bbfd90c3f..3e71b0771ea 100644 --- a/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb +++ b/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb @@ -675,7 +675,7 @@ "outputs": [], "source": [ "system.thermostat.turn_off()\n", - "system.actors.add(lb)\n", + "system.lb = lb\n", "system.thermostat.set_lb(LB_fluid=lb, seed=123, gamma=20.0)" ] }, diff --git a/maintainer/benchmarks/lb.py b/maintainer/benchmarks/lb.py index 42f960d4e15..1815e444437 100644 --- a/maintainer/benchmarks/lb.py +++ b/maintainer/benchmarks/lb.py @@ -138,7 +138,7 @@ lb_class = espressomd.lb.LBFluidWalberlaGPU lbf = lb_class(agrid=agrid, tau=system.time_step, kinematic_viscosity=1., density=1., single_precision=args.single_precision) -system.actors.add(lbf) +system.lb = lbf # time integration loop diff --git a/samples/immersed_boundary/sampleImmersedBoundary.py b/samples/immersed_boundary/sampleImmersedBoundary.py index 11b63a45753..09fa59b2ef5 100644 --- a/samples/immersed_boundary/sampleImmersedBoundary.py +++ b/samples/immersed_boundary/sampleImmersedBoundary.py @@ -78,8 +78,7 @@ lbf = espressomd.lb.LBFluidWalberla( agrid=1, density=1, kinematic_viscosity=1, tau=system.time_step, ext_force_density=[force, 0, 0]) -system.actors.add(lbf) - +system.lb = lbf system.thermostat.set_lb(LB_fluid=lbf, gamma=1.0, act_on_virtual=False) # Setup boundaries diff --git a/samples/lb_circular_couette.py b/samples/lb_circular_couette.py index 9b30fbac048..0d915237ec6 100644 --- a/samples/lb_circular_couette.py +++ b/samples/lb_circular_couette.py @@ -51,7 +51,7 @@ system.periodicity = [False, False, True] lb_fluid = espressomd.lb.LBFluidWalberla( agrid=agrid, density=0.5, kinematic_viscosity=3.2, tau=system.time_step) -system.actors.add(lb_fluid) +system.lb = lb_fluid # set up cylinders cyl_center = agrid * (grid_size // 2 + 0.5) * [1, 1, 0] diff --git a/samples/lb_four_roller_mill.py b/samples/lb_four_roller_mill.py index dc7c33cbc23..c6cda80dff3 100644 --- a/samples/lb_four_roller_mill.py +++ b/samples/lb_four_roller_mill.py @@ -55,7 +55,7 @@ system.cell_system.skin = 0.1 lb_fluid = espressomd.lb.LBFluidWalberla( agrid=agrid, density=0.5, kinematic_viscosity=3.2, tau=system.time_step) -system.actors.add(lb_fluid) +system.lb = lb_fluid # set up rollers by adding tangential slip velocities to cylinders logging.info('Setting up the rollers') diff --git a/samples/lb_planar_couette.py b/samples/lb_planar_couette.py index 8ffbc8f0680..f65a037f11c 100644 --- a/samples/lb_planar_couette.py +++ b/samples/lb_planar_couette.py @@ -78,7 +78,7 @@ def analytical(x, t, nu, v, h, k_max): lbf = espressomd.lb.LBFluidWalberla( agrid=1., density=1., kinematic_viscosity=nu, tau=1.) -system.actors.add(lbf) +system.lb = lbf # sampling time_breakpoints = [50, 200, 500, 2000] diff --git a/samples/lb_profile.py b/samples/lb_profile.py index 6043a649fce..763bd7f95a0 100644 --- a/samples/lb_profile.py +++ b/samples/lb_profile.py @@ -44,7 +44,7 @@ lb_fluid = espressomd.lb.LBFluidWalberla( agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.01, ext_force_density=[0, 0, 0.15], kT=0.0) -system.actors.add(lb_fluid) +system.lb = lb_fluid system.thermostat.set_lb(LB_fluid=lb_fluid, seed=23) ctp = espressomd.math.CylindricalTransformationParameters( center=[5.0, 5.0, 0.0], diff --git a/samples/lbf.py b/samples/lbf.py index 643fb172469..b418df1a960 100644 --- a/samples/lbf.py +++ b/samples/lbf.py @@ -66,7 +66,7 @@ lbf = espressomd.lb.LBFluidWalberlaGPU(**lb_params) else: lbf = espressomd.lb.LBFluidWalberla(**lb_params) -system.actors.add(lbf) +system.lb = lbf system.thermostat.set_lb(LB_fluid=lbf, gamma=1.0) print(lbf.get_params()) diff --git a/samples/object_in_fluid/motivation.py b/samples/object_in_fluid/motivation.py index 52f5316ae82..5e5a4e2d076 100644 --- a/samples/object_in_fluid/motivation.py +++ b/samples/object_in_fluid/motivation.py @@ -20,6 +20,7 @@ """ import espressomd +import espressomd.lb import espressomd.shapes required_features = ["WALBERLA", "EXTERNAL_FORCES", "SOFT_SPHERE", "MASS"] @@ -73,7 +74,7 @@ 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.lb = lbf system.thermostat.set_lb(LB_fluid=lbf, gamma=1.5) # creating boundaries and obstacles in the channel diff --git a/samples/visualization_lbboundaries.py b/samples/visualization_lbboundaries.py index 92ca3b59b7f..c0fea719076 100644 --- a/samples/visualization_lbboundaries.py +++ b/samples/visualization_lbboundaries.py @@ -35,7 +35,7 @@ lb_fluid = espressomd.lb.LBFluidWalberla( agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.01, ext_force_density=[0, 0, 0.15]) -system.actors.add(lb_fluid) +system.lb = lb_fluid cylinder_shape = espressomd.shapes.Cylinder( center=[5.0, 5.0, 5.0], diff --git a/samples/visualization_poiseuille.py b/samples/visualization_poiseuille.py index 94c7f614cc2..809c3f63238 100644 --- a/samples/visualization_poiseuille.py +++ b/samples/visualization_poiseuille.py @@ -55,7 +55,7 @@ lbf = espressomd.lb.LBFluidWalberla(kT=0, agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.1, ext_force_density=[0, 0.003, 0]) -system.actors.add(lbf) +system.lb = lbf system.thermostat.set_lb(LB_fluid=lbf, gamma=1.5) # Setup boundaries diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 45d3942899d..6aa9c4aa872 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -95,6 +95,11 @@ target_link_libraries( target_include_directories(espresso_core PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) +if(ESPRESSO_BUILD_WITH_WALBERLA) + target_link_libraries(espresso_core PRIVATE espresso::walberla + ${WALBERLA_LIBS}) +endif() + add_subdirectory(accumulators) add_subdirectory(analysis) add_subdirectory(bond_breakage) @@ -102,13 +107,14 @@ add_subdirectory(bonded_interactions) add_subdirectory(cell_system) add_subdirectory(cluster_analysis) add_subdirectory(constraints) +add_subdirectory(ek) add_subdirectory(electrostatics) add_subdirectory(error_handling) add_subdirectory(galilei) -add_subdirectory(grid_based_algorithms) add_subdirectory(immersed_boundary) add_subdirectory(integrators) add_subdirectory(io) +add_subdirectory(lb) add_subdirectory(magnetostatics) add_subdirectory(nonbonded_interactions) add_subdirectory(object-in-fluid) diff --git a/src/core/analysis/statistics.cpp b/src/core/analysis/statistics.cpp index b48c0619676..c8effaf24d9 100644 --- a/src/core/analysis/statistics.cpp +++ b/src/core/analysis/statistics.cpp @@ -30,8 +30,8 @@ #include "cells.hpp" #include "errorhandling.hpp" #include "grid.hpp" -#include "grid_based_algorithms/lb_interface.hpp" #include "partCfg_global.hpp" +#include "system/System.hpp" #include #include @@ -83,8 +83,9 @@ Utils::Vector3d calc_linear_momentum(bool include_particles, return m + p.mass() * p.v(); }); } - if (include_lbfluid and lattice_switch != ActiveLB::NONE) { - momentum += LB::calc_fluid_momentum() * LB::get_lattice_speed(); + auto const &system = System::get_system(); + if (include_lbfluid and system.lb.is_solver_set()) { + momentum += system.lb.get_momentum() * system.lb.get_lattice_speed(); } return momentum; } diff --git a/src/core/grid_based_algorithms/CMakeLists.txt b/src/core/ek/CMakeLists.txt similarity index 53% rename from src/core/grid_based_algorithms/CMakeLists.txt rename to src/core/ek/CMakeLists.txt index 1f9c29d3a76..2bf9fab665d 100644 --- a/src/core/grid_based_algorithms/CMakeLists.txt +++ b/src/core/ek/CMakeLists.txt @@ -1,5 +1,5 @@ # -# Copyright (C) 2018-2022 The ESPResSo project +# Copyright (C) 2023 The ESPResSo project # # This file is part of ESPResSo. # @@ -17,17 +17,8 @@ # along with this program. If not, see . # -target_sources( - espresso_core - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/lb_interface.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/lb_interpolation.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/lb_particle_coupling.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/ek_container.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/ek_reactions.cpp) +target_sources(espresso_core PRIVATE Solver.cpp) if(ESPRESSO_BUILD_WITH_WALBERLA) - target_link_libraries(espresso_core PRIVATE espresso::walberla - ${WALBERLA_LIBS}) - target_sources(espresso_core - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/lb_walberla_instance.cpp) + target_sources(espresso_core PRIVATE EKWalberla.cpp) endif() diff --git a/src/core/ek/EKNone.hpp b/src/core/ek/EKNone.hpp new file mode 100644 index 00000000000..450822216dd --- /dev/null +++ b/src/core/ek/EKNone.hpp @@ -0,0 +1,39 @@ +/* + * Copyright (C) 2023 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 . + */ + +#pragma once + +#include "utils.hpp" + +namespace EK { + +struct EKNone { + bool is_ready_for_propagation() const { throw NoEKActive{}; } + void propagate() { throw NoEKActive{}; } + double get_tau() const { throw NoEKActive{}; } + void veto_time_step(double) const { throw NoEKActive{}; } + void sanity_checks() const { throw NoEKActive{}; } + void on_cell_structure_change() const { throw NoEKActive{}; } + void on_boxl_change() const { throw NoEKActive{}; } + void on_node_grid_change() const { throw NoEKActive{}; } + void on_timestep_change() const { throw NoEKActive{}; } + void on_temperature_change() const { throw NoEKActive{}; } +}; + +} // namespace EK diff --git a/src/core/grid_based_algorithms/EKReactions.hpp b/src/core/ek/EKReactions.hpp similarity index 72% rename from src/core/grid_based_algorithms/EKReactions.hpp rename to src/core/ek/EKReactions.hpp index b1b7905e920..12d8869ca50 100644 --- a/src/core/grid_based_algorithms/EKReactions.hpp +++ b/src/core/ek/EKReactions.hpp @@ -17,13 +17,15 @@ * along with this program. If not, see . */ -#ifndef ESPRESSO_EKREACTIONS_HPP -#define ESPRESSO_EKREACTIONS_HPP +#pragma once #include +#include #include #include +namespace EK { + template class EKReactions { using container_type = std::vector>; @@ -36,18 +38,16 @@ template class EKReactions { container_type m_ekreactions; public: - void add(std::shared_ptr const &c) { - assert(std::find(m_ekreactions.begin(), m_ekreactions.end(), c) == - m_ekreactions.end()); - - m_ekreactions.emplace_back(c); + bool contains(std::shared_ptr const &ek_reaction) const noexcept { + return std::find(begin(), end(), ek_reaction) != end(); + } + void add(std::shared_ptr const &ek_reaction) { + assert(not contains(ek_reaction)); + m_ekreactions.emplace_back(ek_reaction); } - void remove(std::shared_ptr const &c) { - assert(std::find(m_ekreactions.begin(), m_ekreactions.end(), c) != - m_ekreactions.end()); - m_ekreactions.erase( - std::remove(m_ekreactions.begin(), m_ekreactions.end(), c), - m_ekreactions.end()); + void remove(std::shared_ptr const &ek_reaction) { + assert(contains(ek_reaction)); + m_ekreactions.erase(std::remove(begin(), end(), ek_reaction), end()); } iterator begin() { return m_ekreactions.begin(); } @@ -57,4 +57,4 @@ template class EKReactions { [[nodiscard]] bool empty() const { return m_ekreactions.empty(); } }; -#endif +} // namespace EK diff --git a/src/core/ek/EKWalberla.cpp b/src/core/ek/EKWalberla.cpp new file mode 100644 index 00000000000..2a23fd24348 --- /dev/null +++ b/src/core/ek/EKWalberla.cpp @@ -0,0 +1,125 @@ +/* + * Copyright (C) 2022-2023 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 . + */ + +#include "config/config.hpp" + +#ifdef WALBERLA + +#include "ek/EKReactions.hpp" +#include "ek/EKWalberla.hpp" +#include "errorhandling.hpp" +#include "grid.hpp" +#include "integrate.hpp" +#include "lb/Implementation.hpp" +#include "lb/LBWalberla.hpp" +#include "system/System.hpp" + +#include +#include +#include +#include + +#include +#include + +namespace EK { + +double EKWalberla::get_tau() const { return ek_container->get_tau(); } + +bool EKWalberla::is_ready_for_propagation() const noexcept { + return not ek_container->empty(); +} + +struct FieldsConnector { + std::size_t velocity_field_id{}; + std::size_t force_field_id{}; + void operator()(LB::Solver::Implementation const &impl) { + using lb_value_type = std::shared_ptr; + if (impl.solver.has_value()) { + if (auto const *ptr = std::get_if(&(*impl.solver))) { + auto const &instance = **ptr; + velocity_field_id = instance.lb_fluid->get_velocity_field_id(); + force_field_id = instance.lb_fluid->get_force_field_id(); + } + } + } +}; + +void EKWalberla::propagate() { + // first calculate the charge for the potential, for that get all the + // field-ids from the ekspecies pass the potential-field-id to the + // flux-kernels of the eks for this the integrate function has to be split + // with a public interface to diffusive and advective-flux this should also + // allow the back-coupling to the LB with a field-id + + if (ek_container->empty()) { + return; + } + + ek_container->reset_charge(); + for (auto const &ek_species : *ek_container) { + ek_container->add_charge(ek_species->get_density_id(), + ek_species->get_valency(), + ek_species->is_double_precision()); + } + ek_container->solve_poisson(); + + FieldsConnector connector{}; + System::get_system().lb.connect(connector); + for (auto const &ek_species : *ek_container) { + try { + ek_species->integrate(ek_container->get_potential_field_id(), + connector.velocity_field_id, + connector.force_field_id); + } catch (std::runtime_error const &e) { + runtimeErrorMsg() << e.what(); + } + } + + perform_reactions(); + + for (auto const &ek_species : *ek_container) { + ek_species->ghost_communication(); + } +} + +void EKWalberla::perform_reactions() { + for (auto const &ek_reaction : *ek_reactions) { + ek_reaction->perform_reaction(); + } +} + +void EKWalberla::veto_time_step(double time_step) const { + walberla_tau_sanity_checks("EK", ek_container->get_tau(), time_step); +} + +void EKWalberla::sanity_checks() const { + auto const &lattice = ek_container->get_lattice(); + auto const agrid = ::box_geo.length()[0] / lattice.get_grid_dimensions()[0]; + auto [my_left, my_right] = lattice.get_local_domain(); + my_left *= agrid; + my_right *= agrid; + walberla_agrid_sanity_checks("EK", my_left, my_right, agrid); + // EK time step and MD time step must agree + walberla_tau_sanity_checks("EK", ek_container->get_tau()); +} + +} // namespace EK + +#endif // WALBERLA diff --git a/src/core/ek/EKWalberla.hpp b/src/core/ek/EKWalberla.hpp new file mode 100644 index 00000000000..dba2d736e29 --- /dev/null +++ b/src/core/ek/EKWalberla.hpp @@ -0,0 +1,77 @@ +/* + * Copyright (C) 2023 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 . + */ + +#pragma once + +#include "config/config.hpp" + +#ifdef WALBERLA + +#include + +#include +#include +#include + +// forward declarations +class EKinWalberlaBase; +template class EKContainer; +template class EKReactions; +namespace walberla { +class EKReactionBase; +} + +namespace EK { + +struct EKWalberla { + using ek_container_type = EKContainer; + using ek_reactions_type = EKReactions; + std::shared_ptr ek_container; + std::shared_ptr ek_reactions; + + EKWalberla(std::shared_ptr ek_container_instance, + std::shared_ptr ek_reactions_instance) + : ek_container{std::move(ek_container_instance)}, + ek_reactions{std::move(ek_reactions_instance)} {} + + double get_tau() const; + void veto_time_step(double time_step) const; + void sanity_checks() const; + bool is_ready_for_propagation() const noexcept; + void propagate(); + void perform_reactions(); + + void on_cell_structure_change() const {} + void on_boxl_change() const { + throw std::runtime_error("MD cell geometry change not supported by EK"); + } + void on_node_grid_change() const { + throw std::runtime_error("MPI topology change not supported by EK"); + } + void on_timestep_change() const { + throw std::runtime_error("Time step change not supported by EK"); + } + void on_temperature_change() const { + throw std::runtime_error("Temperature change not supported by EK"); + } +}; + +} // namespace EK + +#endif // WALBERLA diff --git a/src/core/grid_based_algorithms/ek_container.hpp b/src/core/ek/Implementation.hpp similarity index 57% rename from src/core/grid_based_algorithms/ek_container.hpp rename to src/core/ek/Implementation.hpp index 1892bb5adda..a987b14a6af 100644 --- a/src/core/grid_based_algorithms/ek_container.hpp +++ b/src/core/ek/Implementation.hpp @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022-2023 The ESPResSo project + * Copyright (C) 2023 The ESPResSo project * * This file is part of ESPResSo. * @@ -17,32 +17,32 @@ * along with this program. If not, see . */ -#ifndef ESPRESSO_EK_CONTAINER_HPP -#define ESPRESSO_EK_CONTAINER_HPP +#pragma once #include "config/config.hpp" -#ifdef WALBERLA -#include -#include -#endif // WALBERLA +#include "ek/Solver.hpp" -#include +#include "ek/EKNone.hpp" +#include "ek/EKWalberla.hpp" -struct NoEKActive : public std::exception { - const char *what() const noexcept override { return "EK not activated"; } -}; +#include +#include +#include namespace EK { +using DiffusionAdvectionReactionActor = std::variant< #ifdef WALBERLA -extern EKContainer ek_container; -#endif // WALBERLA + std::shared_ptr, +#endif + std::shared_ptr>; -double get_tau(); -int get_steps_per_md_step(double md_timestep); -void propagate(); +struct Solver::Implementation { + /// @brief Main diffusion-advection-reaction solver. + std::optional solver; -} // namespace EK + Implementation() : solver{} {} +}; -#endif +} // namespace EK diff --git a/src/core/ek/Solver.cpp b/src/core/ek/Solver.cpp new file mode 100644 index 00000000000..54ab3175bb0 --- /dev/null +++ b/src/core/ek/Solver.cpp @@ -0,0 +1,139 @@ +/* + * Copyright (C) 2023 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 . + */ + +#include "config/config.hpp" + +#include "ek/Implementation.hpp" +#include "ek/Solver.hpp" +#include "ek/utils.hpp" + +#include "ek/EKNone.hpp" +#include "ek/EKWalberla.hpp" + +#include "integrate.hpp" +#include "system/System.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace EK { + +Solver::Solver() { impl = std::make_unique(); } + +static auto is_solver_set(std::unique_ptr const &ptr) { + return ptr != nullptr and ptr->solver.has_value(); +} + +static void check_solver(std::unique_ptr const &ptr) { + if (not is_solver_set(ptr)) { + throw NoEKActive{}; + } +} + +bool Solver::is_solver_set() const { return EK::is_solver_set(impl); } + +void Solver::reset() { System::get_system().ek.impl->solver = std::nullopt; } + +bool Solver::is_ready_for_propagation() const { + return is_solver_set() and + std::visit([](auto &ptr) { return ptr->is_ready_for_propagation(); }, + *impl->solver); +} + +void Solver::propagate() { + check_solver(impl); + std::visit([](auto &ptr) { ptr->propagate(); }, *impl->solver); +} + +void Solver::sanity_checks() const { + if (impl->solver) { + std::visit([](auto &ptr) { ptr->sanity_checks(); }, *impl->solver); + } +} + +void Solver::veto_time_step(double time_step) const { + if (impl->solver) { + std::visit([=](auto &ptr) { ptr->veto_time_step(time_step); }, + *impl->solver); + } +} + +void Solver::on_cell_structure_change() { + if (impl->solver) { + auto &solver = *impl->solver; + std::visit([](auto &ptr) { ptr->on_cell_structure_change(); }, solver); + } +} + +void Solver::on_boxl_change() { + if (impl->solver) { + std::visit([](auto &ptr) { ptr->on_boxl_change(); }, *impl->solver); + } +} + +void Solver::on_node_grid_change() { + if (impl->solver) { + std::visit([](auto &ptr) { ptr->on_node_grid_change(); }, *impl->solver); + } +} + +void Solver::on_timestep_change() { + if (impl->solver) { + std::visit([](auto &ptr) { ptr->on_timestep_change(); }, *impl->solver); + } +} + +void Solver::on_temperature_change() { + if (impl->solver) { + std::visit([](auto &ptr) { ptr->on_temperature_change(); }, *impl->solver); + } +} + +double Solver::get_tau() const { + check_solver(impl); + return std::visit([](auto &ptr) { return ptr->get_tau(); }, *impl->solver); +} + +template <> void Solver::set(std::shared_ptr ek_instance) { + assert(impl); + assert(not impl->solver.has_value()); + impl->solver = ek_instance; +} + +#ifdef WALBERLA +template <> +void Solver::set(std::shared_ptr ek_instance) { + assert(impl); + assert(not impl->solver.has_value()); + ek_instance->sanity_checks(); + impl->solver = ek_instance; +} +#endif // WALBERLA + +} // namespace EK diff --git a/src/core/ek/Solver.hpp b/src/core/ek/Solver.hpp new file mode 100644 index 00000000000..f11f835ce3c --- /dev/null +++ b/src/core/ek/Solver.hpp @@ -0,0 +1,101 @@ +/* + * Copyright (C) 2023 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 . + */ + +#pragma once + +#include "utils.hpp" + +#include +#include +#include +#include + +namespace EK { + +struct Solver { + struct Implementation; + + Solver(); + + /** @brief Return true if an @c EK solver is active. */ + [[nodiscard]] bool is_solver_set() const; + + /** @brief Return true if an @c EK solver can be propagated. */ + bool is_ready_for_propagation() const; + + /** @brief Remove the @c EK solver. */ + void reset(); + + /** + * @brief Set the @c EK solver. + * For developers: a specialization must exist for every new @c EK type. + */ + template void set(Args... args); + + /** + * @brief Connector to the implementation internal state. + * For developers: use this mechanism to access the underlying variant. + */ + template void connect(Connector &&connector) const { + assert(impl != nullptr); + connector(*this->impl); + } + + /** + * @brief Propagate the @c EK species. + */ + void propagate(); + + /** + * @brief Perform a full initialization of the lattice-Boltzmann system. + * All derived parameters and the fluid are reset to their default values. + */ + void init() const {} + + /** + * @brief Get the @c EK time step. + */ + double get_tau() const; + + auto get_steps_per_md_step(double time_step) const { + return static_cast(std::round(get_tau() / time_step)); + } + + /** + * @brief Perform @c EK parameter checks. + */ + void sanity_checks() const; + + /** + * @brief Check if a MD time step is compatible with the @c EK tau. + */ + void veto_time_step(double time_step) const; + + void on_boxl_change(); + void on_node_grid_change(); + void on_cell_structure_change(); + void on_timestep_change(); + void on_temperature_change(); + +private: + /** @brief Pointer-to-implementation. */ + std::unique_ptr impl; +}; + +} // namespace EK diff --git a/src/core/grid_based_algorithms/ek_reactions.cpp b/src/core/ek/ek_reactions.cpp similarity index 97% rename from src/core/grid_based_algorithms/ek_reactions.cpp rename to src/core/ek/ek_reactions.cpp index d65fc813984..c6689e3ed57 100644 --- a/src/core/grid_based_algorithms/ek_reactions.cpp +++ b/src/core/ek/ek_reactions.cpp @@ -21,7 +21,7 @@ #ifdef WALBERLA -#include "ek_reactions.hpp" +#include "ek/ek_reactions.hpp" #include diff --git a/src/core/grid_based_algorithms/ek_reactions.hpp b/src/core/ek/utils.hpp similarity index 66% rename from src/core/grid_based_algorithms/ek_reactions.hpp rename to src/core/ek/utils.hpp index 91958140dcc..8b2fc0097de 100644 --- a/src/core/grid_based_algorithms/ek_reactions.hpp +++ b/src/core/ek/utils.hpp @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 The ESPResSo project + * Copyright (C) 2023 The ESPResSo project * * This file is part of ESPResSo. * @@ -17,23 +17,14 @@ * along with this program. If not, see . */ -#ifndef ESPRESSO_EK_REACTIONS_HPP -#define ESPRESSO_EK_REACTIONS_HPP +#pragma once -#include "config/config.hpp" - -#ifdef WALBERLA - -#include "EKReactions.hpp" -#include "walberla_bridge/electrokinetics/reactions/EKReactionBase.hpp" +#include namespace EK { -extern EKReactions ek_reactions; - -void perform_reactions(); +struct NoEKActive : public std::exception { + const char *what() const noexcept override { return "EK not activated"; } +}; } // namespace EK - -#endif // WALBERLA -#endif diff --git a/src/core/event.cpp b/src/core/event.cpp index 9541c21235a..14b6957fda4 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -37,7 +37,6 @@ #include "electrostatics/icc.hpp" #include "errorhandling.hpp" #include "grid.hpp" -#include "grid_based_algorithms/lb_interface.hpp" #include "immersed_boundaries.hpp" #include "integrate.hpp" #include "interactions.hpp" @@ -82,12 +81,14 @@ void on_integration_start(double time_step) { /* sanity checks */ /********************************************/ + auto &system = System::get_system(); integrator_sanity_checks(); #ifdef NPT integrator_npt_sanity_checks(); #endif long_range_interactions_sanity_checks(); - LB::sanity_checks(time_step); + system.lb.sanity_checks(); + system.ek.sanity_checks(); /********************************************/ /* end sanity checks */ @@ -113,7 +114,7 @@ void on_integration_start(double time_step) { } #ifdef ELECTROSTATICS { - auto const &actor = System::get_system().coulomb.impl->solver; + auto const &actor = system.coulomb.impl->solver; if (not Utils::Mpi::all_compare(comm_cart, static_cast(actor)) or (actor and not Utils::Mpi::all_compare(comm_cart, (*actor).index()))) runtimeErrorMsg() << "Nodes disagree about Coulomb long-range method"; @@ -121,7 +122,7 @@ void on_integration_start(double time_step) { #endif #ifdef DIPOLES { - auto const &actor = System::get_system().dipoles.impl->solver; + auto const &actor = system.dipoles.impl->solver; if (not Utils::Mpi::all_compare(comm_cart, static_cast(actor)) or (actor and not Utils::Mpi::all_compare(comm_cart, (*actor).index()))) runtimeErrorMsg() << "Nodes disagree about dipolar long-range method"; @@ -228,49 +229,43 @@ void on_boxl_change(bool skip_method_adaption) { * therefore recalculate them. */ cells_re_init(cell_structure.decomposition_type()); + /* Now give methods a chance to react to the change in box length */ if (not skip_method_adaption) { - /* Now give methods a chance to react to the change in box length */ + auto &system = System::get_system(); + system.lb.on_boxl_change(); + system.ek.on_boxl_change(); #ifdef ELECTROSTATICS - System::get_system().coulomb.on_boxl_change(); + system.coulomb.on_boxl_change(); #endif - #ifdef DIPOLES - System::get_system().dipoles.on_boxl_change(); + system.dipoles.on_boxl_change(); #endif - - LB::init(); } } void on_cell_structure_change() { clear_particle_node(); - if (lattice_switch == ActiveLB::WALBERLA_LB) { - throw std::runtime_error( - "LB does not currently support handling changes of the MD cell " - "geometry. Setup the cell system, skin and interactions before " - "activating the CPU LB."); - } - /* Now give methods a chance to react to the change in cell structure. * Most ES methods need to reinitialize, as they depend on skin, * node grid and so on. */ -#if defined(ELECTROSTATICS) or defined(DIPOLES) if (System::is_system_set()) { + auto &system = System::get_system(); + system.lb.on_cell_structure_change(); + system.ek.on_cell_structure_change(); #ifdef ELECTROSTATICS - System::get_system().coulomb.on_cell_structure_change(); + system.coulomb.on_cell_structure_change(); #endif #ifdef DIPOLES - System::get_system().dipoles.on_cell_structure_change(); + system.dipoles.on_cell_structure_change(); #endif } -#endif } void on_temperature_change() { - if (lattice_switch != ActiveLB::NONE) { - throw std::runtime_error("Temperature change not supported by LB"); - } + auto &system = System::get_system(); + system.lb.on_temperature_change(); + system.ek.on_temperature_change(); } void on_periodicity_change() { @@ -300,9 +295,9 @@ void on_skin_change() { void on_thermostat_param_change() { reinit_thermo = true; } void on_timestep_change() { - if (lattice_switch != ActiveLB::NONE) { - throw std::runtime_error("Time step change not supported by LB"); - } + auto &system = System::get_system(); + system.lb.on_timestep_change(); + system.ek.on_timestep_change(); on_thermostat_param_change(); } @@ -310,11 +305,14 @@ void on_forcecap_change() { recalc_forces = true; } void on_node_grid_change() { grid_changed_n_nodes(); + auto &system = System::get_system(); + system.lb.on_node_grid_change(); + system.ek.on_node_grid_change(); #ifdef ELECTROSTATICS - System::get_system().coulomb.on_node_grid_change(); + system.coulomb.on_node_grid_change(); #endif #ifdef DIPOLES - System::get_system().dipoles.on_node_grid_change(); + system.dipoles.on_node_grid_change(); #endif cells_re_init(cell_structure.decomposition_type()); } @@ -328,7 +326,7 @@ unsigned global_ghost_flags() { /* Position and Properties are always requested. */ unsigned data_parts = Cells::DATA_PART_POSITION | Cells::DATA_PART_PROPERTIES; - if (lattice_switch == ActiveLB::WALBERLA_LB) + if (System::get_system().lb.is_solver_set()) data_parts |= Cells::DATA_PART_MOMENTUM; if (thermo_switch & THERMO_DPD) diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 351bc605ad5..b83b55b4dc0 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -35,11 +35,10 @@ #include "forcecap.hpp" #include "forces_inline.hpp" #include "galilei/ComFixed.hpp" -#include "grid_based_algorithms/lb_interface.hpp" -#include "grid_based_algorithms/lb_particle_coupling.hpp" #include "immersed_boundaries.hpp" #include "integrate.hpp" #include "interactions.hpp" +#include "lb/particle_coupling.hpp" #include "magnetostatics/dipoles.hpp" #include "nonbonded_interactions/VerletCriterion.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" @@ -239,7 +238,7 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) { // Must be done here. Forces need to be ghost-communicated immersed_boundaries.volume_conservation(cell_structure); - if (lattice_switch != ActiveLB::NONE) { + if (System::get_system().lb.is_solver_set()) { LB::couple_particles(thermo_virtual, particles, ghost_particles, time_step); } diff --git a/src/core/grid_based_algorithms/ek_container.cpp b/src/core/grid_based_algorithms/ek_container.cpp deleted file mode 100644 index 780a8eee574..00000000000 --- a/src/core/grid_based_algorithms/ek_container.cpp +++ /dev/null @@ -1,110 +0,0 @@ -/* - * Copyright (C) 2022-2023 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 . - */ - -#include "config/config.hpp" - -#include "ek_container.hpp" -#include "ek_reactions.hpp" -#include "errorhandling.hpp" -#include "lb_interface.hpp" -#include "lb_walberla_instance.hpp" - -#ifdef WALBERLA -#include -#endif // WALBERLA - -#include - -#ifdef WALBERLA -#include -#include -#include -#endif // WALBERLA - -namespace EK { - -#ifdef WALBERLA -EKContainer ek_container; -#endif // WALBERLA - -double get_tau() { -#ifdef WALBERLA - return ek_container.get_tau(); -#else - throw NoEKActive(); -#endif // WALBERLA -} - -int get_steps_per_md_step(double md_timestep) { - return static_cast(std::round(get_tau() / md_timestep)); -} - -void propagate() { -#ifdef WALBERLA - // first calculate the charge for the potential, for that get all the - // field-ids from the ekspecies pass the potential-field-id to the - // flux-kernels of the eks for this the integrate function has to be split - // with a public interface to diffusive and advective-flux this should also - // allow the back-coupling to the LB with a field-id - - if (ek_container.empty()) { - return; - } - - if (!ek_container.is_poisson_solver_set()) { - runtimeErrorMsg() << "EK requires a Poisson solver."; - return; - } - - ek_container.reset_charge(); - std::for_each(ek_container.begin(), ek_container.end(), [](auto const &ek) { - ek_container.add_charge(ek->get_density_id(), ek->get_valency(), - ek->is_double_precision()); - }); - - ek_container.solve_poisson(); - - auto velocity_field_id = std::size_t{}; - auto force_field_id = std::size_t{}; - try { - auto const lbf = ::lb_walberla(); - velocity_field_id = lbf->get_velocity_field_id(); - force_field_id = lbf->get_force_field_id(); - } catch (std::runtime_error const &) { - } - - std::for_each(ek_container.begin(), ek_container.end(), - [velocity_field_id, force_field_id](auto const &ek) { - try { - ek->integrate(ek_container.get_potential_field_id(), - velocity_field_id, force_field_id); - } catch (std::runtime_error const &e) { - runtimeErrorMsg() << e.what(); - } - }); - - EK::perform_reactions(); - - for (auto const &species : ek_container) { - species->ghost_communication(); - } -#endif // WALBERLA -} - -} // namespace EK diff --git a/src/core/grid_based_algorithms/lb_interface.cpp b/src/core/grid_based_algorithms/lb_interface.cpp deleted file mode 100644 index 076ede64683..00000000000 --- a/src/core/grid_based_algorithms/lb_interface.cpp +++ /dev/null @@ -1,170 +0,0 @@ -/* - * Copyright (C) 2010-2022 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 . - */ - -#include "grid_based_algorithms/lb_interface.hpp" -#include "grid_based_algorithms/lb_walberla_instance.hpp" - -#include "BoxGeometry.hpp" -#include "config/config.hpp" -#include "errorhandling.hpp" -#include "grid.hpp" - -#include - -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -ActiveLB lattice_switch = ActiveLB::NONE; - -namespace LB { - -ActiveLB get_lattice_switch() { return lattice_switch; } - -int get_steps_per_md_step(double md_timestep) { - return static_cast(std::round(get_tau() / md_timestep)); -} - -void init() {} - -void propagate() { - if (lattice_switch == ActiveLB::WALBERLA_LB) { -#ifdef WALBERLA - lb_walberla()->integrate(); -#endif - } -} - -void sanity_checks(double time_step) { - if (lattice_switch == ActiveLB::NONE) - return; - - if (lattice_switch == ActiveLB::WALBERLA_LB) { -#ifdef WALBERLA - lb_sanity_checks(*lb_walberla(), *lb_walberla_params(), time_step); -#endif - } -} - -void lebc_sanity_checks(unsigned int shear_direction, - unsigned int shear_plane_normal) { - if (lattice_switch == ActiveLB::WALBERLA_LB) { -#ifdef WALBERLA - lb_walberla()->check_lebc(shear_direction, shear_plane_normal); -#endif - } -} - -double get_agrid() { - if (lattice_switch == ActiveLB::WALBERLA_LB) { -#ifdef WALBERLA - return lb_walberla_params()->get_agrid(); -#endif - } - throw NoLBActive(); -} - -void check_tau_time_step_consistency(double tau, double time_step) { - // use float epsilon since tau may be a float - auto const eps = static_cast(std::numeric_limits::epsilon()); - if ((tau - time_step) / (tau + time_step) < -eps) - throw std::invalid_argument("LB tau (" + std::to_string(tau) + - ") must be >= MD time_step (" + - std::to_string(time_step) + ")"); - auto const factor = tau / time_step; - if (fabs(round(factor) - factor) / factor > eps) - throw std::invalid_argument("LB tau (" + std::to_string(tau) + - ") must be an integer multiple of the " - "MD time_step (" + - std::to_string(time_step) + "). Factor is " + - std::to_string(factor)); -} - -double get_tau() { -#ifdef WALBERLA - if (lattice_switch == ActiveLB::WALBERLA_LB) { - return lb_walberla_params()->get_tau(); - } -#endif - throw NoLBActive(); -} - -double get_kT() { - if (lattice_switch == ActiveLB::WALBERLA_LB) { -#ifdef WALBERLA - return lb_walberla()->get_kT(); -#endif - } - throw NoLBActive(); -} - -double get_lattice_speed() { return get_agrid() / get_tau(); } - -Utils::VectorXd<9> const -get_pressure_tensor(boost::mpi::communicator const &comm) { - if (lattice_switch == ActiveLB::WALBERLA_LB) { -#ifdef WALBERLA - auto const local_pressure_tensor = lb_walberla()->get_pressure_tensor(); - std::remove_const_t pressure_tensor; - boost::mpi::reduce(comm, local_pressure_tensor, pressure_tensor, - std::plus<>(), 0); - return pressure_tensor; -#endif - } - throw NoLBActive(); -} - -Utils::Vector3d calc_fluid_momentum() { - if (lattice_switch == ActiveLB::WALBERLA_LB) { -#ifdef WALBERLA - return lb_walberla()->get_momentum(); -#endif - } - throw NoLBActive(); -} - -std::optional -get_interpolated_velocity(Utils::Vector3d const &pos) { - if (lattice_switch == ActiveLB::WALBERLA_LB) { -#ifdef WALBERLA - auto const folded_pos = folded_position(pos, box_geo) / get_agrid(); - return lb_walberla()->get_velocity_at_pos(folded_pos); -#endif - } - throw NoLBActive(); -} - -std::optional get_interpolated_density(Utils::Vector3d const &pos) { - if (lattice_switch == ActiveLB::WALBERLA_LB) { -#ifdef WALBERLA - auto const folded_pos = folded_position(pos, box_geo) / get_agrid(); - return lb_walberla()->get_density_at_pos(folded_pos); -#endif - } - throw NoLBActive(); -} - -} // namespace LB diff --git a/src/core/grid_based_algorithms/lb_interface.hpp b/src/core/grid_based_algorithms/lb_interface.hpp deleted file mode 100644 index 892ded2ec37..00000000000 --- a/src/core/grid_based_algorithms/lb_interface.hpp +++ /dev/null @@ -1,131 +0,0 @@ -/* - * Copyright (C) 2010-2022 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 . - */ -#ifndef CORE_LB_INTERFACE -#define CORE_LB_INTERFACE - -#include - -#include -#include - -// Forward Declarations -namespace boost::mpi { -class communicator; -} // namespace boost::mpi - -/** @brief LB implementation currently active. */ -enum class ActiveLB : int { NONE, WALBERLA_LB }; - -/** @brief Switch determining the type of lattice dynamics. */ -extern ActiveLB lattice_switch; - -struct NoLBActive : public std::exception { - const char *what() const noexcept override { return "LB not activated"; } -}; - -namespace LB { - -/** - * @brief Get the global variable @ref lattice_switch. - */ -ActiveLB get_lattice_switch(); - -int get_steps_per_md_step(double md_timestep); - -/** - * @brief Propagate the LB fluid. - */ -void propagate(); - -/** - * @brief Perform a full initialization of the lattice-Boltzmann system. - * All derived parameters and the fluid are reset to their default values. - */ -void init(); - -/** - * @brief Check if tau is an integer multiple of time_step, throws if not - */ -void check_tau_time_step_consistency(double tau, double time_step); - -/** - * @brief Perform LB parameter and boundary velocity checks. - */ -void sanity_checks(double time_step); - -/** - * @brief Perform LB LEbc parameter checks. - */ -void lebc_sanity_checks(unsigned int shear_direction, - unsigned int shear_plane_normal); - -/** - * @brief Set the LB fluid velocity for a single node. - */ -void set_velocity(Utils::Vector3i const &ind, Utils::Vector3d const &u); - -/** - * @brief Get the LB time step. - */ -double get_tau(); - -/** - * @brief Get the LB grid spacing. - */ -double get_agrid(); - -/** - * @brief Get the thermal energy parameter of the LB fluid. - */ -double get_kT(); - -/** - * @brief Get the lattice speed (agrid/tau). - */ -double get_lattice_speed(); - -/** @brief Calculate the average pressure tensor of all nodes by accumulating - * over all nodes and dividing by the number of nodes. - * Returns the lower triangle of the LB pressure tensor. - */ -Utils::VectorXd<9> const -get_pressure_tensor(boost::mpi::communicator const &comm); - -Utils::Vector3d calc_fluid_momentum(); - -/** - * @brief Calculates the interpolated fluid velocity on the responsible node - * process. - * @param pos Position at which the velocity is to be calculated. - * @retval optional interpolated fluid velocity. - */ -std::optional -get_interpolated_velocity(Utils::Vector3d const &pos); - -/** - * @brief Calculates the interpolated fluid density on the responsible node - * process. - * @param pos Position at which the density is to be calculated. - * @retval optional interpolated fluid density. - */ -std::optional get_interpolated_density(Utils::Vector3d const &pos); - -} // namespace LB - -#endif diff --git a/src/core/grid_based_algorithms/lb_interpolation.cpp b/src/core/grid_based_algorithms/lb_interpolation.cpp deleted file mode 100644 index 05cacf3b2f4..00000000000 --- a/src/core/grid_based_algorithms/lb_interpolation.cpp +++ /dev/null @@ -1,60 +0,0 @@ -/* - * Copyright (C) 2010-2022 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 . - */ - -#include "grid_based_algorithms/lb_interpolation.hpp" -#include "grid_based_algorithms/lb_interface.hpp" -#include "grid_based_algorithms/lb_walberla_instance.hpp" - -#include "communication.hpp" -#include "config/config.hpp" - -#include - -#include -#include - -const Utils::Vector3d -lb_lbinterpolation_get_interpolated_velocity(const Utils::Vector3d &pos) { - /* calculate fluid velocity at particle's position - this is done by linear interpolation - (Eq. (11) Ahlrichs and Duenweg, JCP 111(17):8225 (1999)) */ - if (lattice_switch == ActiveLB::WALBERLA_LB) { -#ifdef WALBERLA - auto res = lb_walberla()->get_velocity_at_pos(pos / LB::get_agrid(), true); - if (!res) { - std::cout << this_node << ": position: [" << pos << "]\n"; - throw std::runtime_error( - "Interpolated velocity could not be obtained from Walberla"); - } - return *res; -#endif - } - throw std::runtime_error("No LB active."); -} - -void lb_lbinterpolation_add_force_density( - const Utils::Vector3d &pos, const Utils::Vector3d &force_density) { - if (lattice_switch == ActiveLB::WALBERLA_LB) { -#ifdef WALBERLA - if (!lb_walberla()->add_force_at_pos(pos / LB::get_agrid(), force_density)) - throw std::runtime_error("Could not apply force to lb."); -#endif - } else - throw std::runtime_error("No LB active."); -} diff --git a/src/core/grid_based_algorithms/lb_walberla_instance.cpp b/src/core/grid_based_algorithms/lb_walberla_instance.cpp deleted file mode 100644 index 18998f14b52..00000000000 --- a/src/core/grid_based_algorithms/lb_walberla_instance.cpp +++ /dev/null @@ -1,104 +0,0 @@ -/* - * Copyright (C) 2019-2023 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 . - */ -#include "config/config.hpp" - -#ifdef WALBERLA -#include "lb_walberla_instance.hpp" - -#include "communication.hpp" -#include "errorhandling.hpp" -#include "grid.hpp" -#include "integrate.hpp" -#include "lb_interface.hpp" - -#include - -#include - -#include -#include - -static std::weak_ptr lb_walberla_instance; -static std::shared_ptr lb_walberla_params_instance; - -std::shared_ptr lb_walberla() { - auto lb_walberla_instance_handle = ::lb_walberla_instance.lock(); - if (not lb_walberla_instance_handle) { - throw std::runtime_error( - "Attempted access to uninitialized LBWalberla instance."); - } - return lb_walberla_instance_handle; -} - -std::shared_ptr lb_walberla_params() { - if (not ::lb_walberla_params_instance) { - throw std::runtime_error( - "Attempted access to uninitialized LBWalberlaParams instance."); - } - return ::lb_walberla_params_instance; -} - -void lb_sanity_checks(LBWalberlaBase const &lb_fluid, - LBWalberlaParams const &lb_params, double md_time_step) { - auto const agrid = lb_params.get_agrid(); - auto const tau = lb_params.get_tau(); - // waLBerla and ESPResSo must agree on domain decomposition - auto [lb_my_left, lb_my_right] = lb_fluid.get_lattice().get_local_domain(); - lb_my_left *= agrid; - lb_my_right *= agrid; - auto const my_left = local_geo.my_left(); - auto const my_right = local_geo.my_right(); - auto const tol = agrid / 1E6; - if ((lb_my_left - my_left).norm2() > tol or - (lb_my_right - my_right).norm2() > tol) { - runtimeErrorMsg() << "\nMPI rank " << this_node << ": " - << "left ESPResSo: [" << my_left << "], " - << "left waLBerla: [" << lb_my_left << "]" - << "\nMPI rank " << this_node << ": " - << "right ESPResSo: [" << my_right << "], " - << "right waLBerla: [" << lb_my_right << "]"; - throw std::runtime_error( - "waLBerla and ESPResSo disagree about domain decomposition."); - } - // LB time step and MD time step must agree - if (md_time_step > 0.) { - LB::check_tau_time_step_consistency(tau, md_time_step); - } -} - -void activate_lb_walberla(std::shared_ptr lb_fluid, - std::shared_ptr lb_params) { - if (::lattice_switch != ActiveLB::NONE) { - throw std::runtime_error("Cannot add a second LB instance"); - } - lb_sanity_checks(*lb_fluid, *lb_params, get_time_step()); - auto const &lebc = ::box_geo.lees_edwards_bc(); - lb_fluid->check_lebc(lebc.shear_direction, lebc.shear_plane_normal); - ::lb_walberla_instance = std::weak_ptr{lb_fluid}; - ::lb_walberla_params_instance = lb_params; - ::lattice_switch = ActiveLB::WALBERLA_LB; -} - -void deactivate_lb_walberla() { - ::lb_walberla_instance.reset(); - ::lb_walberla_params_instance.reset(); - ::lattice_switch = ActiveLB::NONE; -} - -#endif diff --git a/src/core/grid_based_algorithms/lb_walberla_instance.hpp b/src/core/grid_based_algorithms/lb_walberla_instance.hpp deleted file mode 100644 index 90bf1c65a8e..00000000000 --- a/src/core/grid_based_algorithms/lb_walberla_instance.hpp +++ /dev/null @@ -1,58 +0,0 @@ -/* - * Copyright (C) 2019-2023 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 . - */ -#ifndef GRID_BASED_ALGORITHMS_LBWALBERLA_INSTANCE_HPP -#define GRID_BASED_ALGORITHMS_LBWALBERLA_INSTANCE_HPP - -#include "config/config.hpp" - -#ifdef WALBERLA - -#include - -#include - -struct LBWalberlaParams { - LBWalberlaParams(double agrid, double tau) : m_agrid(agrid), m_tau(tau) {} - double get_agrid() const { return m_agrid; }; - double get_tau() const { return m_tau; }; - -private: - double m_agrid; - double m_tau; -}; - -/** @brief Access the per-MPI-node waLBerla LB instance */ -std::shared_ptr lb_walberla(); - -/** @brief Access the waLBerla parameters */ -std::shared_ptr lb_walberla_params(); - -void lb_sanity_checks(LBWalberlaBase const &lb_fluid, - LBWalberlaParams const &lb_params, double md_time_step); - -/** @brief Register a waLBerla LB instance and update lattice switch. */ -void activate_lb_walberla(std::shared_ptr lb_fluid, - std::shared_ptr lb_params); - -/** @brief De-register a waLBerla LB instance and update lattice switch. */ -void deactivate_lb_walberla(); - -#endif // WALBERLA - -#endif diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 23b7f323263..2b28a058d3b 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -44,16 +44,16 @@ #include "event.hpp" #include "forces.hpp" #include "grid.hpp" -#include "grid_based_algorithms/ek_container.hpp" -#include "grid_based_algorithms/lb_interface.hpp" -#include "grid_based_algorithms/lb_particle_coupling.hpp" #include "interactions.hpp" +#include "lb/particle_coupling.hpp" +#include "lb/utils.hpp" #include "lees_edwards/lees_edwards.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "npt.hpp" #include "rattle.hpp" #include "rotation.hpp" #include "signalhandling.hpp" +#include "system/System.hpp" #include "thermostat.hpp" #include "virtual_sites.hpp" @@ -74,7 +74,9 @@ #include #include #include +#include #include +#include #include #ifdef WALBERLA @@ -189,6 +191,53 @@ void integrator_sanity_checks() { } } +#ifdef WALBERLA +void walberla_tau_sanity_checks(std::string method, double tau, + double time_step) { + if (time_step <= 0.) { + return; + } + // use float epsilon since tau may be a float + auto const eps = static_cast(std::numeric_limits::epsilon()); + if ((tau - time_step) / (tau + time_step) < -eps) + throw std::invalid_argument(method + " tau (" + std::to_string(tau) + + ") must be >= MD time_step (" + + std::to_string(time_step) + ")"); + auto const factor = tau / time_step; + if (std::fabs(std::round(factor) - factor) / factor > eps) + throw std::invalid_argument(method + " tau (" + std::to_string(tau) + + ") must be an integer multiple of the " + "MD time_step (" + + std::to_string(time_step) + "). Factor is " + + std::to_string(factor)); +} + +void walberla_tau_sanity_checks(std::string method, double tau) { + walberla_tau_sanity_checks(std::move(method), tau, ::time_step); +} + +void walberla_agrid_sanity_checks(std::string method, + Utils::Vector3d const &lattice_left, + Utils::Vector3d const &lattice_right, + double agrid) { + // waLBerla and ESPResSo must agree on domain decomposition + auto const &my_left = ::local_geo.my_left(); + auto const &my_right = ::local_geo.my_right(); + auto const tol = agrid / 1E6; + if ((lattice_left - my_left).norm2() > tol or + (lattice_right - my_right).norm2() > tol) { + runtimeErrorMsg() << "\nMPI rank " << ::this_node << ": " + << "left ESPResSo: [" << my_left << "], " + << "left waLBerla: [" << lattice_left << "]" + << "\nMPI rank " << ::this_node << ": " + << "right ESPResSo: [" << my_right << "], " + << "right waLBerla: [" << lattice_right << "]"; + throw std::runtime_error( + "waLBerla and ESPResSo disagree about domain decomposition."); + } +} +#endif + static void resort_particles_if_needed(ParticleRange const &particles) { auto const offset = LeesEdwards::verlet_list_offset( box_geo, cell_structure.get_le_pos_offset_at_last_resort()); @@ -263,6 +312,8 @@ int integrate(int n_steps, int reuse_forces) { CALI_CXX_MARK_FUNCTION; #endif + auto &system = System::get_system(); + // Prepare particle structure and run sanity checks of all active algorithms on_integration_start(time_step); @@ -312,6 +363,13 @@ int integrate(int n_steps, int reuse_forces) { auto caught_sigint = false; auto caught_error = false; + auto lb_active = false; + auto ek_active = false; + if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) { + lb_active = system.lb.is_solver_set(); + ek_active = system.ek.is_ready_for_propagation(); + } + #ifdef VALGRIND CALLGRIND_START_INSTRUMENTATION; #endif @@ -384,17 +442,12 @@ int integrate(int n_steps, int reuse_forces) { // propagate one-step functionalities if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) { - auto const lb_active = LB::get_lattice_switch() != ActiveLB::NONE; -#ifdef WALBERLA - auto const ek_active = not EK::ek_container.empty(); -#else - auto constexpr ek_active = false; -#endif - if (lb_active and ek_active) { + auto &lb = system.lb; + auto &ek = system.ek; // assume that they are coupled, which is not necessarily true - auto const lb_steps_per_md_step = LB::get_steps_per_md_step(time_step); - auto const ek_steps_per_md_step = EK::get_steps_per_md_step(time_step); + auto const lb_steps_per_md_step = lb.get_steps_per_md_step(time_step); + auto const ek_steps_per_md_step = ek.get_steps_per_md_step(time_step); if (lb_steps_per_md_step != ek_steps_per_md_step) { runtimeErrorMsg() @@ -407,24 +460,26 @@ int integrate(int n_steps, int reuse_forces) { fluid_step += 1; if (fluid_step >= lb_steps_per_md_step) { fluid_step = 0; - LB::propagate(); - EK::propagate(); + lb.propagate(); + ek.propagate(); } lb_lbcoupling_propagate(); } else if (lb_active) { - auto const lb_steps_per_md_step = LB::get_steps_per_md_step(time_step); + auto &lb = system.lb; + auto const lb_steps_per_md_step = lb.get_steps_per_md_step(time_step); fluid_step += 1; if (fluid_step >= lb_steps_per_md_step) { fluid_step = 0; - LB::propagate(); + lb.propagate(); } lb_lbcoupling_propagate(); } else if (ek_active) { - auto const ek_steps_per_md_step = EK::get_steps_per_md_step(time_step); + auto &ek = system.ek; + auto const ek_steps_per_md_step = ek.get_steps_per_md_step(time_step); ek_step += 1; if (ek_step >= ek_steps_per_md_step) { ek_step = 0; - EK::propagate(); + ek.propagate(); } } @@ -562,8 +617,12 @@ void increment_sim_time(double amount) { sim_time += amount; } void set_time_step(double value) { if (value <= 0.) throw std::domain_error("time_step must be > 0."); - if (LB::get_lattice_switch() != ActiveLB::NONE) { - LB::check_tau_time_step_consistency(LB::get_tau(), value); + auto &system = System::get_system(); + if (system.lb.is_solver_set()) { + system.lb.veto_time_step(value); + } + if (system.ek.is_solver_set()) { + system.ek.veto_time_step(value); } ::time_step = value; on_timestep_change(); diff --git a/src/core/integrate.hpp b/src/core/integrate.hpp index 50658357631..9ba99e40dcb 100644 --- a/src/core/integrate.hpp +++ b/src/core/integrate.hpp @@ -27,6 +27,14 @@ * Implementation in \ref integrate.cpp. */ +#include "config/config.hpp" + +#ifdef WALBERLA +#include + +#include +#endif + /** \name Integrator switches */ /**@{*/ #define INTEG_METHOD_NPT_ISO 0 @@ -68,6 +76,16 @@ double interaction_range(); */ void integrator_sanity_checks(); +#ifdef WALBERLA +void walberla_tau_sanity_checks(std::string method, double tau, + double time_step); +void walberla_tau_sanity_checks(std::string method, double tau); +void walberla_agrid_sanity_checks(std::string method, + Utils::Vector3d const &lattice_left, + Utils::Vector3d const &lattice_right, + double agrid); +#endif // WALBERLA + /** Integrate equations of motion * @param n_steps Number of integration steps, can be zero * @param reuse_forces Decide when to re-calculate forces diff --git a/src/core/lb/CMakeLists.txt b/src/core/lb/CMakeLists.txt new file mode 100644 index 00000000000..a876f219f92 --- /dev/null +++ b/src/core/lb/CMakeLists.txt @@ -0,0 +1,24 @@ +# +# Copyright (C) 2023 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 . +# + +target_sources(espresso_core PRIVATE Solver.cpp particle_coupling.cpp) + +if(ESPRESSO_BUILD_WITH_WALBERLA) + target_sources(espresso_core PRIVATE LBWalberla.cpp) +endif() diff --git a/src/core/grid_based_algorithms/lb_interpolation.hpp b/src/core/lb/Implementation.hpp similarity index 51% rename from src/core/grid_based_algorithms/lb_interpolation.hpp rename to src/core/lb/Implementation.hpp index 26de6289db6..35b50dafc15 100644 --- a/src/core/grid_based_algorithms/lb_interpolation.hpp +++ b/src/core/lb/Implementation.hpp @@ -1,5 +1,5 @@ /* - * Copyright (C) 2010-2022 The ESPResSo project + * Copyright (C) 2023 The ESPResSo project * * This file is part of ESPResSo. * @@ -16,23 +16,33 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef LATTICE_INTERPOLATION_HPP -#define LATTICE_INTERPOLATION_HPP -#include +#pragma once -/** - * @brief Calculates the fluid velocity at a given position of the - * lattice. - * @note It can lead to undefined behaviour if the - * position is not within the local lattice. - */ -const Utils::Vector3d -lb_lbinterpolation_get_interpolated_velocity(const Utils::Vector3d &pos); +#include "config/config.hpp" -/** - * @brief Add a force density to the fluid at the given position. - */ -void lb_lbinterpolation_add_force_density(const Utils::Vector3d &pos, - const Utils::Vector3d &force_density); +#include "lb/Solver.hpp" + +#include "lb/LBNone.hpp" +#include "lb/LBWalberla.hpp" + +#include +#include +#include + +namespace LB { + +using HydrodynamicsActor = std::variant< +#ifdef WALBERLA + std::shared_ptr, #endif + std::shared_ptr>; + +struct Solver::Implementation { + /// @brief Main hydrodynamics solver. + std::optional solver; + + Implementation() : solver{} {} +}; + +} // namespace LB diff --git a/src/core/lb/LBNone.hpp b/src/core/lb/LBNone.hpp new file mode 100644 index 00000000000..4920f1e4d08 --- /dev/null +++ b/src/core/lb/LBNone.hpp @@ -0,0 +1,57 @@ +/* + * Copyright (C) 2023 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 . + */ + +#pragma once + +#include "utils.hpp" + +namespace LB { + +struct LBNone { + void propagate() { throw NoLBActive{}; } + double get_agrid() const { throw NoLBActive{}; } + double get_tau() const { throw NoLBActive{}; } + double get_kT() const { throw NoLBActive{}; } + Utils::VectorXd<9> get_pressure_tensor() const { throw NoLBActive{}; } + std::optional get_velocity_at_pos(Utils::Vector3d const &, + bool) const { + throw NoLBActive{}; + } + std::optional get_density_at_pos(Utils::Vector3d const &, + bool) const { + throw NoLBActive{}; + } + bool add_force_at_pos(Utils::Vector3d const &pos, + Utils::Vector3d const &force) const { + throw NoLBActive{}; + } + Utils::Vector3d get_momentum() const { throw NoLBActive{}; } + void veto_time_step(double) const { throw NoLBActive{}; } + void sanity_checks() const { throw NoLBActive{}; } + void lebc_sanity_checks(unsigned int, unsigned int) const { + throw NoLBActive{}; + } + void on_cell_structure_change() const { throw NoLBActive{}; } + void on_boxl_change() const { throw NoLBActive{}; } + void on_node_grid_change() const { throw NoLBActive{}; } + void on_timestep_change() const { throw NoLBActive{}; } + void on_temperature_change() const { throw NoLBActive{}; } +}; + +} // namespace LB diff --git a/src/core/lb/LBWalberla.cpp b/src/core/lb/LBWalberla.cpp new file mode 100644 index 00000000000..9aad8f4635e --- /dev/null +++ b/src/core/lb/LBWalberla.cpp @@ -0,0 +1,89 @@ +/* + * Copyright (C) 2019-2023 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 . + */ +#include "config/config.hpp" + +#ifdef WALBERLA + +#include "LBWalberla.hpp" + +#include "communication.hpp" +#include "errorhandling.hpp" +#include "grid.hpp" +#include "integrate.hpp" +#include "system/System.hpp" + +#include + +#include + +#include + +namespace LB { + +double LBWalberla::get_kT() const { return lb_fluid->get_kT(); } + +Utils::VectorXd<9> LBWalberla::get_pressure_tensor() const { + return lb_fluid->get_pressure_tensor(); +} + +void LBWalberla::propagate() { lb_fluid->integrate(); } + +void LBWalberla::lebc_sanity_checks(unsigned int shear_direction, + unsigned int shear_plane_normal) const { + lb_fluid->check_lebc(shear_direction, shear_plane_normal); +} + +std::optional +LBWalberla::get_velocity_at_pos(Utils::Vector3d const &pos, + bool consider_points_in_halo) const { + return lb_fluid->get_velocity_at_pos(pos, consider_points_in_halo); +} + +std::optional +LBWalberla::get_density_at_pos(Utils::Vector3d const &pos, + bool consider_points_in_halo) const { + return lb_fluid->get_density_at_pos(pos, consider_points_in_halo); +} + +Utils::Vector3d LBWalberla::get_momentum() const { + return lb_fluid->get_momentum(); +} + +bool LBWalberla::add_force_at_pos(Utils::Vector3d const &pos, + Utils::Vector3d const &force) { + return lb_fluid->add_force_at_pos(pos, force); +} + +void LBWalberla::veto_time_step(double time_step) const { + walberla_tau_sanity_checks("LB", lb_params->get_tau(), time_step); +} + +void LBWalberla::sanity_checks() const { + auto const agrid = lb_params->get_agrid(); + auto [my_left, my_right] = lb_fluid->get_lattice().get_local_domain(); + my_left *= agrid; + my_right *= agrid; + walberla_agrid_sanity_checks("LB", my_left, my_right, agrid); + // LB time step and MD time step must agree + walberla_tau_sanity_checks("LB", lb_params->get_tau()); +} + +} // namespace LB + +#endif // WALBERLA diff --git a/src/core/lb/LBWalberla.hpp b/src/core/lb/LBWalberla.hpp new file mode 100644 index 00000000000..56919928e97 --- /dev/null +++ b/src/core/lb/LBWalberla.hpp @@ -0,0 +1,89 @@ +/* + * Copyright (C) 2023 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 . + */ + +#pragma once + +#include "config/config.hpp" + +#ifdef WALBERLA + +#include + +#include +#include +#include +#include + +class LBWalberlaBase; + +namespace LB { + +struct LBWalberlaParams { + LBWalberlaParams(double agrid, double tau) : m_agrid(agrid), m_tau(tau) {} + double get_agrid() const { return m_agrid; }; + double get_tau() const { return m_tau; }; + +private: + double m_agrid; + double m_tau; +}; + +struct LBWalberla { + std::shared_ptr lb_fluid; + std::shared_ptr lb_params; + LBWalberla(std::shared_ptr lb_fluid_, + std::shared_ptr lb_params_) + : lb_fluid{std::move(lb_fluid_)}, lb_params{std::move(lb_params_)} {} + double get_kT() const; + auto get_tau() const { return lb_params->get_tau(); } + auto get_agrid() const { return lb_params->get_agrid(); } + auto get_lattice_speed() const { return get_agrid() / get_tau(); } + Utils::VectorXd<9> get_pressure_tensor() const; + std::optional + get_velocity_at_pos(Utils::Vector3d const &pos, + bool consider_points_in_halo) const; + std::optional get_density_at_pos(Utils::Vector3d const &pos, + bool consider_points_in_halo) const; + Utils::Vector3d get_momentum() const; + bool add_force_at_pos(Utils::Vector3d const &pos, + Utils::Vector3d const &force); + void propagate(); + void veto_time_step(double time_step) const; + void sanity_checks() const; + void lebc_sanity_checks(unsigned int shear_direction, + unsigned int shear_plane_normal) const; + + void on_cell_structure_change() const {} + void on_boxl_change() const { + throw std::runtime_error("MD cell geometry change not supported by LB"); + } + void on_node_grid_change() const { + throw std::runtime_error("MPI topology change not supported by LB"); + } + void on_timestep_change() const { + throw std::runtime_error("Time step change not supported by LB"); + } + void on_temperature_change() const { + throw std::runtime_error("Temperature change not supported by LB"); + } +}; + +} // namespace LB + +#endif // WALBERLA diff --git a/src/core/lb/Solver.cpp b/src/core/lb/Solver.cpp new file mode 100644 index 00000000000..5e7e7c9fce6 --- /dev/null +++ b/src/core/lb/Solver.cpp @@ -0,0 +1,222 @@ +/* + * Copyright (C) 2023 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 . + */ + +#include "config/config.hpp" + +#include "lb/Implementation.hpp" +#include "lb/Solver.hpp" +#include "lb/utils.hpp" + +#include "lb/LBNone.hpp" +#include "lb/LBWalberla.hpp" + +#include "BoxGeometry.hpp" +#include "grid.hpp" +#include "system/System.hpp" + +#ifdef WALBERLA +#include +#endif + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace LB { + +Solver::Solver() { impl = std::make_unique(); } + +static auto is_solver_set(std::unique_ptr const &ptr) { + return ptr != nullptr and ptr->solver.has_value(); +} + +static void check_solver(std::unique_ptr const &ptr) { + if (not is_solver_set(ptr)) { + throw NoLBActive{}; + } +} + +bool Solver::is_solver_set() const { return LB::is_solver_set(impl); } + +void Solver::reset() { System::get_system().lb.impl->solver = std::nullopt; } + +void Solver::propagate() { + check_solver(impl); + std::visit([](auto &ptr) { ptr->propagate(); }, *impl->solver); +} + +void Solver::sanity_checks() const { + if (impl->solver) { + std::visit([](auto &ptr) { ptr->sanity_checks(); }, *impl->solver); + } +} + +void Solver::veto_time_step(double time_step) const { + if (impl->solver) { + std::visit([=](auto &ptr) { ptr->veto_time_step(time_step); }, + *impl->solver); + } +} + +void Solver::lebc_sanity_checks(unsigned int shear_direction, + unsigned int shear_plane_normal) const { + if (impl->solver) { + auto const callback = [=](auto &ptr) { + ptr->lebc_sanity_checks(shear_direction, shear_plane_normal); + }; + std::visit(callback, *impl->solver); + } +} + +void Solver::on_cell_structure_change() { + if (impl->solver) { + auto &solver = *impl->solver; + std::visit([](auto &ptr) { ptr->on_cell_structure_change(); }, solver); + } +} + +void Solver::on_boxl_change() { + if (impl->solver) { + std::visit([](auto &ptr) { ptr->on_boxl_change(); }, *impl->solver); + } +} + +void Solver::on_node_grid_change() { + if (impl->solver) { + std::visit([](auto &ptr) { ptr->on_node_grid_change(); }, *impl->solver); + } +} + +void Solver::on_timestep_change() { + if (impl->solver) { + std::visit([](auto &ptr) { ptr->on_timestep_change(); }, *impl->solver); + } +} + +void Solver::on_temperature_change() { + if (impl->solver) { + std::visit([](auto &ptr) { ptr->on_temperature_change(); }, *impl->solver); + } +} + +double Solver::get_agrid() const { + check_solver(impl); + return std::visit([](auto &ptr) { return ptr->get_agrid(); }, *impl->solver); +} + +double Solver::get_tau() const { + check_solver(impl); + return std::visit([](auto &ptr) { return ptr->get_tau(); }, *impl->solver); +} + +double Solver::get_kT() const { + check_solver(impl); + return std::visit([](auto &ptr) { return ptr->get_kT(); }, *impl->solver); +} + +Utils::VectorXd<9> Solver::get_pressure_tensor() const { + check_solver(impl); + return std::visit([](auto &ptr) { return ptr->get_pressure_tensor(); }, + *impl->solver); +} + +std::optional +Solver::get_interpolated_velocity(Utils::Vector3d const &pos) const { + /* calculate fluid velocity at particle's position + this is done by linear interpolation + (Eq. (11) Ahlrichs and Duenweg, JCP 111(17):8225 (1999)) */ + return std::visit( + [&](auto &ptr) { + auto const agrid = ptr->get_agrid(); + auto const lb_pos = folded_position(pos, ::box_geo) / agrid; + return ptr->get_velocity_at_pos(lb_pos, false); + }, + *impl->solver); +} + +std::optional +Solver::get_interpolated_density(Utils::Vector3d const &pos) const { + return std::visit( + [&](auto &ptr) { + auto const agrid = ptr->get_agrid(); + auto const lb_pos = folded_position(pos, ::box_geo) / agrid; + return ptr->get_density_at_pos(lb_pos, false); + }, + *impl->solver); +} + +Utils::Vector3d +Solver::get_coupling_interpolated_velocity(Utils::Vector3d const &pos) const { + return std::visit( + [&](auto &ptr) { + auto const agrid = ptr->get_agrid(); + auto const res = ptr->get_velocity_at_pos(pos / agrid, true); + assert(res); + return *res * (ptr->get_agrid() / ptr->get_tau()); + }, + *impl->solver); +} + +void Solver::add_force_density(Utils::Vector3d const &pos, + Utils::Vector3d const &force_density) { + std::visit( + [&](auto &ptr) { + if (not ptr->add_force_at_pos(pos / ptr->get_agrid(), force_density)) { + throw std::runtime_error("Cannot apply force to LB"); + } + }, + *impl->solver); +} + +Utils::Vector3d Solver::get_momentum() const { + check_solver(impl); + return std::visit([](auto const &ptr) { return ptr->get_momentum(); }, + *impl->solver); +} + +template <> void Solver::set(std::shared_ptr lb_instance) { + assert(impl); + assert(not impl->solver.has_value()); + impl->solver = lb_instance; +} + +#ifdef WALBERLA +template <> +void Solver::set(std::shared_ptr lb_fluid, + std::shared_ptr lb_params) { + assert(impl); + assert(not impl->solver.has_value()); + auto lb_instance = std::make_shared(lb_fluid, lb_params); + lb_instance->sanity_checks(); + auto const &lebc = ::box_geo.lees_edwards_bc(); + lb_fluid->check_lebc(lebc.shear_direction, lebc.shear_plane_normal); + impl->solver = lb_instance; +} +#endif // WALBERLA + +} // namespace LB diff --git a/src/core/lb/Solver.hpp b/src/core/lb/Solver.hpp new file mode 100644 index 00000000000..ef71e7a6154 --- /dev/null +++ b/src/core/lb/Solver.hpp @@ -0,0 +1,160 @@ +/* + * Copyright (C) 2023 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 . + */ + +#pragma once + +#include "utils.hpp" + +#include + +#include +#include +#include +#include + +namespace LB { + +struct Solver { + struct Implementation; + + Solver(); + + /** @brief Return true if a @c LB solver is active. */ + [[nodiscard]] bool is_solver_set() const; + + /** @brief Remove the @c LB solver. */ + void reset(); + + /** + * @brief Set the @c LB solver. + * For developers: a specialization must exist for every new @c LB type. + */ + template void set(Args... args); + + /** + * @brief Connector to the implementation internal state. + * For developers: use this mechanism to access the underlying variant. + */ + template void connect(Connector &&connector) const { + assert(impl != nullptr); + connector(*this->impl); + } + + /** + * @brief Propagate the @c LB fluid. + */ + void propagate(); + + /** + * @brief Perform a full initialization of the lattice-Boltzmann system. + * All derived parameters and the fluid are reset to their default values. + */ + void init() const {} + + /** + * @brief Perform @c LB parameter and boundary velocity checks. + */ + void sanity_checks() const; + + /** + * @brief Check if a MD time step is compatible with the @c LB tau. + */ + void veto_time_step(double time_step) const; + + /** + * @brief Perform @c LB LEbc parameter checks. + */ + void lebc_sanity_checks(unsigned int shear_direction, + unsigned int shear_plane_normal) const; + + /** + * @brief Get the @c LB time step. + */ + double get_tau() const; + + /** + * @brief Get the @c LB grid spacing. + */ + double get_agrid() const; + + /** + * @brief Get the thermal energy parameter of the @c LB fluid. + */ + double get_kT() const; + + /** + * @brief Get the lattice speed (agrid/tau). + */ + auto get_lattice_speed() const { return get_agrid() / get_tau(); } + + Utils::VectorXd<9> get_pressure_tensor() const; + + Utils::Vector3d get_momentum() const; + + /** + * @brief Calculate the interpolated fluid velocity in @c LB units. + * Use this function in MPI-parallel code. The @c LB ghost layer is ignored. + * @param pos Position in MD units at which the velocity is to be calculated. + * @retval interpolated fluid velocity. + */ + std::optional + get_interpolated_velocity(Utils::Vector3d const &pos) const; + + /** + * @brief Calculate the interpolated fluid density in @c LB units. + * Use this function in MPI-parallel code. The @c LB ghost layer is ignored. + * @param pos Position in MD units at which the density is to be calculated. + * @retval interpolated fluid density. + */ + std::optional + get_interpolated_density(Utils::Vector3d const &pos) const; + + /** + * @brief Calculate the interpolated fluid velocity in MD units. + * Special method used only for particle coupling. Uses the @c LB ghost layer. + * @param pos Position in MD units at which the velocity is to be calculated. + * @retval interpolated fluid velocity. + */ + Utils::Vector3d + get_coupling_interpolated_velocity(Utils::Vector3d const &pos) const; + + /** + * @brief Add a force density to the fluid at the given position. + * @param pos Position at which the force density is to be applied. + * @param force_density Force density to apply. + */ + void add_force_density(Utils::Vector3d const &pos, + Utils::Vector3d const &force_density); + + auto get_steps_per_md_step(double time_step) const { + return static_cast(std::round(get_tau() / time_step)); + } + + void on_boxl_change(); + void on_node_grid_change(); + void on_cell_structure_change(); + void on_timestep_change(); + void on_temperature_change(); + +private: + /** @brief Pointer-to-implementation. */ + std::unique_ptr impl; +}; + +} // namespace LB diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/lb/particle_coupling.cpp similarity index 79% rename from src/core/grid_based_algorithms/lb_particle_coupling.cpp rename to src/core/lb/particle_coupling.cpp index 153497bf576..5a03f99e2e2 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/lb/particle_coupling.cpp @@ -16,7 +16,9 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#include "LocalBox.hpp" + +#include "lb/particle_coupling.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "cells.hpp" #include "communication.hpp" @@ -24,16 +26,14 @@ #include "errorhandling.hpp" #include "grid.hpp" #include "random.hpp" +#include "system/System.hpp" #include "thermostat.hpp" -#include "grid_based_algorithms/lb_interface.hpp" -#include "grid_based_algorithms/lb_interpolation.hpp" -#include "grid_based_algorithms/lb_particle_coupling.hpp" - #include #include #include +#include #ifdef CALIPER #include @@ -48,6 +48,8 @@ LB::ParticleCouplingConfig lb_particle_coupling; +static auto is_lb_active() { return System::get_system().lb.is_solver_set(); } + void mpi_bcast_lb_particle_coupling_local() { boost::mpi::broadcast(comm_cart, lb_particle_coupling, 0); } @@ -61,8 +63,7 @@ void mpi_bcast_lb_particle_coupling() { void lb_lbcoupling_activate() { lb_particle_coupling.couple_to_md = true; } void lb_lbcoupling_deactivate() { - if (lattice_switch != ActiveLB::NONE && this_node == 0 && - lb_particle_coupling.gamma > 0.) { + if (is_lb_active() and this_node == 0 and lb_particle_coupling.gamma > 0.) { runtimeWarningMsg() << "Recalculating forces, so the LB coupling forces are not " "included in the particle force the first time step. This " @@ -79,7 +80,7 @@ void lb_lbcoupling_set_gamma(double gamma) { double lb_lbcoupling_get_gamma() { return lb_particle_coupling.gamma; } bool lb_lbcoupling_is_seed_required() { - if (lattice_switch == ActiveLB::WALBERLA_LB) { + if (is_lb_active()) { return not lb_particle_coupling.rng_counter_coupling.is_initialized(); } return false; @@ -90,35 +91,26 @@ uint64_t lb_coupling_get_rng_state_cpu() { } uint64_t lb_lbcoupling_get_rng_state() { - if (lattice_switch == ActiveLB::WALBERLA_LB) { + if (is_lb_active()) { return lb_coupling_get_rng_state_cpu(); } throw std::runtime_error("No LB active"); } void lb_lbcoupling_set_rng_state(uint64_t counter) { - if (lattice_switch == ActiveLB::WALBERLA_LB) { + if (is_lb_active()) { lb_particle_coupling.rng_counter_coupling = Utils::Counter(counter); } else throw std::runtime_error("No LB active"); } -void add_md_force(Utils::Vector3d const &pos, Utils::Vector3d const &force, - double time_step) { +void add_md_force(LB::Solver &lb, Utils::Vector3d const &pos, + Utils::Vector3d const &force, double time_step) { /* transform momentum transfer to lattice units (eq. (12) @cite ahlrichs99a) */ - auto const delta_j = (time_step / LB::get_lattice_speed()) * force; - lb_lbinterpolation_add_force_density(pos, delta_j); -} - -Utils::Vector3d lb_particle_coupling_drift_vel_offset(const Particle &p) { - Utils::Vector3d vel_offset{}; - -#ifdef LB_ELECTROHYDRODYNAMICS - vel_offset += p.mu_E(); -#endif - return vel_offset; + auto const delta_j = (time_step / lb.get_lattice_speed()) * force; + lb.add_force_density(pos, delta_j); } static Thermostat::GammaType lb_handle_particle_anisotropy(Particle const &p) { @@ -136,16 +128,13 @@ static Thermostat::GammaType lb_handle_particle_anisotropy(Particle const &p) { #endif // THERMOSTAT_PER_PARTICLE } -Utils::Vector3d lb_drag_force(Particle const &p, +Utils::Vector3d lb_drag_force(LB::Solver const &lb, Particle const &p, Utils::Vector3d const &shifted_pos, Utils::Vector3d const &vel_offset) { /* calculate fluid velocity at particle's position this is done by linear interpolation (eq. (11) @cite ahlrichs99a) */ - auto const interpolated_u = - lb_lbinterpolation_get_interpolated_velocity(shifted_pos) * - LB::get_lattice_speed(); - - Utils::Vector3d v_drift = interpolated_u + vel_offset; + auto const v_fluid = lb.get_coupling_interpolated_velocity(shifted_pos); + auto const v_drift = v_fluid + vel_offset; auto const gamma = lb_handle_particle_anisotropy(p); /* calculate viscous force (eq. (9) @cite ahlrichs99a) */ @@ -160,7 +149,7 @@ Utils::Vector3d lb_drag_force(Particle const &p, * * @return True iff the point is inside of the box up to halo. */ -inline bool in_local_domain(Utils::Vector3d const &pos, double halo = 0.) { +static bool in_local_domain(Utils::Vector3d const &pos, double halo = 0.) { auto const halo_vec = Utils::Vector3d::broadcast(halo); auto const lower_corner = local_geo.my_left() - halo_vec; auto const upper_corner = local_geo.my_right() + halo_vec; @@ -168,18 +157,21 @@ inline bool in_local_domain(Utils::Vector3d const &pos, double halo = 0.) { return pos >= lower_corner and pos < upper_corner; } -bool in_local_halo(Utils::Vector3d const &pos) { - auto const halo = 0.5 * LB::get_agrid(); - +static bool in_local_halo(Utils::Vector3d const &pos, double agrid) { + auto const halo = 0.5 * agrid; return in_local_domain(pos, halo); } +bool in_local_halo(Utils::Vector3d const &pos) { + return in_local_domain(pos, System::get_system().lb.get_agrid()); +} + /** * @brief Return a vector of positions shifted by +,- box length in each * coordinate */ -std::vector positions_in_halo(Utils::Vector3d pos, - const BoxGeometry &box) { +std::vector +positions_in_halo(Utils::Vector3d pos, BoxGeometry const &box, double agrid) { std::vector res; for (int i : {-1, 0, 1}) { for (int j : {-1, 0, 1}) { @@ -197,7 +189,7 @@ std::vector positions_in_halo(Utils::Vector3d pos, pos_shifted[le.shear_direction] -= le.pos_offset; } - if (in_local_halo(pos_shifted)) { + if (in_local_halo(pos_shifted, agrid)) { res.push_back(pos_shifted); } } @@ -231,16 +223,18 @@ void ParticleCoupling::kernel(Particle &p) { if (p.is_virtual() and not m_couple_virtual) return; + auto const agrid = m_lb.get_agrid(); + // Calculate coupling force Utils::Vector3d force_on_particle = {}; #ifdef ENGINE if (not p.swimming().is_engine_force_on_fluid) #endif - for (auto pos : positions_in_halo(p.pos(), box_geo)) { - if (in_local_halo(pos)) { - auto const vel_offset = lb_particle_coupling_drift_vel_offset(p); - auto const drag_force = lb_drag_force(p, pos, vel_offset); + for (auto pos : positions_in_halo(p.pos(), box_geo, agrid)) { + if (in_local_halo(pos, agrid)) { + auto const vel_offset = lb_drift_velocity_offset(p); + auto const drag_force = lb_drag_force(m_lb, p, pos, vel_offset); auto const random_force = get_noise_term(p); force_on_particle = drag_force + random_force; break; @@ -256,13 +250,13 @@ void ParticleCoupling::kernel(Particle &p) { // couple positions including shifts by one box length to add // forces to ghost layers - for (auto pos : positions_in_halo(p.pos(), box_geo)) { + for (auto pos : positions_in_halo(p.pos(), box_geo, agrid)) { if (in_local_domain(pos)) { /* Particle is in our LB volume, so this node * is responsible to adding its force */ p.force() += force_on_particle; } - add_md_force(pos, force_on_fluid, m_time_step); + add_md_force(m_lb, pos, force_on_fluid, m_time_step); } } @@ -289,9 +283,10 @@ void couple_particles(bool couple_virtual, ParticleRange const &real_particles, #ifdef CALIPER CALI_CXX_MARK_FUNCTION; #endif - if (lattice_switch == ActiveLB::WALBERLA_LB) { - if (lb_particle_coupling.couple_to_md) { - ParticleCoupling coupling{couple_virtual, time_step}; + if (lb_particle_coupling.couple_to_md) { + auto &lb = System::get_system().lb; + if (lb.is_solver_set()) { + ParticleCoupling coupling{lb, couple_virtual, time_step}; CouplingBookkeeping bookkeeping{}; for (auto const &particle_range : {real_particles, ghost_particles}) { for (auto &p : particle_range) { @@ -310,11 +305,8 @@ void couple_particles(bool couple_virtual, ParticleRange const &real_particles, } // namespace LB void lb_lbcoupling_propagate() { - if (lattice_switch != ActiveLB::NONE) { - if (LB::get_kT() > 0.0) { - if (lattice_switch == ActiveLB::WALBERLA_LB) { - lb_particle_coupling.rng_counter_coupling->increment(); - } - } + auto const &lb = System::get_system().lb; + if (lb.is_solver_set() and lb.get_kT() > 0.0) { + lb_particle_coupling.rng_counter_coupling->increment(); } } diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.hpp b/src/core/lb/particle_coupling.hpp similarity index 77% rename from src/core/grid_based_algorithms/lb_particle_coupling.hpp rename to src/core/lb/particle_coupling.hpp index f0b0a6d2821..1486ef20605 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.hpp +++ b/src/core/lb/particle_coupling.hpp @@ -16,12 +16,13 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef LB_PARTICLE_COUPLING_HPP -#define LB_PARTICLE_COUPLING_HPP +#pragma once + +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" -#include "grid.hpp" +#include "lb/Solver.hpp" #include #include @@ -29,7 +30,6 @@ #include #include -#include #include #include @@ -71,40 +71,36 @@ bool in_local_halo(Utils::Vector3d const &pos); /** * @brief Add a force to the lattice force density. - * @param pos Position of the force + * @param lb Hydrodynamics solver + * @param pos Position of the force in LB units. * @param force Force in MD units. * @param time_step MD time step. */ -void add_md_force(Utils::Vector3d const &pos, Utils::Vector3d const &force, - double time_step); +void add_md_force(LB::Solver &lb, Utils::Vector3d const &pos, + Utils::Vector3d const &force, double time_step); // internal function exposed for unit testing -std::vector positions_in_halo(Utils::Vector3d pos, - const BoxGeometry &box); +std::vector +positions_in_halo(Utils::Vector3d pos, BoxGeometry const &box, double agrid); // internal function exposed for unit testing void add_swimmer_force(Particle const &p, double time_step); -/** - * @brief Calculate particle drift velocity offset due to ENGINE and - * ELECTROHYDRODYNAMICS. - */ -Utils::Vector3d lb_particle_coupling_drift_vel_offset(const Particle &p); - void mpi_bcast_lb_particle_coupling(); /** @brief Calculate drag force on a single particle. * * See section II.C. @cite ahlrichs99a * + * @param[in] lb The coupled fluid * @param[in] p The coupled particle - * @param[in] shifted_pos The particle position with optional shift - * @param[in] vel_offset Velocity offset to be added to interpolated LB - * velocity before calculating the force + * @param[in] shifted_pos The particle position in LB units with optional shift + * @param[in] vel_offset Velocity offset in MD units to be added to + * interpolated LB velocity before calculating the force * * @return The viscous coupling force */ -Utils::Vector3d lb_drag_force(Particle const &p, +Utils::Vector3d lb_drag_force(LB::Solver const &lb, Particle const &p, Utils::Vector3d const &shifted_pos, Utils::Vector3d const &vel_offset); @@ -136,14 +132,16 @@ void couple_particles(bool couple_virtual, ParticleRange const &real_particles, ParticleRange const &ghost_particles, double time_step); class ParticleCoupling { + LB::Solver &m_lb; bool m_couple_virtual; bool m_thermalized; double m_time_step; double m_noise_pref_wo_gamma; public: - ParticleCoupling(bool couple_virtual, double time_step, double kT) - : m_couple_virtual{couple_virtual}, m_thermalized{kT != 0.}, + ParticleCoupling(LB::Solver &lb, bool couple_virtual, double time_step, + double kT) + : m_lb{lb}, m_couple_virtual{couple_virtual}, m_thermalized{kT != 0.}, m_time_step{time_step} { assert(kT >= 0.); /* Eq. (16) @cite ahlrichs99a, without the gamma term. @@ -155,12 +153,24 @@ class ParticleCoupling { m_noise_pref_wo_gamma = std::sqrt(variance_inv * 2. * kT / time_step); } - ParticleCoupling(bool couple_virtual, double time_step) - : ParticleCoupling(couple_virtual, time_step, - LB::get_kT() * Utils::sqr(LB::get_lattice_speed())) {} + ParticleCoupling(LB::Solver &lb, bool couple_virtual, double time_step) + : ParticleCoupling(lb, couple_virtual, time_step, + lb.get_kT() * Utils::sqr(lb.get_lattice_speed())) {} Utils::Vector3d get_noise_term(Particle const &p) const; void kernel(Particle &p); + + /** + * @brief Calculate particle drift velocity offset due to ENGINE and + * ELECTROHYDRODYNAMICS. + */ + auto lb_drift_velocity_offset(Particle const &p) const { + Utils::Vector3d vel_offset{}; +#ifdef LB_ELECTROHYDRODYNAMICS + vel_offset += p.mu_E(); +#endif + return vel_offset; + } }; /** @@ -192,5 +202,3 @@ class CouplingBookkeeping { }; } // namespace LB - -#endif diff --git a/src/core/lb/utils.hpp b/src/core/lb/utils.hpp new file mode 100644 index 00000000000..fc793910f63 --- /dev/null +++ b/src/core/lb/utils.hpp @@ -0,0 +1,30 @@ +/* + * Copyright (C) 2023 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 . + */ + +#pragma once + +#include + +namespace LB { + +struct NoLBActive : public std::exception { + const char *what() const noexcept override { return "LB not activated"; } +}; + +} // namespace LB diff --git a/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp b/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp index a8360630999..9c4b06e8ed5 100644 --- a/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp +++ b/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp @@ -20,7 +20,7 @@ #include "BoxGeometry.hpp" #include "grid.hpp" -#include "grid_based_algorithms/lb_interface.hpp" +#include "system/System.hpp" #include "utils_histogram.hpp" #include @@ -46,13 +46,14 @@ CylindricalLBFluxDensityProfileAtParticlePositions::evaluate( local_folded_positions.reserve(local_particles.size()); local_flux_densities.reserve(local_particles.size()); - auto const vel_conv = LB::get_lattice_speed(); + auto const &lb = System::get_system().lb; + auto const vel_conv = lb.get_lattice_speed(); for (auto const &p : local_particles) { auto const pos = folded_position(traits.position(p), box_geo); auto const pos_shifted = pos - transform_params->center(); - auto const vel = *(LB::get_interpolated_velocity(pos)); - auto const dens = *(LB::get_interpolated_density(pos)); + auto const vel = *lb.get_interpolated_velocity(pos); + auto const dens = *lb.get_interpolated_density(pos); auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder( pos_shifted, transform_params->axis(), transform_params->orientation()); auto const flux_cyl = Utils::transform_vector_cartesian_to_cylinder( diff --git a/src/core/observables/CylindricalLBVelocityProfile.cpp b/src/core/observables/CylindricalLBVelocityProfile.cpp index 6a12d4285ef..e2606321975 100644 --- a/src/core/observables/CylindricalLBVelocityProfile.cpp +++ b/src/core/observables/CylindricalLBVelocityProfile.cpp @@ -19,7 +19,7 @@ #include "CylindricalLBVelocityProfile.hpp" -#include "grid_based_algorithms/lb_interface.hpp" +#include "system/System.hpp" #include "utils_histogram.hpp" #include @@ -36,10 +36,11 @@ std::vector CylindricalLBVelocityProfile::operator()( decltype(sampling_positions) local_positions{}; std::vector local_velocities{}; - auto const vel_conv = LB::get_lattice_speed(); + auto const &lb = System::get_system().lb; + auto const vel_conv = lb.get_lattice_speed(); for (auto const &pos : sampling_positions) { - if (auto const vel = LB::get_interpolated_velocity(pos)) { + if (auto const vel = lb.get_interpolated_velocity(pos)) { auto const pos_shifted = pos - transform_params->center(); auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder( pos_shifted, transform_params->axis(), diff --git a/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp b/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp index f498b168f07..cce02d6f13e 100644 --- a/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp +++ b/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp @@ -20,7 +20,7 @@ #include "BoxGeometry.hpp" #include "grid.hpp" -#include "grid_based_algorithms/lb_interface.hpp" +#include "system/System.hpp" #include "utils_histogram.hpp" #include @@ -42,12 +42,13 @@ std::vector CylindricalLBVelocityProfileAtParticlePositions::evaluate( local_folded_positions.reserve(local_particles.size()); local_velocities.reserve(local_particles.size()); - auto const vel_conv = LB::get_lattice_speed(); + auto const &lb = System::get_system().lb; + auto const vel_conv = lb.get_lattice_speed(); for (auto const &p : local_particles) { auto const pos = folded_position(traits.position(p), box_geo); auto const pos_shifted = pos - transform_params->center(); - auto const vel = *(LB::get_interpolated_velocity(pos)); + auto const vel = *lb.get_interpolated_velocity(pos); auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder( pos_shifted, transform_params->axis(), transform_params->orientation()); auto const vel_cyl = Utils::transform_vector_cartesian_to_cylinder( diff --git a/src/core/observables/LBFluidPressureTensor.hpp b/src/core/observables/LBFluidPressureTensor.hpp index cafeaa643d2..3a7fe78bcf5 100644 --- a/src/core/observables/LBFluidPressureTensor.hpp +++ b/src/core/observables/LBFluidPressureTensor.hpp @@ -20,11 +20,15 @@ #define OBSERVABLES_LB_FLUID_STRESS_HPP #include "Observable.hpp" -#include "grid_based_algorithms/lb_interface.hpp" +#include "system/System.hpp" #include +#include + +#include #include +#include #include namespace Observables { @@ -33,9 +37,11 @@ class LBFluidPressureTensor : public Observable { std::vector shape() const override { return {3, 3}; } std::vector operator()(boost::mpi::communicator const &comm) const override { - auto const unit_conversion = - 1. / (LB::get_agrid() * Utils::sqr(LB::get_tau())); - auto const tensor = LB::get_pressure_tensor(comm) * unit_conversion; + auto const &lb = System::get_system().lb; + auto const pressure_conv = 1. / (lb.get_agrid() * Utils::sqr(lb.get_tau())); + auto const local_tensor = lb.get_pressure_tensor() * pressure_conv; + std::remove_const_t tensor; + boost::mpi::reduce(comm, local_tensor, tensor, std::plus<>(), 0); return tensor.as_vector(); } }; diff --git a/src/core/observables/LBVelocityProfile.cpp b/src/core/observables/LBVelocityProfile.cpp index 706d824afee..13db451211c 100644 --- a/src/core/observables/LBVelocityProfile.cpp +++ b/src/core/observables/LBVelocityProfile.cpp @@ -18,7 +18,7 @@ */ #include "LBVelocityProfile.hpp" -#include "grid_based_algorithms/lb_interface.hpp" +#include "system/System.hpp" #include "utils_histogram.hpp" #include @@ -36,10 +36,11 @@ LBVelocityProfile::operator()(boost::mpi::communicator const &comm) const { decltype(sampling_positions) local_positions{}; std::vector local_velocities{}; - auto const vel_conv = LB::get_lattice_speed(); + auto const &lb = System::get_system().lb; + auto const vel_conv = lb.get_lattice_speed(); for (auto const &pos : sampling_positions) { - if (auto const vel = LB::get_interpolated_velocity(pos)) { + if (auto const vel = lb.get_interpolated_velocity(pos)) { local_positions.emplace_back(pos); local_velocities.emplace_back((*vel) * vel_conv); } diff --git a/src/core/system/System.hpp b/src/core/system/System.hpp index 2b999f197c2..0b4ad5bcb14 100644 --- a/src/core/system/System.hpp +++ b/src/core/system/System.hpp @@ -27,9 +27,11 @@ #include "electrostatics/solver.hpp" #include "magnetostatics/solver.hpp" +#include "ek/Solver.hpp" +#include "lb/Solver.hpp" + #include -#include #include namespace System { @@ -49,6 +51,8 @@ class System { } Coulomb::Solver coulomb; Dipoles::Solver dipoles; + LB::Solver lb; + EK::Solver ek; }; System &get_system(); diff --git a/src/core/system/System.impl.hpp b/src/core/system/System.impl.hpp index 2af2e24eea6..03b364d525d 100644 --- a/src/core/system/System.impl.hpp +++ b/src/core/system/System.impl.hpp @@ -21,3 +21,6 @@ #include "electrostatics/coulomb.hpp" #include "magnetostatics/dipoles.hpp" + +#include "ek/Implementation.hpp" +#include "lb/Implementation.hpp" diff --git a/src/core/unit_tests/ek_interface_test.cpp b/src/core/unit_tests/ek_interface_test.cpp index 1d01ef76b59..0a938ba5b05 100644 --- a/src/core/unit_tests/ek_interface_test.cpp +++ b/src/core/unit_tests/ek_interface_test.cpp @@ -17,16 +17,31 @@ * along with this program. If not, see . */ -#define BOOST_TEST_MODULE LB particle coupling test +#define BOOST_TEST_MODULE EK interface test #define BOOST_TEST_DYN_LINK #define BOOST_TEST_NO_MAIN #include +#include "ParticleFactory.hpp" + #include "EspressoSystemStandAlone.hpp" #include "config/config.hpp" +#include "ek/EKReactions.hpp" +#include "ek/EKWalberla.hpp" +#include "ek/Implementation.hpp" #include "errorhandling.hpp" -#include "grid_based_algorithms/ek_container.hpp" -#include "grid_based_algorithms/ek_reactions.hpp" +#include "grid.hpp" +#include "system/System.hpp" + +#ifdef WALBERLA +#include +#include +#include +#include +#include +#include +#include +#endif #include @@ -55,82 +70,161 @@ static struct { namespace espresso { // ESPResSo system instance static std::unique_ptr system; -} // namespace espresso +// ESPResSo actors +#ifdef WALBERLA +static std::shared_ptr ek_container; +static std::shared_ptr ek_reactions; +static std::shared_ptr ek_instance; +static std::shared_ptr ek_lattice; +#endif -static auto get_n_runtime_errors() { return check_runtime_errors_local(); } +static auto make_ek_actor() { +#ifdef WALBERLA + auto constexpr n_ghost_layers = 1u; + auto constexpr single_precision = true; + ek_lattice = std::make_shared(params.grid_dimensions, + ::node_grid, n_ghost_layers); + ek_container = std::make_shared( + params.tau, new_ek_poisson_none(ek_lattice, single_precision)); + ek_reactions = std::make_shared(); + ek_instance = std::make_shared(ek_container, ek_reactions); +#endif +} +static void add_ek_actor() { #ifdef WALBERLA + System::get_system().ek.set<::EK::EKWalberla>(ek_instance); +#endif +} -#include "grid.hpp" +static void remove_ek_actor() { System::get_system().ek.reset(); } +} // namespace espresso -#include -#include -#include -#include -#include +#ifdef WALBERLA +namespace walberla { +class EKReactionImpl : public EKReactionBase { +public: + using EKReactionBase::EKReactionBase; + using EKReactionBase::get_coefficient; + using EKReactionBase::get_lattice; + using EKReactionBase::get_reactants; + + void perform_reaction() override {} + ~EKReactionImpl() override = default; +}; +} // namespace walberla +#endif // WALBERLA +/** Fixture to manage the lifetime of the EK actor. */ +struct CleanupActorEK : public ParticleFactory { + CleanupActorEK() : ParticleFactory() { + espresso::make_ek_actor(); + espresso::add_ek_actor(); + } + + // NOLINTNEXTLINE(bugprone-exception-escape) + ~CleanupActorEK() { espresso::remove_ek_actor(); } +}; + +BOOST_FIXTURE_TEST_SUITE(suite, CleanupActorEK) + +static auto get_n_runtime_errors() { return check_runtime_errors_local(); } + +#ifdef WALBERLA BOOST_AUTO_TEST_CASE(ek_interface_walberla) { + auto &ek = System::get_system().ek; + { // tau setters and getters - BOOST_CHECK_EQUAL(EK::ek_container.get_tau(), 0.); - BOOST_CHECK_EQUAL(EK::get_tau(), 0.); - BOOST_CHECK_EQUAL(EK::get_steps_per_md_step(1.), 0); - EK::ek_container.set_tau(2.); - BOOST_CHECK_EQUAL(EK::ek_container.get_tau(), 2.); - BOOST_CHECK_EQUAL(EK::get_tau(), 2.); - BOOST_CHECK_EQUAL(EK::get_steps_per_md_step(1.), 2); - BOOST_CHECK_EQUAL(EK::get_steps_per_md_step(2.), 1); - BOOST_CHECK_EQUAL(EK::get_steps_per_md_step(5.), 0); + espresso::ek_container->set_tau(2.); + BOOST_CHECK_EQUAL(ek.get_tau(), 2.); + BOOST_CHECK_EQUAL(ek.get_steps_per_md_step(1.), 2); + BOOST_CHECK_EQUAL(ek.get_steps_per_md_step(2.), 1); + BOOST_CHECK_EQUAL(ek.get_steps_per_md_step(5.), 0); } { // setup a minimal EK model without coupling to LB - auto constexpr n_ghost_layers = 1u; + using walberla::EKReactant; auto constexpr single_precision = true; - auto ek_lattice = std::make_shared( - params.grid_dimensions, ::node_grid, n_ghost_layers); + auto constexpr stoich = 1.; + auto constexpr order = 2.; auto ek_species = new_ek_walberla( - ek_lattice, params.diffusion, params.kT, params.valency, + espresso::ek_lattice, params.diffusion, params.kT, params.valency, params.ext_efield, params.density, false, false, single_precision); - auto ek_solver_none = new_ek_poisson_none(ek_lattice, single_precision); - - BOOST_REQUIRE(EK::ek_reactions.empty()); - BOOST_REQUIRE(EK::ek_container.empty()); - BOOST_REQUIRE(not EK::ek_container.is_poisson_solver_set()); - EK::propagate(); // no-op + auto ek_reactant = std::make_shared(ek_species, stoich, order); + auto ek_reaction = std::make_shared( + espresso::ek_lattice, + std::vector>{ + {ek_reactant, ek_reactant, ek_reactant}}, + 1.); + ek_reaction->perform_reaction(); + + BOOST_REQUIRE(espresso::ek_reactions->empty()); + BOOST_REQUIRE(espresso::ek_container->empty()); + ek.propagate(); // no-op BOOST_REQUIRE_EQUAL(get_n_runtime_errors(), 0); - EK::ek_container.set_poisson_solver(ek_solver_none); - BOOST_REQUIRE(EK::ek_container.is_poisson_solver_set()); - BOOST_REQUIRE(EK::ek_container.empty()); - EK::ek_container.set_tau(0.); - BOOST_CHECK_THROW(EK::ek_container.add(ek_species), std::runtime_error); - EK::ek_container.set_tau(2.); - EK::ek_container.add(ek_species); - BOOST_REQUIRE(not EK::ek_container.empty()); - EK::propagate(); // no-op + BOOST_REQUIRE(espresso::ek_container->empty()); + BOOST_REQUIRE(not espresso::ek_container->contains(ek_species)); + BOOST_REQUIRE(not espresso::ek_reactions->contains(ek_reaction)); + espresso::ek_container->set_tau(2.); + espresso::ek_container->add(ek_species); + BOOST_REQUIRE(not espresso::ek_container->empty()); + BOOST_REQUIRE(espresso::ek_container->contains(ek_species)); + ek.propagate(); // no-op BOOST_REQUIRE_EQUAL(get_n_runtime_errors(), 0); - EK::ek_container.remove(ek_species); - BOOST_REQUIRE(EK::ek_container.empty()); - EK::propagate(); // no-op + espresso::ek_reactions->add(ek_reaction); + BOOST_REQUIRE(espresso::ek_reactions->contains(ek_reaction)); + BOOST_REQUIRE(not espresso::ek_reactions->empty()); + espresso::ek_reactions->remove(ek_reaction); + BOOST_REQUIRE(not espresso::ek_reactions->contains(ek_reaction)); + BOOST_REQUIRE(espresso::ek_reactions->empty()); + espresso::ek_container->remove(ek_species); + BOOST_REQUIRE(espresso::ek_container->empty()); + BOOST_REQUIRE(not espresso::ek_container->contains(ek_species)); + ek.propagate(); // no-op BOOST_REQUIRE_EQUAL(get_n_runtime_errors(), 0); } + + { + // EK prevents changing most of the system state + BOOST_CHECK_THROW(ek.on_boxl_change(), std::runtime_error); + BOOST_CHECK_THROW(ek.on_timestep_change(), std::runtime_error); + BOOST_CHECK_THROW(ek.on_temperature_change(), std::runtime_error); + BOOST_CHECK_THROW(ek.on_node_grid_change(), std::runtime_error); + } } +#endif // WALBERLA -#else // WALBERLA +BOOST_AUTO_TEST_CASE(ek_interface_none) { + auto &ek = System::get_system().ek; -BOOST_AUTO_TEST_CASE(ek_interface) { { - EK::propagate(); // no-op - BOOST_CHECK_THROW(EK::get_tau(), NoEKActive); - BOOST_CHECK_THROW(EK::get_tau(), std::exception); - BOOST_CHECK_THROW(EK::get_steps_per_md_step(1.), std::exception); + using EK::NoEKActive; + ek.reset(); + BOOST_CHECK_THROW(ek.get_tau(), NoEKActive); + BOOST_CHECK_THROW(ek.propagate(), NoEKActive); + auto ek_impl = std::make_shared(); + ek.set(ek_impl); + BOOST_CHECK_THROW(ek.is_ready_for_propagation(), NoEKActive); + BOOST_CHECK_THROW(ek.propagate(), NoEKActive); + BOOST_CHECK_THROW(ek.get_tau(), NoEKActive); + BOOST_CHECK_THROW(ek.sanity_checks(), NoEKActive); + BOOST_CHECK_THROW(ek.veto_time_step(0.), NoEKActive); + BOOST_CHECK_THROW(ek.on_cell_structure_change(), NoEKActive); + BOOST_CHECK_THROW(ek.on_boxl_change(), NoEKActive); + BOOST_CHECK_THROW(ek.on_node_grid_change(), NoEKActive); + BOOST_CHECK_THROW(ek.on_timestep_change(), NoEKActive); + BOOST_CHECK_THROW(ek.on_temperature_change(), NoEKActive); + BOOST_CHECK_THROW(ek.get_steps_per_md_step(1.), std::exception); auto const err_msg = std::string(NoEKActive().what()); auto const ref_msg = std::string("EK not activated"); BOOST_CHECK_EQUAL(err_msg, ref_msg); + ek.reset(); } } -#endif // WALBERLA +BOOST_AUTO_TEST_SUITE_END() int main(int argc, char **argv) { espresso::system = std::make_unique(argc, argv); diff --git a/src/core/unit_tests/lb_particle_coupling_test.cpp b/src/core/unit_tests/lb_particle_coupling_test.cpp index 25378093921..1f5d2456295 100644 --- a/src/core/unit_tests/lb_particle_coupling_test.cpp +++ b/src/core/unit_tests/lb_particle_coupling_test.cpp @@ -39,12 +39,12 @@ namespace utf = boost::unit_test; #include "errorhandling.hpp" #include "event.hpp" #include "grid.hpp" -#include "grid_based_algorithms/lb_interface.hpp" -#include "grid_based_algorithms/lb_interpolation.hpp" -#include "grid_based_algorithms/lb_particle_coupling.hpp" -#include "grid_based_algorithms/lb_walberla_instance.hpp" +#include "lb/LBNone.hpp" +#include "lb/LBWalberla.hpp" +#include "lb/particle_coupling.hpp" #include "particle_node.hpp" #include "random.hpp" +#include "system/System.hpp" #include "thermostat.hpp" #include @@ -55,6 +55,7 @@ namespace utf = boost::unit_test; #include #include +#include #include #include @@ -92,14 +93,14 @@ namespace espresso { // ESPResSo system instance static std::unique_ptr system; // ESPResSo actors -static std::shared_ptr lb_params; +static std::shared_ptr lb_params; static std::shared_ptr lb_lattice; static std::shared_ptr lb_fluid; static auto make_lb_actor() { auto constexpr n_ghost_layers = 1u; auto constexpr single_precision = false; - lb_params = std::make_shared(params.agrid, params.tau); + lb_params = std::make_shared(params.agrid, params.tau); lb_lattice = std::make_shared(params.grid_dimensions, ::node_grid, n_ghost_layers); lb_fluid = new_lb_walberla(lb_lattice, params.viscosity, params.density, @@ -108,9 +109,11 @@ static auto make_lb_actor() { lb_fluid->ghost_communication(); } -static void add_lb_actor() { activate_lb_walberla(lb_fluid, lb_params); } +static void add_lb_actor() { + System::get_system().lb.set<::LB::LBWalberla>(lb_fluid, lb_params); +} -static void remove_lb_actor() { deactivate_lb_walberla(); } +static void remove_lb_actor() { System::get_system().lb.reset(); } static void set_lb_kT(double kT) { lb_fluid->set_collision_model(kT, params.seed); @@ -178,8 +181,9 @@ BOOST_AUTO_TEST_CASE(rng_initial_state) { BOOST_AUTO_TEST_CASE(rng) { lb_lbcoupling_set_rng_state(17); lb_particle_coupling.gamma = 0.2; + auto &lb = System::get_system().lb; - LB::ParticleCoupling coupling{true, params.time_step, 1.}; + LB::ParticleCoupling coupling{lb, true, params.time_step, 1.}; BOOST_REQUIRE(lb_particle_coupling.rng_counter_coupling); BOOST_CHECK_EQUAL(lb_lbcoupling_get_rng_state(), 17); BOOST_CHECK(not lb_lbcoupling_is_seed_required()); @@ -205,43 +209,38 @@ BOOST_AUTO_TEST_CASE(rng) { BOOST_CHECK(step1_random1 != step2_random1); BOOST_CHECK(step1_random1 != step2_random2); - LB::ParticleCoupling coupling_unthermalized{true, params.time_step, 0.}; + LB::ParticleCoupling coupling_unthermalized{lb, true, params.time_step, 0.}; auto const step3_norandom = coupling_unthermalized.get_noise_term(test_partcl_2); BOOST_CHECK((step3_norandom == Utils::Vector3d{0., 0., 0.})); } -BOOST_AUTO_TEST_CASE(access_outside_domain) { - auto const invalid_pos = 2 * params.box_dimensions; - BOOST_CHECK_THROW(lb_lbinterpolation_get_interpolated_velocity(invalid_pos), - std::runtime_error); - BOOST_CHECK_THROW(lb_lbinterpolation_add_force_density(invalid_pos, {}), - std::runtime_error); -} - BOOST_AUTO_TEST_CASE(drift_vel_offset) { Particle p{}; - BOOST_CHECK_EQUAL(lb_particle_coupling_drift_vel_offset(p).norm(), 0); + auto &lb = System::get_system().lb; + LB::ParticleCoupling coupling{lb, false, params.time_step}; + BOOST_CHECK_EQUAL(coupling.lb_drift_velocity_offset(p).norm(), 0); Utils::Vector3d expected{}; #ifdef LB_ELECTROHYDRODYNAMICS p.mu_E() = Utils::Vector3d{-2., 1.5, 1.}; expected += p.mu_E(); #endif - BOOST_CHECK_SMALL( - (lb_particle_coupling_drift_vel_offset(p) - expected).norm(), eps); + BOOST_CHECK_SMALL((coupling.lb_drift_velocity_offset(p) - expected).norm(), + eps); } BOOST_DATA_TEST_CASE(drag_force, bdata::make(kTs), kT) { espresso::set_lb_kT(kT); + auto &lb = System::get_system().lb; Particle p{}; p.v() = {-2.5, 1.5, 2.}; - p.pos() = lb_walberla()->get_lattice().get_local_domain().first; + p.pos() = espresso::lb_fluid->get_lattice().get_local_domain().first; lb_lbcoupling_set_gamma(0.2); Utils::Vector3d drift_offset{-1., 1., 1.}; // Drag force in quiescent fluid { - auto const observed = lb_drag_force(p, p.pos(), drift_offset); + auto const observed = lb_drag_force(lb, p, p.pos(), drift_offset); const Utils::Vector3d expected{0.3, -0.1, -.2}; BOOST_CHECK_SMALL((observed - expected).norm(), eps); } @@ -250,8 +249,9 @@ BOOST_DATA_TEST_CASE(drag_force, bdata::make(kTs), kT) { #ifdef ENGINE BOOST_DATA_TEST_CASE(swimmer_force, bdata::make(kTs), kT) { espresso::set_lb_kT(kT); + auto &lb = System::get_system().lb; auto const first_lb_node = - lb_walberla()->get_lattice().get_local_domain().first; + espresso::lb_fluid->get_lattice().get_local_domain().first; Particle p{}; p.swimming().swimming = true; p.swimming().f_swim = 2.; @@ -261,7 +261,7 @@ BOOST_DATA_TEST_CASE(swimmer_force, bdata::make(kTs), kT) { // swimmer coupling { if (in_local_halo(p.pos())) { - LB::ParticleCoupling coupling{true, params.time_step}; + LB::ParticleCoupling coupling{lb, true, params.time_step}; coupling.kernel(p); auto const interpolated = LB::get_force_to_be_applied(p.pos()); auto const expected = @@ -294,7 +294,7 @@ BOOST_DATA_TEST_CASE(swimmer_force, bdata::make(kTs), kT) { // remove force of the particle from the fluid { if (in_local_halo(p.pos())) { - add_md_force(p.pos(), -Utils::Vector3d{0., 0., p.swimming().f_swim}, + add_md_force(lb, p.pos(), -Utils::Vector3d{0., 0., p.swimming().f_swim}, params.time_step); auto const reset = LB::get_force_to_be_applied(p.pos()); BOOST_REQUIRE_SMALL(reset.norm(), eps); @@ -306,12 +306,13 @@ BOOST_DATA_TEST_CASE(swimmer_force, bdata::make(kTs), kT) { BOOST_DATA_TEST_CASE(particle_coupling, bdata::make(kTs), kT) { espresso::set_lb_kT(kT); lb_lbcoupling_set_rng_state(17); + auto &lb = System::get_system().lb; auto const first_lb_node = - lb_walberla()->get_lattice().get_local_domain().first; + espresso::lb_fluid->get_lattice().get_local_domain().first; auto const gamma = 0.2; lb_lbcoupling_set_gamma(gamma); Particle p{}; - LB::ParticleCoupling coupling{false, params.time_step}; + LB::ParticleCoupling coupling{lb, false, params.time_step}; auto expected = coupling.get_noise_term(p); #ifdef LB_ELECTROHYDRODYNAMICS p.mu_E() = Utils::Vector3d{-2., 1.5, 1.}; @@ -334,7 +335,7 @@ BOOST_DATA_TEST_CASE(particle_coupling, bdata::make(kTs), kT) { // remove force of the particle from the fluid { if (in_local_halo(p.pos())) { - add_md_force(p.pos(), -expected, params.time_step); + add_md_force(lb, p.pos(), -expected, params.time_step); } } } @@ -345,8 +346,9 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, auto const rank = comm.rank(); espresso::set_lb_kT(kT); lb_lbcoupling_set_rng_state(17); + auto &lb = System::get_system().lb; auto const first_lb_node = - lb_walberla()->get_lattice().get_local_domain().first; + espresso::lb_fluid->get_lattice().get_local_domain().first; auto const gamma = 0.2; auto const pid = 0; auto const skin = params.skin; @@ -366,7 +368,7 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, set_particle_property(pid, &Particle::mu_E, Utils::Vector3d{-2., 1.5, 1.}); #endif - LB::ParticleCoupling coupling{thermo_virtual, params.time_step}; + LB::ParticleCoupling coupling{lb, thermo_virtual, params.time_step}; auto const p_opt = copy_particle_to_head_node(comm, pid); auto expected = Utils::Vector3d{}; if (rank == 0) { @@ -408,8 +410,9 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, boost::mpi::communicator world; assert(world.size() <= 4); auto const cutoff = 8 / world.size(); + auto const agrid = params.agrid; { - auto const shifts = positions_in_halo({0., 0., 0.}, box_geo); + auto const shifts = positions_in_halo({0., 0., 0.}, box_geo, agrid); BOOST_REQUIRE_EQUAL(shifts.size(), cutoff); for (std::size_t i = 0; i < shifts.size(); ++i) { BOOST_REQUIRE_EQUAL(shifts[i], reference_shifts[i]); @@ -417,14 +420,14 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, } { auto const reference_shift = Utils::Vector3d{1., 1., 1.}; - auto const shifts = positions_in_halo({1., 1., 1.}, box_geo); + auto const shifts = positions_in_halo({1., 1., 1.}, box_geo, agrid); BOOST_REQUIRE_EQUAL(shifts.size(), 1); BOOST_REQUIRE_EQUAL(shifts[0], reference_shift); } { auto const reference_origin = Utils::Vector3d{1., 2., 0.}; auto const reference_shift = Utils::Vector3d{1., 2., 8.}; - auto const shifts = positions_in_halo({1., 2., 0.}, box_geo); + auto const shifts = positions_in_halo({1., 2., 0.}, box_geo, agrid); BOOST_REQUIRE_EQUAL(shifts.size(), 2); BOOST_REQUIRE_EQUAL(shifts[0], reference_origin); BOOST_REQUIRE_EQUAL(shifts[1], reference_shift); @@ -478,7 +481,7 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, } // remove force of the particle from the fluid set_particle_property(pid, &Particle::force, Utils::Vector3d{}); - add_md_force(p_pos, -expected, params.time_step); + add_md_force(lb, p_pos, -expected, params.time_step); } } @@ -504,6 +507,27 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, } } +BOOST_AUTO_TEST_CASE(runtime_exceptions) { + boost::mpi::communicator world; + auto &lb = System::get_system().lb; + // LB prevents changing most of the system state + { + BOOST_CHECK_THROW(lb.on_boxl_change(), std::runtime_error); + BOOST_CHECK_THROW(lb.on_timestep_change(), std::runtime_error); + BOOST_CHECK_THROW(lb.on_temperature_change(), std::runtime_error); + BOOST_CHECK_THROW(lb.on_node_grid_change(), std::runtime_error); + } + // check access out of the LB local domain + { + if (world.rank() != 0) { + BOOST_CHECK(not lb.get_interpolated_density({0., 0., 0.}).has_value()); + BOOST_CHECK(not lb.get_interpolated_velocity({0., 0., 0.}).has_value()); + BOOST_CHECK_THROW(lb.add_force_density({0., 0., 0.}, {0., 0., 0.}), + std::runtime_error); + } + } +} + BOOST_AUTO_TEST_SUITE_END() bool test_lb_domain_mismatch_local() { @@ -512,16 +536,17 @@ bool test_lb_domain_mismatch_local() { auto const node_grid_reversed = Utils::Vector3i{{::node_grid[2], ::node_grid[1], ::node_grid[0]}}; auto const n_ghost_layers = 1u; - auto const params = LBWalberlaParams(0.5, 0.01); + auto const params = std::make_shared(0.5, 0.01); ::node_grid = node_grid_reversed; auto const lattice = std::make_shared( Utils::Vector3i{12, 12, 12}, node_grid_original, n_ghost_layers); auto const ptr = new_lb_walberla(lattice, 1.0, 1.0, false); ptr->set_collision_model(0.0, 0); ::node_grid = node_grid_original; + auto lb_instance = std::make_shared(ptr, params); if (world.rank() == 0) { try { - lb_sanity_checks(*ptr, params, params.get_tau()); + lb_instance->sanity_checks(); } catch (std::runtime_error const &err) { auto const what_ref = std::string("waLBerla and ESPResSo disagree " "about domain decomposition."); @@ -531,35 +556,27 @@ bool test_lb_domain_mismatch_local() { return false; } -BOOST_AUTO_TEST_CASE(exceptions) { +BOOST_AUTO_TEST_CASE(lb_exceptions) { + boost::mpi::communicator world; + auto &lb = System::get_system().lb; + // LB exceptions mechanism { - boost::mpi::communicator world; using std::exception; - // accessing uninitialized pointers is not allowed - BOOST_CHECK_THROW(lb_walberla(), std::runtime_error); - BOOST_CHECK_THROW(lb_walberla_params(), std::runtime_error); // getters and setters - BOOST_CHECK_THROW(LB::get_agrid(), exception); - BOOST_CHECK_THROW(LB::get_tau(), exception); - BOOST_CHECK_THROW(LB::get_kT(), exception); - BOOST_CHECK_THROW(LB::get_pressure_tensor(world), exception); + BOOST_CHECK_THROW(lb.get_agrid(), exception); + BOOST_CHECK_THROW(lb.get_tau(), exception); + BOOST_CHECK_THROW(lb.get_kT(), exception); + BOOST_CHECK_THROW(lb.get_pressure_tensor(), exception); BOOST_CHECK_THROW(LB::get_force_to_be_applied({-10., -10., -10.}), std::runtime_error); // coupling, interpolation, boundaries BOOST_CHECK_THROW(lb_lbcoupling_get_rng_state(), std::runtime_error); BOOST_CHECK_THROW(lb_lbcoupling_set_rng_state(0ul), std::runtime_error); - BOOST_CHECK_THROW(lb_lbinterpolation_get_interpolated_velocity({}), - std::runtime_error); - BOOST_CHECK_THROW(lb_lbinterpolation_add_force_density({}, {}), - std::runtime_error); - BOOST_CHECK_THROW(LB::get_interpolated_velocity({}), exception); - BOOST_CHECK_THROW(LB::get_interpolated_density({}), exception); - BOOST_CHECK_THROW(LB::calc_fluid_momentum(), exception); + BOOST_CHECK_THROW(lb.get_momentum(), exception); } // waLBerla and ESPResSo must agree on domain decomposition { - boost::mpi::communicator world; auto const has_thrown_correct_exception = test_lb_domain_mismatch_local(); auto const n_errors = check_runtime_errors_local(); auto const error_queue = @@ -574,8 +591,36 @@ BOOST_AUTO_TEST_CASE(exceptions) { auto const error_what = error.what().substr(1, what_ref.size()); BOOST_CHECK_EQUAL(error_what, what_ref); } + } else { + BOOST_TEST_REQUIRE(not has_thrown_correct_exception); } } + + // LBNone exceptions mechanism + { + using LB::NoLBActive; + auto const vec = Utils::Vector3d{}; + auto lb_impl = std::make_shared(); + lb.set(lb_impl); + BOOST_CHECK_THROW(lb.get_agrid(), NoLBActive); + BOOST_CHECK_THROW(lb.get_tau(), NoLBActive); + BOOST_CHECK_THROW(lb.get_kT(), NoLBActive); + BOOST_CHECK_THROW(lb.get_pressure_tensor(), NoLBActive); + BOOST_CHECK_THROW(lb.get_momentum(), NoLBActive); + BOOST_CHECK_THROW(lb.sanity_checks(), NoLBActive); + BOOST_CHECK_THROW(lb.veto_time_step(0.), NoLBActive); + BOOST_CHECK_THROW(lb.lebc_sanity_checks(0u, 1u), NoLBActive); + BOOST_CHECK_THROW(lb.propagate(), NoLBActive); + BOOST_CHECK_THROW(lb.on_cell_structure_change(), NoLBActive); + BOOST_CHECK_THROW(lb.on_boxl_change(), NoLBActive); + BOOST_CHECK_THROW(lb.on_node_grid_change(), NoLBActive); + BOOST_CHECK_THROW(lb.on_timestep_change(), NoLBActive); + BOOST_CHECK_THROW(lb.on_temperature_change(), NoLBActive); + BOOST_CHECK_THROW(lb_impl->get_density_at_pos(vec, true), NoLBActive); + BOOST_CHECK_THROW(lb_impl->get_velocity_at_pos(vec, true), NoLBActive); + BOOST_CHECK_THROW(lb_impl->add_force_at_pos(vec, vec), NoLBActive); + lb.reset(); + } } int main(int argc, char **argv) { diff --git a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp b/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp index 9daad17461f..ab54e7b0530 100644 --- a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp +++ b/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp @@ -25,15 +25,22 @@ #include "cells.hpp" #include "errorhandling.hpp" #include "forces.hpp" -#include "grid_based_algorithms/lb_interface.hpp" -#include "grid_based_algorithms/lb_interpolation.hpp" -#include "grid_based_algorithms/lb_particle_coupling.hpp" +#include "grid.hpp" #include "integrate.hpp" +#include "lb/particle_coupling.hpp" +#include "system/System.hpp" -#include +struct DeferredActiveLBChecks { + DeferredActiveLBChecks() { + auto const &lb = System::get_system().lb; + m_value = lb.is_solver_set(); + } + auto operator()() const { return m_value; } + bool m_value; +}; -static bool lb_active_check() { - if (lattice_switch == ActiveLB::NONE) { +static bool lb_active_check(DeferredActiveLBChecks const &check) { + if (not check()) { runtimeErrorMsg() << "LB needs to be active for inertialess tracers."; return false; } @@ -41,8 +48,10 @@ static bool lb_active_check() { } void VirtualSitesInertialessTracers::after_force_calc(double time_step) { - auto const to_lb_units = - (lattice_switch == ActiveLB::NONE) ? 0. : 1. / LB::get_agrid(); + DeferredActiveLBChecks check_lb_solver_set{}; + auto &system = System::get_system(); + auto const agrid = (check_lb_solver_set()) ? system.lb.get_agrid() : 0.; + auto const to_lb_units = (check_lb_solver_set()) ? 1. / agrid : 0.; // Distribute summed-up forces from physical particles to ghosts init_forces_ghosts(cell_structure.ghost_particles()); @@ -57,12 +66,12 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) { for (auto const &p : particle_range) { if (!p.is_virtual()) continue; - if (!lb_active_check()) { + if (!lb_active_check(check_lb_solver_set)) { return; } if (bookkeeping.should_be_coupled(p)) { - for (auto pos : positions_in_halo(p.pos(), box_geo)) { - add_md_force(pos * to_lb_units, p.force(), time_step); + for (auto pos : positions_in_halo(p.pos(), ::box_geo, agrid)) { + add_md_force(system.lb, pos * to_lb_units, p.force(), time_step); } } } @@ -73,18 +82,18 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) { } void VirtualSitesInertialessTracers::after_lb_propagation(double time_step) { - auto const to_md_units = - (lattice_switch == ActiveLB::NONE) ? 0. : LB::get_lattice_speed(); + DeferredActiveLBChecks check_lb_solver_set{}; + auto const &lb = System::get_system().lb; // Advect particles for (auto &p : cell_structure.local_particles()) { if (!p.is_virtual()) continue; - if (!lb_active_check()) { + if (!lb_active_check(check_lb_solver_set)) { return; } - p.v() = lb_lbinterpolation_get_interpolated_velocity(p.pos()) * to_md_units; - for (unsigned int i = 0; i < 3; i++) { + p.v() = lb.get_coupling_interpolated_velocity(p.pos()); + for (unsigned int i = 0u; i < 3u; i++) { if (!p.is_fixed_along(i)) { p.pos()[i] += p.v()[i] * time_step; } diff --git a/src/python/espressomd/actors.py b/src/python/espressomd/actors.py deleted file mode 100644 index 9ab4abcba7a..00000000000 --- a/src/python/espressomd/actors.py +++ /dev/null @@ -1,99 +0,0 @@ -# Copyright (C) 2010-2022 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 . - -from . import highlander -from . import script_interface - - -class Actors: - - """ - Container for actor objects. - """ - - active_actors = [] - - def __del__(self): - self.clear() - - def __getstate__(self): - return self.active_actors - - def __setstate__(self, active_actors): - self.active_actors[:] = [] - for actor in active_actors: - self.active_actors.append(actor) - actor._activate() - - def add(self, actor): - """ - Parameters - ---------- - actor : - Actor to add to this container. - - """ - if actor in Actors.active_actors: - raise highlander.ThereCanOnlyBeOne(actor) - - if isinstance(actor, script_interface.ScriptInterfaceHelper): - actor._activate() - - self.active_actors.append(actor) - - if not isinstance(actor, script_interface.ScriptInterfaceHelper): - actor._activate() - - def remove(self, actor): - """ - Parameters - ---------- - actor : - Actor to remove from this container. - - """ - if actor not in self.active_actors: - raise Exception("Actor is not active") - actor._deactivate() - self.active_actors.remove(actor) - - def clear(self): - """Remove all actors.""" - # The order in which actors are removed matters. Some actors set up - # global bitfields that activate sanity checks on the MD cellsystem, - # and reset these bitfields when removed. Actors need to be removed - # in the reverse order they were inserted to guarantee pre-conditions - # and post-conditions are always met. - while len(self.active_actors): - self.remove(self.active_actors[-1]) - - def __str__(self): - return str(self.active_actors) - - def __getitem__(self, key): - return self.active_actors[key] - - def __len__(self): - return len(self.active_actors) - - def __iter__(self): - for actor in self.active_actors: - yield actor - - def __delitem__(self, idx): - actor = self[idx] - self.remove(actor) diff --git a/src/python/espressomd/electrokinetics.py b/src/python/espressomd/electrokinetics.py index f7585cd6cc7..1fb2a59ae9a 100644 --- a/src/python/espressomd/electrokinetics.py +++ b/src/python/espressomd/electrokinetics.py @@ -50,21 +50,6 @@ class EKNone(ScriptInterfaceHelper): _so_creation_policy = "GLOBAL" -@script_interface_register -class EKContainer(ScriptObjectList): - _so_name = "walberla::EKContainer" - _so_features = ("WALBERLA",) - - def add(self, ekspecies): - self.call_method("add", object=ekspecies) - - def remove(self, ekspecies): - self.call_method("remove", object=ekspecies) - - def clear(self): - self.call_method("clear") - - @script_interface_register class EKSpecies(ScriptInterfaceHelper, espressomd.detail.walberla.LatticeModel): @@ -665,17 +650,33 @@ def __setitem__(self, key, values): class EKReactions(ScriptObjectList): _so_name = "walberla::EKReactions" _so_creation_policy = "GLOBAL" + _so_bind_methods = ("clear",) def add(self, reaction): if not isinstance(reaction, (EKBulkReaction, EKIndexedReaction)): raise TypeError("reaction object is not of correct type.") - self.call_method("add", object=reaction) - return reaction - def remove(self, reaction): self.call_method("remove", object=reaction) - def clear(self): - self.call_method("clear") + +@script_interface_register +class EKContainer(ScriptObjectList): + _so_name = "walberla::EKContainer" + _so_creation_policy = "GLOBAL" + _so_features = ("WALBERLA",) + _so_bind_methods = ("clear",) + + def __init__(self, *args, **kwargs): + if "sip" not in kwargs: + kwargs["reactions"] = EKReactions() + super().__init__(*args, **kwargs) + else: + super().__init__(**kwargs) + + def add(self, ekspecies): + self.call_method("add", object=ekspecies) + + def remove(self, ekspecies): + self.call_method("remove", object=ekspecies) diff --git a/src/python/espressomd/system.py b/src/python/espressomd/system.py index a8bc698a247..e5cb9e2ec0a 100644 --- a/src/python/espressomd/system.py +++ b/src/python/espressomd/system.py @@ -21,7 +21,6 @@ import collections from . import accumulators -from . import actors from . import analyze from . import bond_breakage from . import cell_system @@ -34,7 +33,6 @@ from . import galilei from . import interactions from . import integrate -from . import electrokinetics from . import lees_edwards from . import particle_data from . import thermostat @@ -86,8 +84,6 @@ class System(ScriptInterfaceHelper): collision_detection: :class:`espressomd.collision_detection.CollisionDetection` comfixed: :class:`espressomd.comfixed.ComFixed` constraints: :class:`espressomd.constraints.Constraints` - ekcontainer: :class:`espressomd.electrokinetics.EKContainer` - ekreactions: :class:`espressomd.electrokinetics.EKReactions` cuda_init_handle: :class:`espressomd.cuda_init.CudaInitHandle` galilei: :class:`espressomd.galilei.GalileiTransform` integrator: :class:`espressomd.integrate.IntegratorHandle` @@ -196,7 +192,6 @@ def __init__(self, **kwargs): raise ValueError( f"Property '{arg}' can not be set via argument to System class.") System.__setattr__(self, arg, kwargs.get(arg)) - self.actors = actors.Actors() self.analysis = analyze.Analysis() self.auto_update_accumulators = accumulators.AutoUpdateAccumulators() self.bonded_inter = interactions.BondedInteractions() @@ -214,8 +209,8 @@ def __init__(self, **kwargs): if has_features("DIPOLES"): self.magnetostatics = magnetostatics.Container() if has_features("WALBERLA"): - self.ekcontainer = electrokinetics.EKContainer() - self.ekreactions = electrokinetics.EKReactions() + self._lb = None + self._ekcontainer = None self.galilei = galilei.GalileiTransform() self.lees_edwards = lees_edwards.LeesEdwards() self.non_bonded_inter = interactions.NonBondedInteractions() @@ -233,7 +228,6 @@ def _setup_atexit(self): import atexit def session_shutdown(): - self.actors.clear() self.call_method("session_shutdown") atexit.register(session_shutdown) @@ -265,9 +259,9 @@ def __getstate__(self): checkpointable_properties.append("magnetostatics") if has_features("ELECTROSTATICS"): checkpointable_properties.append("electrostatics") - checkpointable_properties += ["actors", "thermostat"] if has_features("WALBERLA"): - checkpointable_properties += ["ekcontainer", "ekreactions"] + checkpointable_properties += ["_lb", "_ekcontainer"] + checkpointable_properties += ["thermostat"] odict = collections.OrderedDict() odict["_system_handle"] = self.call_method("get_system_handle") @@ -280,6 +274,14 @@ def __setstate__(self, params): self.call_method("set_system_handle", obj=params.pop("_system_handle")) for property_name in params.keys(): System.__setattr__(self, property_name, params[property_name]) + # note: several members can only be instantiated once + if has_features("WALBERLA"): + if self._lb is not None: + lb, self._lb = self._lb, None + self.lb = lb + if self._ekcontainer is not None: + ekcontainer, self._ekcontainer = self._ekcontainer, None + self.ekcontainer = ekcontainer self.call_method("lock_system_creation") @property @@ -363,6 +365,48 @@ def virtual_sites(self, value): assert_features("VIRTUAL_SITES") self._active_virtual_sites_handle.implementation = value + @property + def lb(self): + """ + LB solver. + + """ + assert_features("WALBERLA") + return self._lb + + @lb.setter + def lb(self, lb): + assert_features("WALBERLA") + if lb != self._lb: + if self._lb is not None: + self._lb.call_method("deactivate") + self._lb = None + if lb is not None: + lb.call_method("activate") + self._lb = lb + + @property + def ekcontainer(self): + """ + EK system (diffusion-advection-reaction models). + + Type: :class:`espressomd.electrokinetics.EKContainer` + + """ + assert_features("WALBERLA") + return self._ekcontainer + + @ekcontainer.setter + def ekcontainer(self, ekcontainer): + assert_features("WALBERLA") + if ekcontainer != self._ekcontainer: + if self._ekcontainer is not None: + self._ekcontainer.call_method("deactivate") + self._ekcontainer = None + if ekcontainer is not None: + ekcontainer.call_method("activate") + self._ekcontainer = ekcontainer + def change_volume_and_rescale_particles(self, d_new, dir="xyz"): """Change box size and rescale particle coordinates. diff --git a/src/python/espressomd/thermostat.pxd b/src/python/espressomd/thermostat.pxd index 4b0394d12a1..ccaceaacfc5 100644 --- a/src/python/espressomd/thermostat.pxd +++ b/src/python/espressomd/thermostat.pxd @@ -115,10 +115,7 @@ cdef extern from "stokesian_dynamics/sd_interface.hpp": void set_sd_kT(double kT) except + double get_sd_kT() -cdef extern from "grid_based_algorithms/lb_interface.hpp": - double lb_lbfluid_get_kT "LB::get_kT"() except + - -cdef extern from "grid_based_algorithms/lb_particle_coupling.hpp": +cdef extern from "lb/particle_coupling.hpp": void lb_lbcoupling_set_rng_state(stdint.uint64_t) except + stdint.uint64_t lb_lbcoupling_get_rng_state() except + void lb_lbcoupling_set_gamma(double) except + diff --git a/src/python/espressomd/thermostat.pyx b/src/python/espressomd/thermostat.pyx index d6ce6e505e8..80baef72f8b 100644 --- a/src/python/espressomd/thermostat.pyx +++ b/src/python/espressomd/thermostat.pyx @@ -19,7 +19,6 @@ import functools include "myconfig.pxi" from . cimport utils -from .lb import HydrodynamicInteraction def AssertThermostatType(*allowed_thermostats): @@ -109,11 +108,19 @@ cdef class Thermostat: seed=thmst["seed"]) langevin_set_rng_counter(thmst["counter"]) if thmst["type"] == "LB": + # The LB fluid object is a deep copy of the original, we must + # attach it to be able to set up the particle coupling global + # variable, and then detach it to allow the original LB + # fluid object to be attached by the System class. + # This workaround will be removed when the thermostat class + # gets rewritten as a ScriptInterface class (#4266, #4594) + thmst["LB_fluid"].call_method("activate") self.set_lb( LB_fluid=thmst["LB_fluid"], act_on_virtual=thmst["act_on_virtual"], gamma=thmst["gamma"], seed=thmst["rng_counter_fluid"]) + thmst["LB_fluid"].call_method("deactivate") if thmst["type"] == "NPT_ISO": if NPT: self.set_npt(kT=thmst["kT"], gamma0=thmst["gamma0"], @@ -537,12 +544,13 @@ cdef class Thermostat: Frictional coupling constant for the MD particle coupling. """ + from .lb import HydrodynamicInteraction if not isinstance(LB_fluid, HydrodynamicInteraction): raise ValueError( "The LB thermostat requires a LB / LBGPU instance as a keyword arg.") self._LB_fluid = LB_fluid - if lb_lbfluid_get_kT() > 0.: + if LB_fluid.kT > 0.: if seed is None and lb_lbcoupling_is_seed_required(): raise ValueError( "seed has to be given as keyword arg") diff --git a/src/python/espressomd/visualization.py b/src/python/espressomd/visualization.py index 1379f9fe81a..ac94b603bdf 100644 --- a/src/python/espressomd/visualization.py +++ b/src/python/espressomd/visualization.py @@ -491,9 +491,17 @@ def update_system_info(self): if len(self.system_info['Bonded interactions']) == 0: self.system_info['Bonded interactions'].append('None') - # collect all actors - for a in self.system.actors: - self.system_info['Actors'].append(str(a)) + # collect all solvers + if self.system.electrostatics.solver: + self.system_info['Actors'].append( + self.system.electrostatics.solver) + if self.system.magnetostatics.solver: + self.system_info['Actors'].append( + self.system.magnetostatics.solver) + if self.system.ekcontainer: + self.system_info['Actors'].append(self.system.ekcontainer) + if self.system.lb: + self.system_info['Actors'].append(self.system.lb) if len(self.system_info['Actors']) == 0: self.system_info['Actors'].append('None') @@ -1594,13 +1602,10 @@ def _init_espresso_visualization(self): self.depth = 0 # LOOK FOR LB ACTOR - lb_types = [espressomd.lb.LBFluidWalberla] - for actor in self.system.actors: - if isinstance(actor, tuple(lb_types)): - self.lb_params = actor.get_params() - self.lb = actor - self.lb_is_cpu = True - break + if self.system.lb: + self.lb = self.system.lb + self.lb_params = self.lb.get_params() + self.lb_is_cpu = True if self.specs['LB_draw_velocity_plane']: if self.specs['LB_plane_axis'] == 0: diff --git a/src/script_interface/lees_edwards/LeesEdwards.hpp b/src/script_interface/lees_edwards/LeesEdwards.hpp index 029931309b5..92a030cf4ca 100644 --- a/src/script_interface/lees_edwards/LeesEdwards.hpp +++ b/src/script_interface/lees_edwards/LeesEdwards.hpp @@ -22,9 +22,9 @@ #include "Protocol.hpp" #include "core/grid.hpp" -#include "core/grid_based_algorithms/lb_interface.hpp" #include "core/lees_edwards/LeesEdwardsBC.hpp" #include "core/lees_edwards/lees_edwards.hpp" +#include "core/system/System.hpp" #include "script_interface/ScriptInterface.hpp" #include "script_interface/auto_parameters/AutoParameters.hpp" @@ -47,7 +47,8 @@ class LeesEdwards : public AutoParameters { if (is_none(value)) { context()->parallel_try_catch([]() { auto constexpr invalid_dir = LeesEdwardsBC::invalid_dir; - LB::lebc_sanity_checks(invalid_dir, invalid_dir); + auto const &system = System::get_system(); + system.lb.lebc_sanity_checks(invalid_dir, invalid_dir); }); m_protocol = nullptr; ::box_geo.set_lees_edwards_bc(LeesEdwardsBC{}); @@ -100,7 +101,8 @@ class LeesEdwards : public AutoParameters { throw std::invalid_argument("Parameters 'shear_direction' and " "'shear_plane_normal' must differ"); } - LB::lebc_sanity_checks(shear_direction, shear_plane_normal); + auto const &system = System::get_system(); + system.lb.lebc_sanity_checks(shear_direction, shear_plane_normal); // update box geometry and cell structure ::box_geo.set_lees_edwards_bc( LeesEdwardsBC{0., 0., shear_direction, shear_plane_normal}); diff --git a/src/script_interface/walberla/EKContainer.hpp b/src/script_interface/walberla/EKContainer.hpp index e5f41554d22..5f33f50c28e 100644 --- a/src/script_interface/walberla/EKContainer.hpp +++ b/src/script_interface/walberla/EKContainer.hpp @@ -25,9 +25,15 @@ #include "EKFFT.hpp" #include "EKNone.hpp" +#include "EKReactions.hpp" #include "EKSpecies.hpp" -#include "core/grid_based_algorithms/ek_container.hpp" +#include +#include + +#include "core/ek/EKWalberla.hpp" +#include "core/ek/Solver.hpp" +#include "core/system/System.hpp" #include #include @@ -36,66 +42,103 @@ #include #include #include +#include namespace ScriptInterface::walberla { class EKContainer : public ObjectList { using Base = ObjectList; - boost::variant< + std::variant< #ifdef WALBERLA_FFT std::shared_ptr, #endif std::shared_ptr> - m_poisson_solver{std::shared_ptr()}; + m_poisson_solver; + + std::shared_ptr m_ek_reactions; + std::shared_ptr<::EK::EKWalberla> m_ek_instance; + std::shared_ptr<::EK::EKWalberla::ek_container_type> m_ek_container; + bool m_is_active; void add_in_core(std::shared_ptr const &obj_ptr) override { context()->parallel_try_catch( - [&obj_ptr]() { EK::ek_container.add(obj_ptr->get_ekinstance()); }); + [this, &obj_ptr]() { m_ek_container->add(obj_ptr->get_ekinstance()); }); } void remove_in_core(std::shared_ptr const &obj_ptr) override { - EK::ek_container.remove(obj_ptr->get_ekinstance()); + m_ek_container->remove(obj_ptr->get_ekinstance()); + } + + struct GetPoissonSolverAsVariant { + template + auto operator()(std::shared_ptr const &solver) const { + return (solver) ? Variant{solver} : Variant{none}; + } + }; + + Variant get_solver() const { + return std::visit(GetPoissonSolverAsVariant(), m_poisson_solver); + } + + struct GetPoissonSolverCoreInstance { + template + std::shared_ptr<::walberla::PoissonSolver> + operator()(std::shared_ptr const &solver) const { + return solver->get_instance(); + } + }; + + auto extract_solver(Variant const &v) { + std::optional solver; + auto so_ptr = get_value(v); + if (auto ptr = std::dynamic_pointer_cast(so_ptr)) { + solver = std::move(ptr); + } +#ifdef WALBERLA_FFT + else if (auto ptr = std::dynamic_pointer_cast(so_ptr)) { + solver = std::move(ptr); + } +#endif + assert(solver.has_value()); + return *solver; + } + + void set_solver(Variant const &v) { + m_poisson_solver = extract_solver(v); + auto handle = std::visit(GetPoissonSolverCoreInstance{}, m_poisson_solver); + context()->parallel_try_catch( + [this, &handle]() { m_ek_container->set_poisson_solver(handle); }); } public: EKContainer() : Base::ObjectList() { add_parameters({ - {"tau", - [this](Variant const &v) { - if (is_none(v)) { - if (get_value(do_call_method("size", {})) == 0) { - EK::ek_container.set_tau(0.); - return; - } - context()->parallel_try_catch([]() { - throw std::domain_error( - "Parameter 'tau' is required when container isn't empty"); - }); - } - auto const tau = get_value(v); - context()->parallel_try_catch([tau]() { - if (tau <= 0.) { - throw std::domain_error("Parameter 'tau' must be > 0"); - } - }); - EK::ek_container.set_tau(get_value(v)); - }, - []() { - auto const tau = EK::ek_container.get_tau(); - return (tau == 0.) ? Variant{none} : Variant{tau}; - }}, + {"tau", AutoParameter::read_only, + [this]() { return m_ek_container->get_tau(); }}, {"solver", [this](Variant const &v) { set_solver(v); }, [this]() { return get_solver(); }}, + {"reactions", AutoParameter::read_only, + [this]() { return m_ek_reactions; }}, + {"is_active", AutoParameter::read_only, + [this]() { return m_is_active; }}, }); } void do_construct(VariantMap const ¶ms) override { - if (params.count("solver")) { - set_solver(params.at("solver")); - } - if (params.count("tau")) { - do_set_parameter("tau", params.at("tau")); - } + m_is_active = false; + auto const tau = get_value(params, "tau"); + context()->parallel_try_catch([tau]() { + if (tau <= 0.) { + throw std::domain_error("Parameter 'tau' must be > 0"); + } + }); + m_poisson_solver = extract_solver( + (params.count("solver") == 1) ? params.at("solver") : Variant{none}); + m_ek_container = std::make_shared<::EK::EKWalberla::ek_container_type>( + tau, std::visit(GetPoissonSolverCoreInstance{}, m_poisson_solver)); + m_ek_reactions = get_value(params, "reactions"); + m_ek_instance = std::make_shared<::EK::EKWalberla>( + m_ek_container, m_ek_reactions->get_handle()); // EK species must be added after tau Base::do_construct(params); } @@ -103,57 +146,22 @@ class EKContainer : public ObjectList { protected: Variant do_call_method(std::string const &method, VariantMap const ¶meters) override { - if (method == "is_poisson_solver_set") { - return EK::ek_container.is_poisson_solver_set(); - } - - return Base::do_call_method(method, parameters); - } - -private: - struct GetPoissonSolverVariant : public boost::static_visitor { - template - auto operator()(std::shared_ptr const &solver) const { - return (solver) ? Variant{solver} : Variant{none}; - } - }; - - struct GetPoissonSolverInstance - : public boost::static_visitor< - std::shared_ptr<::walberla::PoissonSolver>> { - template - auto operator()(std::shared_ptr const &solver) const { - return (solver) ? solver->get_instance() - : std::shared_ptr<::walberla::PoissonSolver>(); + if (method == "activate") { + context()->parallel_try_catch([this]() { + ::System::get_system().ek.set<::EK::EKWalberla>(m_ek_instance); + }); + m_is_active = true; + return {}; } - }; - - Variant get_solver() const { - auto const visitor = GetPoissonSolverVariant(); - return boost::apply_visitor(visitor, m_poisson_solver); - } - - void set_solver(Variant const &solver_variant) { - std::optional solver; - if (is_none(solver_variant)) { - solver = std::shared_ptr(); - } else { -#ifdef WALBERLA_FFT - try { - solver = get_value>(solver_variant); - } catch (...) { - } -#endif - if (not solver) { - solver = get_value>(solver_variant); + if (method == "deactivate") { + if (m_is_active) { + ::System::get_system().ek.reset(); + m_is_active = false; } + return {}; } - assert(solver.has_value()); - m_poisson_solver = *solver; - auto const visitor = GetPoissonSolverInstance(); - auto const instance = boost::apply_visitor(visitor, m_poisson_solver); - context()->parallel_try_catch( - [&instance]() { EK::ek_container.set_poisson_solver(instance); }); + + return Base::do_call_method(method, parameters); } }; diff --git a/src/script_interface/walberla/EKFFT.hpp b/src/script_interface/walberla/EKFFT.hpp index d091b001469..c02fc1bae68 100644 --- a/src/script_interface/walberla/EKFFT.hpp +++ b/src/script_interface/walberla/EKFFT.hpp @@ -25,8 +25,8 @@ #ifdef WALBERLA_FFT #include "EKPoissonSolver.hpp" +#include "LatticeWalberla.hpp" -#include #include #include @@ -47,7 +47,7 @@ class EKFFT : public EKPoissonSolver { public: void do_construct(VariantMap const &args) override { m_single_precision = get_value_or(args, "single_precision", false); - m_lattice = get_value>(args, "lattice"); + m_lattice = get_value(args, "lattice"); // unit conversions auto const agrid = get_value(m_lattice->get_parameter("agrid")); diff --git a/src/script_interface/walberla/EKReactions.hpp b/src/script_interface/walberla/EKReactions.hpp index cb15a361109..84104049454 100644 --- a/src/script_interface/walberla/EKReactions.hpp +++ b/src/script_interface/walberla/EKReactions.hpp @@ -23,9 +23,12 @@ #ifdef WALBERLA -#include "EKReaction.hpp" +#include + +#include "core/ek/EKReactions.hpp" +#include "core/ek/EKWalberla.hpp" -#include "core/grid_based_algorithms/ek_reactions.hpp" +#include "EKReaction.hpp" #include #include @@ -35,13 +38,24 @@ namespace ScriptInterface::walberla { class EKReactions : public ObjectList { + std::shared_ptr<::EK::EKWalberla::ek_reactions_type> m_ek_reactions; + void add_in_core(std::shared_ptr const &obj_ptr) override { - EK::ek_reactions.add(obj_ptr->get_instance()); + m_ek_reactions->add(obj_ptr->get_instance()); } void remove_in_core(std::shared_ptr const &obj_ptr) override { - EK::ek_reactions.remove(obj_ptr->get_instance()); + m_ek_reactions->remove(obj_ptr->get_instance()); } + +protected: + void do_construct(VariantMap const ¶ms) override { + m_ek_reactions = std::make_shared<::EK::EKWalberla::ek_reactions_type>(); + } + +public: + auto &get_handle() { return m_ek_reactions; } }; + } // namespace ScriptInterface::walberla #endif // WALBERLA diff --git a/src/script_interface/walberla/LBFluid.cpp b/src/script_interface/walberla/LBFluid.cpp index 39bd39f0385..7c1a5a93a44 100644 --- a/src/script_interface/walberla/LBFluid.cpp +++ b/src/script_interface/walberla/LBFluid.cpp @@ -27,10 +27,11 @@ #include "core/BoxGeometry.hpp" #include "core/event.hpp" #include "core/grid.hpp" -#include "core/grid_based_algorithms/lb_walberla_instance.hpp" #include "core/integrate.hpp" +#include "core/lb/LBWalberla.hpp" #include "core/lees_edwards/lees_edwards.hpp" #include "core/lees_edwards/protocols.hpp" +#include "core/system/System.hpp" #include @@ -68,14 +69,17 @@ std::unordered_map const LBVTKHandle::obs_map = { Variant LBFluid::do_call_method(std::string const &name, VariantMap const ¶ms) { if (name == "activate") { - context()->parallel_try_catch( - [&]() { ::activate_lb_walberla(m_instance, m_lb_params); }); + context()->parallel_try_catch([this]() { + ::System::get_system().lb.set<::LB::LBWalberla>(m_instance, m_lb_params); + }); m_is_active = true; return {}; } if (name == "deactivate") { - ::deactivate_lb_walberla(); - m_is_active = false; + if (m_is_active) { + ::System::get_system().lb.reset(); + m_is_active = false; + } return {}; } if (name == "add_force_at_pos") { @@ -134,7 +138,7 @@ void LBFluid::do_construct(VariantMap const ¶ms) { auto const kT = get_value(params, "kT"); auto const ext_f = get_value(params, "ext_force_density"); auto const single_precision = get_value(params, "single_precision"); - m_lb_params = std::make_shared<::LBWalberlaParams>(agrid, tau); + m_lb_params = std::make_shared<::LB::LBWalberlaParams>(agrid, tau); m_is_active = false; m_seed = get_value(params, "seed"); context()->parallel_try_catch([&]() { diff --git a/src/script_interface/walberla/LBFluid.hpp b/src/script_interface/walberla/LBFluid.hpp index a726108624a..093e7159334 100644 --- a/src/script_interface/walberla/LBFluid.hpp +++ b/src/script_interface/walberla/LBFluid.hpp @@ -27,7 +27,7 @@ #include "LatticeWalberla.hpp" #include "VTKHandle.hpp" -#include "core/grid_based_algorithms/lb_walberla_instance.hpp" +#include "core/lb/LBWalberla.hpp" #include @@ -48,7 +48,7 @@ class LBVTKHandle; class LBFluid : public LatticeModel<::LBWalberlaBase, LBVTKHandle> { using Base = LatticeModel<::LBWalberlaBase, LBVTKHandle>; - std::shared_ptr<::LBWalberlaParams> m_lb_params; + std::shared_ptr<::LB::LBWalberlaParams> m_lb_params; bool m_is_active; int m_seed; double m_conv_dist; diff --git a/src/script_interface/walberla/LatticeWalberla.hpp b/src/script_interface/walberla/LatticeWalberla.hpp index 2a0ba74861e..c20fab5fffd 100644 --- a/src/script_interface/walberla/LatticeWalberla.hpp +++ b/src/script_interface/walberla/LatticeWalberla.hpp @@ -40,6 +40,7 @@ namespace ScriptInterface::walberla { class LatticeWalberla : public AutoParameters { std::shared_ptr<::LatticeWalberla> m_lattice; double m_agrid; + Utils::Vector3d m_box_l; public: LatticeWalberla() { @@ -49,11 +50,13 @@ class LatticeWalberla : public AutoParameters { [this]() { return static_cast(m_lattice->get_ghost_layers()); }}, {"shape", AutoParameter::read_only, [this]() { return m_lattice->get_grid_dimensions(); }}, + {"_box_l", AutoParameter::read_only, [this]() { return m_box_l; }}, }); } void do_construct(VariantMap const &args) override { m_agrid = get_value(args, "agrid"); + m_box_l = get_value_or(args, "_box_l", ::box_geo.length()); auto const n_ghost_layers = get_value(args, "n_ghost_layers"); context()->parallel_try_catch([&]() { @@ -63,9 +66,8 @@ class LatticeWalberla : public AutoParameters { if (n_ghost_layers < 0) { throw std::domain_error("Parameter 'n_ghost_layers' must be >= 0"); } - auto const box_size = ::box_geo.length(); auto const grid_dim = - ::LatticeWalberla::calc_grid_dimensions(box_size, m_agrid); + ::LatticeWalberla::calc_grid_dimensions(m_box_l, m_agrid); m_lattice = std::make_shared<::LatticeWalberla>( grid_dim, node_grid, static_cast(n_ghost_layers)); }); diff --git a/src/walberla_bridge/include/walberla_bridge/electrokinetics/EKContainer.hpp b/src/walberla_bridge/include/walberla_bridge/electrokinetics/EKContainer.hpp index 83a5b861c84..80597203083 100644 --- a/src/walberla_bridge/include/walberla_bridge/electrokinetics/EKContainer.hpp +++ b/src/walberla_bridge/include/walberla_bridge/electrokinetics/EKContainer.hpp @@ -38,9 +38,9 @@ template class EKContainer { using const_iterator = typename container_type::const_iterator; private: - container_type m_ekcontainer; - double m_tau{}; + double m_tau; std::shared_ptr m_poisson_solver; + container_type m_ekcontainer; bool lattice_equal(LatticeWalberla const &lhs, LatticeWalberla const &rhs) const { @@ -48,32 +48,19 @@ template class EKContainer { (lhs.get_grid_dimensions() == rhs.get_grid_dimensions()); } - void sanity_checks(std::shared_ptr const &new_ek_epecies) const { - if (m_tau == 0.) { - throw std::runtime_error("EKContainer parameter 'tau' needs to be set"); - } - if (is_poisson_solver_set()) { - if (not lattice_equal(new_ek_epecies->get_lattice(), - m_poisson_solver->get_lattice())) { - throw std::runtime_error("EKSpecies lattice incompatible with existing " - "Poisson solver lattice"); - } - } - if (not m_ekcontainer.empty()) { - auto const &old_ek_species = m_ekcontainer.front(); - if (not lattice_equal(new_ek_epecies->get_lattice(), - old_ek_species->get_lattice())) { - throw std::runtime_error( - "EKSpecies lattice incompatible with existing EKSpecies lattice"); - } + void sanity_checks(std::shared_ptr const &new_ek_species) const { + if (not lattice_equal(new_ek_species->get_lattice(), + m_poisson_solver->get_lattice())) { + throw std::runtime_error("EKSpecies lattice incompatible with existing " + "Poisson solver lattice"); } } void sanity_checks( - std::shared_ptr const &new_solver) const { + std::shared_ptr const &new_ek_solver) const { if (not m_ekcontainer.empty()) { auto const &old_ek_species = m_ekcontainer.front(); - if (not lattice_equal(new_solver->get_lattice(), + if (not lattice_equal(new_ek_solver->get_lattice(), old_ek_species->get_lattice())) { throw std::runtime_error("Poisson solver lattice incompatible with " "existing EKSpecies lattice"); @@ -82,17 +69,22 @@ template class EKContainer { } public: - void add(std::shared_ptr const &ek_species) { - assert(std::find(m_ekcontainer.begin(), m_ekcontainer.end(), ek_species) == - m_ekcontainer.end()); + EKContainer(double tau, std::shared_ptr solver) + : m_tau{tau}, m_poisson_solver{std::move(solver)}, m_ekcontainer{} {} + + bool contains(std::shared_ptr const &ek_species) const noexcept { + return std::find(m_ekcontainer.begin(), m_ekcontainer.end(), ek_species) != + m_ekcontainer.end(); + } + void add(std::shared_ptr const &ek_species) { + assert(not contains(ek_species)); sanity_checks(ek_species); m_ekcontainer.emplace_back(ek_species); } void remove(std::shared_ptr const &ek_species) { - assert(std::find(m_ekcontainer.begin(), m_ekcontainer.end(), ek_species) != - m_ekcontainer.end()); + assert(contains(ek_species)); m_ekcontainer.erase( std::remove(m_ekcontainer.begin(), m_ekcontainer.end(), ek_species), m_ekcontainer.end()); @@ -106,16 +98,11 @@ template class EKContainer { void set_poisson_solver(std::shared_ptr const &solver) { - if (solver != nullptr) { - sanity_checks(solver); - } + assert(solver != nullptr); + sanity_checks(solver); m_poisson_solver = solver; } - [[nodiscard]] bool is_poisson_solver_set() const noexcept { - return m_poisson_solver != nullptr; - } - [[nodiscard]] double get_tau() const noexcept { return m_tau; } void set_tau(double tau) noexcept { m_tau = tau; } @@ -132,4 +119,8 @@ template class EKContainer { [[nodiscard]] std::size_t get_potential_field_id() const { return m_poisson_solver->get_potential_field_id(); } + + LatticeWalberla const &get_lattice() const noexcept { + return m_poisson_solver->get_lattice(); + } }; diff --git a/src/walberla_bridge/tests/CMakeLists.txt b/src/walberla_bridge/tests/CMakeLists.txt index 581a123f449..4cbf9ce15d7 100644 --- a/src/walberla_bridge/tests/CMakeLists.txt +++ b/src/walberla_bridge/tests/CMakeLists.txt @@ -33,7 +33,11 @@ endfunction() espresso_walberla_unit_test(NAME ResourceManager_test SRC ResourceManager_test.cpp) -espresso_walberla_unit_test(NAME kernels_unit_tests SRC kernels_unit_tests.cpp) +espresso_walberla_unit_test(NAME lb_kernels_unit_tests SRC + lb_kernels_unit_tests.cpp) + +espresso_walberla_unit_test(NAME ek_kernels_unit_tests SRC + ek_kernels_unit_tests.cpp) espresso_walberla_unit_test( NAME LatticeWalberla_unit_tests SRC LatticeWalberla_unit_tests.cpp DEPENDS diff --git a/src/walberla_bridge/tests/ek_kernels_unit_tests.cpp b/src/walberla_bridge/tests/ek_kernels_unit_tests.cpp new file mode 100644 index 00000000000..c209257ba4d --- /dev/null +++ b/src/walberla_bridge/tests/ek_kernels_unit_tests.cpp @@ -0,0 +1,167 @@ +/* + * Copyright (C) 2023 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 . + */ +#define BOOST_TEST_MODULE waLBerla EK kernels +#define BOOST_TEST_DYN_LINK + +#include "config/config.hpp" + +#ifdef WALBERLA + +#include + +#include "../src/electrokinetics/generated_kernels/Dirichlet_double_precision.h" +#include "../src/electrokinetics/generated_kernels/Dirichlet_single_precision.h" +#include "../src/electrokinetics/generated_kernels/FixedFlux_double_precision.h" +#include "../src/electrokinetics/generated_kernels/FixedFlux_single_precision.h" + +bool operator!=( + walberla::pystencils::Dirichlet_single_precision::IndexInfo const &lhs, + walberla::pystencils::Dirichlet_single_precision::IndexInfo const &rhs) { + return not(lhs == rhs); +} + +bool operator!=( + walberla::pystencils::Dirichlet_double_precision::IndexInfo const &lhs, + walberla::pystencils::Dirichlet_double_precision::IndexInfo const &rhs) { + return not(lhs == rhs); +} + +bool operator!=( + walberla::pystencils::Dirichlet_single_precision::IndexVectors const &lhs, + walberla::pystencils::Dirichlet_single_precision::IndexVectors const &rhs) { + return not(lhs == rhs); +} + +bool operator!=( + walberla::pystencils::Dirichlet_double_precision::IndexVectors const &lhs, + walberla::pystencils::Dirichlet_double_precision::IndexVectors const &rhs) { + return not(lhs == rhs); +} + +bool operator!=( + walberla::pystencils::FixedFlux_single_precision::IndexInfo const &lhs, + walberla::pystencils::FixedFlux_single_precision::IndexInfo const &rhs) { + return not(lhs == rhs); +} + +bool operator!=( + walberla::pystencils::FixedFlux_double_precision::IndexInfo const &lhs, + walberla::pystencils::FixedFlux_double_precision::IndexInfo const &rhs) { + return not(lhs == rhs); +} + +bool operator!=( + walberla::pystencils::FixedFlux_single_precision::IndexVectors const &lhs, + walberla::pystencils::FixedFlux_single_precision::IndexVectors const &rhs) { + return not(lhs == rhs); +} + +bool operator!=( + walberla::pystencils::FixedFlux_double_precision::IndexVectors const &lhs, + walberla::pystencils::FixedFlux_double_precision::IndexVectors const &rhs) { + return not(lhs == rhs); +} + +BOOST_AUTO_TEST_CASE(dirichlet_density) { + using Dirichlet_f = walberla::pystencils::Dirichlet_single_precision; + using Dirichlet_d = walberla::pystencils::Dirichlet_double_precision; + + // check IndexInfo + auto dens1_f = Dirichlet_f::IndexInfo(1, 2, 3, 0); + auto dens2_f = Dirichlet_f::IndexInfo(1, 2, 3, 0); + auto dens1_d = Dirichlet_d::IndexInfo(1, 2, 3, 0); + auto dens2_d = Dirichlet_d::IndexInfo(1, 2, 3, 0); + dens1_f.value = dens2_f.value = 2.0f; + dens1_d.value = dens2_d.value = 2.0; + BOOST_TEST((dens1_f == dens2_f)); + BOOST_TEST((dens1_d == dens2_d)); + dens2_f.value += 1.0f; + dens2_d.value += 1.0; + BOOST_TEST((dens1_f != dens2_f)); + BOOST_TEST((dens1_d != dens2_d)); + dens2_f.value += 1.0f; + dens2_d.value += 1.0; + BOOST_TEST((dens1_f != dens2_f)); + BOOST_TEST((dens1_d != dens2_d)); + + // check IndexVector + auto vec1_f = Dirichlet_f::IndexVectors(); + auto vec2_f = Dirichlet_f::IndexVectors(); + auto vec1_d = Dirichlet_d::IndexVectors(); + auto vec2_d = Dirichlet_d::IndexVectors(); + vec1_f.indexVector(Dirichlet_f::IndexVectors::Type::ALL).push_back(dens1_f); + vec2_f.indexVector(Dirichlet_f::IndexVectors::Type::ALL).push_back(dens1_f); + vec1_d.indexVector(Dirichlet_d::IndexVectors::Type::ALL).push_back(dens1_d); + vec2_d.indexVector(Dirichlet_d::IndexVectors::Type::ALL).push_back(dens1_d); + BOOST_TEST((vec1_f == vec2_f)); + BOOST_TEST((vec1_d == vec2_d)); + vec1_f.indexVector(Dirichlet_f::IndexVectors::Type::ALL).push_back(dens1_f); + vec2_f.indexVector(Dirichlet_f::IndexVectors::Type::ALL).push_back(dens2_f); + vec1_d.indexVector(Dirichlet_d::IndexVectors::Type::ALL).push_back(dens1_d); + vec2_d.indexVector(Dirichlet_d::IndexVectors::Type::ALL).push_back(dens2_d); + BOOST_TEST((vec1_f != vec2_f)); + BOOST_TEST((vec1_d != vec2_d)); +} + +BOOST_AUTO_TEST_CASE(dirichlet_flux) { + using FixedFlux_f = walberla::pystencils::FixedFlux_single_precision; + using FixedFlux_d = walberla::pystencils::FixedFlux_double_precision; + + // check IndexInfo + auto flux1_f = FixedFlux_f::IndexInfo(1, 2, 3, 0); + auto flux2_f = FixedFlux_f::IndexInfo(1, 2, 3, 0); + auto flux1_d = FixedFlux_d::IndexInfo(1, 2, 3, 0); + auto flux2_d = FixedFlux_d::IndexInfo(1, 2, 3, 0); + flux1_f.flux_0 = flux2_f.flux_0 = 1.0f; + flux1_f.flux_1 = flux2_f.flux_1 = 2.0f; + flux1_f.flux_2 = flux2_f.flux_2 = 3.0f; + flux1_d.flux_0 = flux2_d.flux_0 = 1.0; + flux1_d.flux_1 = flux2_d.flux_1 = 2.0; + flux1_d.flux_2 = flux2_d.flux_2 = 3.0; + BOOST_TEST((flux1_f == flux2_f)); + BOOST_TEST((flux1_d == flux2_d)); + flux2_f.flux_2 += 1.0f; + flux2_d.flux_2 += 1.0; + BOOST_TEST((flux1_f != flux2_f)); + BOOST_TEST((flux1_d != flux2_d)); + flux2_f.flux_2 += 1.0f; + flux2_d.flux_2 += 1.0; + BOOST_TEST((flux1_f != flux2_f)); + BOOST_TEST((flux1_d != flux2_d)); + + // check IndexVector + auto vec1_f = FixedFlux_f::IndexVectors(); + auto vec2_f = FixedFlux_f::IndexVectors(); + auto vec1_d = FixedFlux_d::IndexVectors(); + auto vec2_d = FixedFlux_d::IndexVectors(); + vec1_f.indexVector(FixedFlux_f::IndexVectors::Type::ALL).push_back(flux1_f); + vec2_f.indexVector(FixedFlux_f::IndexVectors::Type::ALL).push_back(flux1_f); + vec1_d.indexVector(FixedFlux_d::IndexVectors::Type::ALL).push_back(flux1_d); + vec2_d.indexVector(FixedFlux_d::IndexVectors::Type::ALL).push_back(flux1_d); + BOOST_TEST((vec1_f == vec2_f)); + BOOST_TEST((vec1_d == vec2_d)); + vec1_f.indexVector(FixedFlux_f::IndexVectors::Type::ALL).push_back(flux1_f); + vec2_f.indexVector(FixedFlux_f::IndexVectors::Type::ALL).push_back(flux2_f); + vec1_d.indexVector(FixedFlux_d::IndexVectors::Type::ALL).push_back(flux1_d); + vec2_d.indexVector(FixedFlux_d::IndexVectors::Type::ALL).push_back(flux2_d); + BOOST_TEST((vec1_f != vec2_f)); + BOOST_TEST((vec1_d != vec2_d)); +} + +#endif diff --git a/src/walberla_bridge/tests/kernels_unit_tests.cpp b/src/walberla_bridge/tests/lb_kernels_unit_tests.cpp similarity index 92% rename from src/walberla_bridge/tests/kernels_unit_tests.cpp rename to src/walberla_bridge/tests/lb_kernels_unit_tests.cpp index fed6eb4f067..2666843c37f 100644 --- a/src/walberla_bridge/tests/kernels_unit_tests.cpp +++ b/src/walberla_bridge/tests/lb_kernels_unit_tests.cpp @@ -16,7 +16,7 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#define BOOST_TEST_MODULE Walberla kernels +#define BOOST_TEST_MODULE waLBerla LB kernels #define BOOST_TEST_DYN_LINK #include "config/config.hpp" @@ -38,26 +38,26 @@ #include bool operator!=( - const walberla::lbm::Dynamic_UBB_single_precision::IndexInfo &lhs, - const walberla::lbm::Dynamic_UBB_single_precision::IndexInfo &rhs) { + walberla::lbm::Dynamic_UBB_single_precision::IndexInfo const &lhs, + walberla::lbm::Dynamic_UBB_single_precision::IndexInfo const &rhs) { return not(lhs == rhs); } bool operator!=( - const walberla::lbm::Dynamic_UBB_double_precision::IndexInfo &lhs, - const walberla::lbm::Dynamic_UBB_double_precision::IndexInfo &rhs) { + walberla::lbm::Dynamic_UBB_double_precision::IndexInfo const &lhs, + walberla::lbm::Dynamic_UBB_double_precision::IndexInfo const &rhs) { return not(lhs == rhs); } bool operator!=( - const walberla::lbm::Dynamic_UBB_single_precision::IndexVectors &lhs, - const walberla::lbm::Dynamic_UBB_single_precision::IndexVectors &rhs) { + walberla::lbm::Dynamic_UBB_single_precision::IndexVectors const &lhs, + walberla::lbm::Dynamic_UBB_single_precision::IndexVectors const &rhs) { return not(lhs == rhs); } bool operator!=( - const walberla::lbm::Dynamic_UBB_double_precision::IndexVectors &lhs, - const walberla::lbm::Dynamic_UBB_double_precision::IndexVectors &rhs) { + walberla::lbm::Dynamic_UBB_double_precision::IndexVectors const &lhs, + walberla::lbm::Dynamic_UBB_double_precision::IndexVectors const &rhs) { return not(lhs == rhs); } diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index f504423711d..ce792497bc6 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -329,7 +329,6 @@ python_test(FILE random_pairs.py MAX_NUM_PROC 4) python_test(FILE lb_electrohydrodynamics.py MAX_NUM_PROC 4 GPU_SLOTS 1) python_test(FILE cluster_analysis.py MAX_NUM_PROC 4) python_test(FILE pair_criteria.py MAX_NUM_PROC 4) -python_test(FILE actor.py MAX_NUM_PROC 1) python_test(FILE drude.py MAX_NUM_PROC 2) python_test(FILE thermostats_anisotropic.py MAX_NUM_PROC 4) python_test(FILE thermalized_bond.py MAX_NUM_PROC 4) diff --git a/testsuite/python/actor.py b/testsuite/python/actor.py deleted file mode 100644 index 9687537a328..00000000000 --- a/testsuite/python/actor.py +++ /dev/null @@ -1,258 +0,0 @@ -# -# Copyright (C) 2018-2022 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 . -# - -""" -Testmodule for the actor base class. - -""" - -import unittest as ut -import espressomd.actors -import espressomd.highlander -import espressomd.utils as utils - - -class BaseActor: - - """ - Abstract base class for interactions affecting particles in the system, - such as LB fluids. Derived classes must implement the interface to the - relevant core objects and global variables. - """ - - # Keys in active_list have to match the method name. - active_list = dict(HydrodynamicInteraction=False) - - def __init__(self, **kwargs): - self._isactive = False - utils.check_valid_keys(self.valid_keys(), kwargs.keys()) - utils.check_required_keys(self.required_keys(), kwargs.keys()) - self._params = self.default_params() - self._params.update(kwargs) - - def _activate(self): - inter = self._get_interaction_type() - if inter in BaseActor.active_list: - if BaseActor.active_list[inter]: - raise espressomd.highlander.ThereCanOnlyBeOne( - self.__class__.__bases__[0]) - BaseActor.active_list[inter] = True - - self.validate_params() - self._activate_method() - utils.handle_errors("Activation of an actor") - self._isactive = True - - def _deactivate(self): - self._deactivate_method() - utils.handle_errors("Deactivation of an actor") - self._isactive = False - inter = self._get_interaction_type() - if inter in BaseActor.active_list: - if not BaseActor.active_list[inter]: - raise Exception( - f"Class not registered in Actor.active_list: {self.__class__.__bases__[0].__name__}") - BaseActor.active_list[inter] = False - - def get_params(self): - """Get interaction parameters""" - # If this instance refers to an actual interaction defined in the es - # core, load current parameters from there - if self.is_active(): - update = self._get_params_from_es_core() - self._params.update(update) - return self._params - - def set_params(self, **p): - """Update the given parameters.""" - # Check if keys are valid - utils.check_valid_keys(self.valid_keys(), p.keys()) - - # When an interaction is newly activated, all required keys must be - # given - if not self.is_active(): - utils.check_required_keys(self.required_keys(), p.keys()) - - self._params.update(p) - # validate updated parameters - self.validate_params() - # Put in values given by the user - if self.is_active(): - self._set_params_in_es_core() - - def is_active(self): - return self._isactive - - -class TestActor(BaseActor): - - def __init__(self, *args, **kwargs): - self._core_args = None - self._activated = False - self._deactivated = False - self._validated = False - - super().__init__(*args, **kwargs) - - def _get_params_from_es_core(self): - return self._core_args - - def _set_params_in_es_core(self): - self._core_args = self._params - - def valid_keys(self): - return {"a", "b", "c"} - - def required_keys(self): - return {"a", "c"} - - def default_params(self): - return {"a": False, "b": False, "c": False} - - def _activate_method(self): - self._activated = True - - def _deactivate_method(self): - self._deactivated = True - - def validate_params(self): - self._validated = True - - def _get_interaction_type(self): - return None - - -class TestHydrodynamicActor(TestActor): - - def _get_interaction_type(self): - return "HydrodynamicInteraction" - - -class ActorTest(ut.TestCase): - - def test_ctor(self): - a = TestActor(a=False, c=False) - self.assertFalse(a.is_active()) - self.assertEqual(a.get_params(), a.default_params()) - - def test_params_non_active(self): - a = TestActor(a=True, c=True) - - a.set_params(a=False, b=True, c=False) - params = a.get_params() - self.assertEqual(params["a"], False) - self.assertEqual(params["b"], True) - self.assertEqual(params["c"], False) - self.assertEqual(a._core_args, None) - - def test_params_active(self): - a = TestActor(a=True, c=True) - a._activate() - - a.set_params(a=False, b=True, c=False) - params = a.get_params() - self.assertEqual(params["a"], False) - self.assertEqual(params["b"], True) - self.assertEqual(params["c"], False) - self.assertEqual(a._core_args, params) - - def test_activation(self): - a = TestActor(a=True, c=True) - a._activate() - self.assertTrue(a.is_active()) - - def test_deactivation(self): - a = TestActor(a=True, c=True) - a._activate() - self.assertTrue(a.is_active()) - a._deactivate() - self.assertFalse(a.is_active()) - - params = a.get_params() - self.assertEqual(params["a"], True) - self.assertEqual(params["b"], False) - self.assertEqual(params["c"], True) - - def test_exception(self): - error_msg_valid = (r"Only the following keys can be given as keyword arguments: " - r"\['a', 'b', 'c'\], got \['a', 'c', 'd'\] \(unknown \['d'\]\)") - error_msg_required = (r"The following keys have to be given as keyword arguments: " - r"\['a', 'c'\], got \['a'\] \(missing \['c'\]\)") - with self.assertRaisesRegex(ValueError, error_msg_valid): - TestActor(a=True, c=True, d=True) - with self.assertRaisesRegex(ValueError, error_msg_required): - TestActor(a=True) - valid_actor = TestActor(a=True, c=True) - with self.assertRaisesRegex(ValueError, error_msg_valid): - valid_actor.set_params(a=True, c=True, d=True) - with self.assertRaisesRegex(ValueError, error_msg_required): - valid_actor.set_params(a=True) - - -class ActorsTest(ut.TestCase): - - actors = espressomd.actors.Actors() - - def tearDown(self): - self.actors.clear() - - def test_clear(self): - # clearing the list of actors removes all of them - for actors_size in range(10): - for _ in range(actors_size): - actor = TestActor(a=False, c=False) - self.actors.add(actor) - self.assertEqual(len(self.actors), actors_size) - self.actors.clear() - self.assertEqual(len(self.actors), 0) - - def test_deactivation(self): - actor = TestActor(a=False, c=False) - self.assertFalse(actor.is_active()) - # adding an actor activates it - self.actors.add(actor) - self.assertTrue(actor.is_active()) - # removing an actor deactivates it - self.actors.clear() - self.assertFalse(actor.is_active()) - # re-adding an actor re-activates it - self.actors.add(actor) - self.assertTrue(actor.is_active()) - # removing an actor deactivates it - del self.actors[0] - self.assertFalse(actor.is_active()) - - def test_unique(self): - # an actor can only be added once - actor = TestHydrodynamicActor(a=False, c=False) - self.actors.add(actor) - with self.assertRaises(espressomd.highlander.ThereCanOnlyBeOne): - self.actors.add(actor) - with self.assertRaises(espressomd.highlander.ThereCanOnlyBeOne): - actor._activate() - # an actor can only be removed once - self.actors.remove(actor) - with self.assertRaisesRegex(Exception, "Actor is not active"): - self.actors.remove(actor) - with self.assertRaisesRegex(Exception, "Class not registered.*: TestActor"): - actor._deactivate() - - -if __name__ == "__main__": - ut.main() diff --git a/testsuite/python/array_properties.py b/testsuite/python/array_properties.py index 8a205bd9e6c..85a036bd3d9 100644 --- a/testsuite/python/array_properties.py +++ b/testsuite/python/array_properties.py @@ -108,7 +108,8 @@ def setUp(self): self.system.box_l = [12.0, 12.0, 12.0] def tearDown(self): - self.system.actors.clear() + if espressomd.has_features("WALBERLA"): + self.system.lb = None def assert_copy_is_writable(self, array): cpy = np.copy(array) @@ -191,7 +192,7 @@ def test_rot_aniso(self): def test_lb(self): lbf = espressomd.lb.LBFluidWalberla( agrid=0.5, density=1., kinematic_viscosity=1., tau=0.01) - self.system.actors.add(lbf) + self.system.lb = lbf self.assert_operator_usage_raises(lbf[0, 0, 0].velocity) self.assert_operator_usage_raises(lbf[0, 0, 0].pressure_tensor) diff --git a/testsuite/python/ek_boundary.py b/testsuite/python/ek_boundary.py index 1c5da29bceb..87780779d1f 100644 --- a/testsuite/python/ek_boundary.py +++ b/testsuite/python/ek_boundary.py @@ -30,6 +30,7 @@ class EKBoundariesBase: system = espressomd.System(box_l=[10.0, 5.0, 5.0]) system.cell_system.skin = 0.1 + system.time_step = 1.0 ek_species_params = {"kT": 1.5, "density": 0.85, "valency": 0.0, @@ -43,9 +44,12 @@ class EKBoundariesBase: def setUp(self): self.lattice = self.ek_lattice_class(agrid=0.5, n_ghost_layers=1) + ek_solver = espressomd.electrokinetics.EKNone(lattice=self.lattice) + self.system.ekcontainer = espressomd.electrokinetics.EKContainer( + tau=self.ek_species_params["tau"], solver=ek_solver) def tearDown(self): - self.system.ekcontainer.clear() + self.system.ekcontainer = None def make_default_ek_species(self): return self.ek_species_class( diff --git a/testsuite/python/ek_bulk_reactions.py b/testsuite/python/ek_bulk_reactions.py index be6c539458f..39305a5b942 100644 --- a/testsuite/python/ek_bulk_reactions.py +++ b/testsuite/python/ek_bulk_reactions.py @@ -37,9 +37,8 @@ class EKReaction(ut.TestCase): system.time_step = TAU system.cell_system.skin = 0.4 - def tearDown(self) -> None: - self.system.ekcontainer.clear() - self.system.ekreactions.clear() + def tearDown(self): + self.system.ekcontainer = None def analytic_density_base( self, time: float, coeffs, rate_constant: float, init_density: float) -> float: @@ -71,8 +70,8 @@ def detail_test_reaction(self, single_precision: bool): n_ghost_layers=1, agrid=self.AGRID) eksolver = espressomd.electrokinetics.EKNone(lattice=lattice) - - self.system.ekcontainer.tau = self.TAU + self.system.ekcontainer = espressomd.electrokinetics.EKContainer( + tau=self.TAU, solver=eksolver) reaction_rate: float = 1e-5 @@ -105,12 +104,10 @@ def detail_test_reaction(self, single_precision: bool): stoech_coeff=product_coeff, order=0.0)) - self.system.ekcontainer.solver = eksolver - reaction = espressomd.electrokinetics.EKBulkReaction( reactants=reactants, coefficient=reaction_rate, lattice=lattice, tau=self.TAU) - self.system.ekreactions.add(reaction) + self.system.ekcontainer.reactions.add(reaction) self.system.integrator.run(self.TIME) diff --git a/testsuite/python/ek_diffusion.py b/testsuite/python/ek_diffusion.py index 0f979ae551b..f60e23a8736 100644 --- a/testsuite/python/ek_diffusion.py +++ b/testsuite/python/ek_diffusion.py @@ -38,8 +38,8 @@ class EKDiffusion(ut.TestCase): system.time_step = TAU system.cell_system.skin = 0.4 - def tearDown(self) -> None: - self.system.ekcontainer.clear() + def tearDown(self): + self.system.ekcontainer = None def analytical_density(self, pos: np.ndarray, time: int, D: float): return (4 * np.pi * D * time)**(-3 / 2) * \ @@ -68,8 +68,8 @@ def detail_test_diffusion(self, single_precision: bool): eksolver = espressomd.electrokinetics.EKNone(lattice=lattice) - self.system.ekcontainer.tau = self.TAU - self.system.ekcontainer.solver = eksolver + self.system.ekcontainer = espressomd.electrokinetics.EKContainer( + tau=self.TAU, solver=eksolver) self.system.ekcontainer.add(ekspecies) center = np.asarray(lattice.shape // 2, dtype=int) diff --git a/testsuite/python/ek_eof.py b/testsuite/python/ek_eof.py index c66f4450524..500b9042baa 100644 --- a/testsuite/python/ek_eof.py +++ b/testsuite/python/ek_eof.py @@ -42,8 +42,8 @@ class EKEOF: system.cell_system.skin = 0.4 def tearDown(self): - self.system.actors.clear() - self.system.ekcontainer.clear() + self.system.lb = None + self.system.ekcontainer = None def test_eof(self): """ @@ -78,15 +78,15 @@ def test_eof(self): eksolver = self.ek_solver_class( lattice=lattice, permittivity=eps0 * epsR, **self.ek_params) - self.system.ekcontainer.tau = self.TAU - self.system.ekcontainer.solver = eksolver + self.system.ekcontainer = espressomd.electrokinetics.EKContainer( + tau=self.TAU, solver=eksolver) self.system.ekcontainer.add(ekspecies) self.system.ekcontainer.add(ekwallcharge) lb_fluid = espressomd.lb.LBFluidWalberla( lattice=lattice, density=1.0, kinematic_viscosity=visc, tau=self.TAU, **self.ek_params) - self.system.actors.add(lb_fluid) + self.system.lb = lb_fluid wall_bot = espressomd.shapes.Wall(normal=[1, 0, 0], dist=offset) wall_top = espressomd.shapes.Wall( diff --git a/testsuite/python/ek_fixeddensity.py b/testsuite/python/ek_fixeddensity.py index c1b3f29d13e..c5c6252a343 100644 --- a/testsuite/python/ek_fixeddensity.py +++ b/testsuite/python/ek_fixeddensity.py @@ -43,7 +43,7 @@ class EKFixedDensity(ut.TestCase): system.cell_system.skin = 0.4 def tearDown(self) -> None: - self.system.ekcontainer.clear() + self.system.ekcontainer = None def test_constant_density_bc_single(self): self.detail_test_constant_density_bc(single_precision=True) @@ -66,8 +66,8 @@ def detail_test_constant_density_bc(self, single_precision: bool): eksolver = espressomd.electrokinetics.EKNone(lattice=lattice) - self.system.ekcontainer.tau = self.TAU - self.system.ekcontainer.solver = eksolver + self.system.ekcontainer = espressomd.electrokinetics.EKContainer( + tau=self.TAU, solver=eksolver) self.system.ekcontainer.add(ekspecies) # left and right no flux diff --git a/testsuite/python/ek_fixedflux.py b/testsuite/python/ek_fixedflux.py index 8275d0b9d91..14ced8aba46 100644 --- a/testsuite/python/ek_fixedflux.py +++ b/testsuite/python/ek_fixedflux.py @@ -38,11 +38,11 @@ class EKFixedFlux(ut.TestCase): INFLOW_FLUX = 0.1 system = espressomd.System(box_l=[BOX_L, BOX_L, BOX_L]) - system.time_step = 1.0 + system.time_step = TAU system.cell_system.skin = 0.4 - def tearDown(self) -> None: - self.system.ekcontainer.clear() + def tearDown(self): + self.system.ekcontainer = None def test_inflow_single(self): self.detail_test_inflow(single_precision=True) @@ -67,8 +67,8 @@ def detail_test_inflow(self, single_precision: bool): eksolver = espressomd.electrokinetics.EKNone(lattice=lattice) - self.system.ekcontainer.tau = self.TAU - self.system.ekcontainer.solver = eksolver + self.system.ekcontainer = espressomd.electrokinetics.EKContainer( + tau=self.TAU, solver=eksolver) self.system.ekcontainer.add(ekspecies) ekspecies[1:-1, 1:-1, 1:-1].density = self.DENSITY diff --git a/testsuite/python/ek_indexed_reactions.py b/testsuite/python/ek_indexed_reactions.py index 01d478665c8..952b0017c89 100644 --- a/testsuite/python/ek_indexed_reactions.py +++ b/testsuite/python/ek_indexed_reactions.py @@ -42,9 +42,8 @@ class EKReaction(ut.TestCase): system.time_step = TAU system.cell_system.skin = 0.4 - def tearDown(self) -> None: - self.system.ekcontainer.clear() - self.system.ekreactions.clear() + def tearDown(self): + self.system.ekcontainer = None def analytic_density_profiles( self, width, reaction_rates, diffusion_coefficients, initial_densities, agrid): @@ -76,9 +75,8 @@ def detail_test_reaction(self, single_precision: bool): eksolver = espressomd.electrokinetics.EKNone(lattice=lattice) - self.system.ekcontainer.tau = self.TAU - - self.system.ekcontainer.solver = eksolver + self.system.ekcontainer = espressomd.electrokinetics.EKContainer( + tau=self.TAU, solver=eksolver) species_A = espressomd.electrokinetics.EKSpecies( lattice=lattice, density=self.INITIAL_DENSITIES[0], @@ -126,8 +124,8 @@ def detail_test_reaction(self, single_precision: bool): lattice=lattice, tau=self.TAU) reaction_right[-2, :, :] = True - self.system.ekreactions.add(reaction_left) - self.system.ekreactions.add(reaction_right) + self.system.ekcontainer.reactions.add(reaction_left) + self.system.ekcontainer.reactions.add(reaction_right) wall_left = espressomd.shapes.Wall( normal=[1, 0, 0], dist=self.PADDING * self.AGRID) diff --git a/testsuite/python/ek_interface.py b/testsuite/python/ek_interface.py index 0ce7352927f..9e1e986b114 100644 --- a/testsuite/python/ek_interface.py +++ b/testsuite/python/ek_interface.py @@ -51,36 +51,24 @@ class EKTest: system.cell_system.skin = 1.0 def setUp(self): + self.system.box_l = 3 * [6.0] self.lattice = self.ek_lattice_class( n_ghost_layers=1, agrid=self.params["agrid"]) ek_solver = espressomd.electrokinetics.EKNone(lattice=self.lattice) - self.system.ekcontainer.solver = ek_solver - self.system.ekcontainer.tau = self.system.time_step + self.system.ekcontainer = espressomd.electrokinetics.EKContainer( + tau=self.system.time_step, solver=ek_solver) def tearDown(self): - self.system.actors.clear() - self.system.ekcontainer.clear() + self.system.lb = None + self.system.ekcontainer = None self.system.part.clear() self.system.thermostat.turn_off() self.system.time_step = self.params["tau"] - self.system.ekcontainer.solver = None def make_default_ek_species(self): return self.ek_species_class( lattice=self.lattice, **self.ek_params, **self.ek_species_params) - def test_ek_container_poisson_solver(self): - ekcontainer = self.system.ekcontainer - ekcontainer.solver = None - self.assertFalse(ekcontainer.call_method("is_poisson_solver_set")) - self.assertIsNone(ekcontainer.solver) - ek_solver = espressomd.electrokinetics.EKNone(lattice=self.lattice) - self.assertFalse(ekcontainer.call_method("is_poisson_solver_set")) - ekcontainer.solver = ek_solver - self.assertTrue(ekcontainer.call_method("is_poisson_solver_set")) - self.assertIsInstance(self.system.ekcontainer.solver, - espressomd.electrokinetics.EKNone) - def test_ek_species(self): # inactive species ek_species = self.make_default_ek_species() @@ -190,8 +178,6 @@ def test_ek_fft_solver(self): self.assertAlmostEqual(ek_solver.permittivity, 0.05, delta=self.atol) self.system.ekcontainer.solver = ek_solver - self.assertTrue( - self.system.ekcontainer.call_method("is_poisson_solver_set")) self.assertIsInstance(self.system.ekcontainer.solver, espressomd.electrokinetics.EKFFT) @@ -204,8 +190,6 @@ def test_ek_none_solver(self): self.ek_params["single_precision"]) self.system.ekcontainer.solver = ek_solver - self.assertTrue( - self.system.ekcontainer.call_method("is_poisson_solver_set")) self.assertIsInstance(self.system.ekcontainer.solver, espressomd.electrokinetics.EKNone) @@ -215,26 +199,38 @@ def test_ek_solver_exceptions(self): self.system.ekcontainer.add(ek_species) incompatible_lattice = self.ek_lattice_class( n_ghost_layers=2, agrid=self.params["agrid"]) + incompatible_ek_solver = espressomd.electrokinetics.EKNone( + lattice=incompatible_lattice) incompatible_ek_species = self.ek_species_class( - lattice=incompatible_lattice, **self.ek_params, **self.ek_species_params) + lattice=incompatible_lattice, **self.ek_params, + **self.ek_species_params) with self.assertRaisesRegex(RuntimeError, "EKSpecies lattice incompatible with existing Poisson solver lattice"): self.system.ekcontainer.add(incompatible_ek_species) - with self.assertRaisesRegex(RuntimeError, "EKSpecies lattice incompatible with existing EKSpecies lattice"): - self.system.ekcontainer.solver = None - self.system.ekcontainer.add(incompatible_ek_species) - with self.assertRaisesRegex(ValueError, "Parameter 'tau' is required when container isn't empty"): - self.system.ekcontainer.tau = None with self.assertRaisesRegex(RuntimeError, "Poisson solver lattice incompatible with existing EKSpecies lattice"): self.system.ekcontainer.clear() + self.system.ekcontainer.solver = incompatible_ek_solver self.system.ekcontainer.add(incompatible_ek_species) self.system.ekcontainer.solver = ek_solver - with self.assertRaisesRegex(ValueError, "Parameter 'tau' must be > 0"): - self.system.ekcontainer.tau = 0. - self.assertAlmostEqual( - self.system.ekcontainer.tau, self.system.time_step, delta=1e-7) - self.system.ekcontainer.clear() - self.system.ekcontainer.tau = None - self.assertIsNone(self.system.ekcontainer.tau) + + def test_parameter_change_exceptions(self): + ek_solver = self.system.ekcontainer.solver + ek_species = self.make_default_ek_species() + self.system.ekcontainer.add(ek_species) + self.system.ekcontainer.solver = ek_solver + with self.assertRaisesRegex(Exception, "Temperature change not supported by EK"): + self.system.thermostat.turn_off() + with self.assertRaisesRegex(Exception, "Time step change not supported by EK"): + self.system.time_step /= 2. + if espressomd.has_features("ELECTROSTATICS"): + self.system.electrostatics.solver = espressomd.electrostatics.DH( + prefactor=1., kappa=1., r_cut=1.) # should not fail + self.system.electrostatics.clear() + with self.assertRaisesRegex(RuntimeError, "MD cell geometry change not supported by EK"): + self.system.box_l = [1., 2., 3.] + np.testing.assert_allclose( + np.copy(self.system.box_l), [1., 2., 3.], atol=1e-7) + with self.assertRaisesRegex(RuntimeError, "MPI topology change not supported by EK"): + self.system.cell_system.node_grid = self.system.cell_system.node_grid def test_ek_reactants(self): ek_species = self.make_default_ek_species() @@ -312,13 +308,6 @@ def test_runtime_exceptions(self): print("\nTesting EK runtime error messages:", file=sys.stderr) sys.stderr.flush() - # check exceptions without EK Poisson solver - with self.assertRaisesRegex(Exception, "EK requires a Poisson solver"): - self.system.ekcontainer.solver = None - self.system.integrator.run(1) - ek_solver = espressomd.electrokinetics.EKNone(lattice=self.lattice) - self.system.ekcontainer.solver = ek_solver - # check exceptions without LB force field with self.assertRaisesRegex(Exception, "friction coupling enabled but no force field accessible"): ek_species.friction_coupling = True @@ -342,7 +331,7 @@ def test_runtime_exceptions(self): lb = self.lb_fluid_class( lattice=self.lattice, density=0.5, kinematic_viscosity=3., tau=2. * self.params["tau"], **self.lb_params) - self.system.actors.add(lb) + self.system.lb = lb self.system.integrator.run(1) print("End of EK runtime error messages", file=sys.stderr) diff --git a/testsuite/python/ek_noflux.py b/testsuite/python/ek_noflux.py index 35ab7784b9c..f54048fc2ea 100644 --- a/testsuite/python/ek_noflux.py +++ b/testsuite/python/ek_noflux.py @@ -39,8 +39,8 @@ class EKNoFlux(ut.TestCase): system.time_step = 1.0 system.cell_system.skin = 0.4 - def tearDown(self) -> None: - self.system.ekcontainer.clear() + def tearDown(self): + self.system.ekcontainer = None def test_noflux_single(self): self.detail_test_noflux(single_precision=True) @@ -65,8 +65,8 @@ def detail_test_noflux(self, single_precision: bool): eksolver = espressomd.electrokinetics.EKNone(lattice=lattice) - self.system.ekcontainer.tau = 1.0 - self.system.ekcontainer.solver = eksolver + self.system.ekcontainer = espressomd.electrokinetics.EKContainer( + tau=1., solver=eksolver) self.system.ekcontainer.add(ekspecies) center = np.asarray(self.system.box_l / 2, dtype=int) diff --git a/testsuite/python/engine_lb.py b/testsuite/python/engine_lb.py index 7112e4c87d0..32d948b209a 100644 --- a/testsuite/python/engine_lb.py +++ b/testsuite/python/engine_lb.py @@ -116,12 +116,12 @@ def setUp(self): ) self.set_cellsystem() self.lbf = self.lb_class(**self.LB_params, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf self.system.thermostat.set_lb(LB_fluid=self.lbf, gamma=self.gamma) def tearDown(self): self.system.part.clear() - self.system.actors.clear() + self.system.lb = None self.system.thermostat.turn_off() def test_dipole_particle_addition(self): diff --git a/testsuite/python/ibm.py b/testsuite/python/ibm.py index 20ddb7a0a0f..ed4bfd801f7 100644 --- a/testsuite/python/ibm.py +++ b/testsuite/python/ibm.py @@ -35,7 +35,6 @@ class IBM(ut.TestCase): def tearDown(self): self.system.part.clear() - self.system.actors.clear() self.system.thermostat.turn_off() def compute_dihedral_angle(self, pos0, pos1, pos2, pos3): diff --git a/testsuite/python/lattice_vtk.py b/testsuite/python/lattice_vtk.py index 0e20ddaa6e6..90c43f788d7 100644 --- a/testsuite/python/lattice_vtk.py +++ b/testsuite/python/lattice_vtk.py @@ -117,11 +117,11 @@ def make_actor(self): def add_actor(self): self.lbf = self.make_actor() - self.system.actors.add(self.lbf) + self.system.lb = self.lbf return self.lbf def clear_actors(self): - self.system.actors.clear() + self.system.lb = None @utx.skipIfMissingModules("espressomd.io.vtk") def test_vtk(self): @@ -240,13 +240,13 @@ def make_actor(self): def add_actor(self): self.solver = self.ek_solver(lattice=self.lattice) self.species = self.make_actor() - self.system.ekcontainer.tau = 0.1 - self.system.ekcontainer.solver = self.solver + self.system.ekcontainer = espressomd.electrokinetics.EKContainer( + tau=0.1, solver=self.solver) self.system.ekcontainer.add(self.species) return self.species def clear_actors(self): - self.system.ekcontainer.clear() + self.system.ekcontainer = None @utx.skipIfMissingModules("espressomd.io.vtk") def test_vtk(self): diff --git a/testsuite/python/lb.py b/testsuite/python/lb.py index ef8dd7b3b5c..592d7f434ee 100644 --- a/testsuite/python/lb.py +++ b/testsuite/python/lb.py @@ -54,8 +54,11 @@ class LBTest: system.cell_system.skin = 1.0 interpolation = False + def setUp(self): + self.system.box_l = 3 * [6.0] + def tearDown(self): - self.system.actors.clear() + self.system.lb = None self.system.part.clear() self.system.thermostat.turn_off() self.system.time_step = self.params['tau'] @@ -68,17 +71,17 @@ def test_properties(self): # activated actor lbf = self.lb_class(kT=1.0, seed=42, **self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf self.system.thermostat.set_lb(LB_fluid=lbf, seed=1) self.assertTrue(lbf.is_active) self.check_properties(lbf) - self.system.actors.remove(lbf) + self.system.lb = None # deactivated actor lbf = self.lb_class(kT=1.0, seed=42, **self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf self.system.thermostat.set_lb(LB_fluid=lbf, seed=1) - self.system.actors.remove(lbf) + self.system.lb = None self.assertFalse(lbf.is_active) self.check_properties(lbf) @@ -206,14 +209,10 @@ def make_kwargs(**kwargs): self.lb_class(**make_kwargs(kT=-1., seed=42)) with self.assertRaisesRegex(ValueError, "Parameter 'seed' must be >= 0"): self.lb_class(**make_kwargs(kT=0., seed=-42)) - with self.assertRaisesRegex(RuntimeError, "Cannot add a second LB instance"): - lbf = self.lb_class(**make_kwargs()) - self.system.actors.add(lbf) - lbf.call_method("activate") def test_node_exceptions(self): lbf = self.lb_class(**self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf lb_node = lbf[0, 0, 0] # check exceptions from LB node with self.assertRaisesRegex(RuntimeError, "Property 'boundary_force' is read-only"): @@ -236,7 +235,7 @@ def test_node_exceptions(self): def test_slice_exceptions(self): lbf = self.lb_class(**self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf lb_slice = lbf[:, :, :] # check exceptions from LB slice with self.assertRaisesRegex(RuntimeError, "Property 'boundary_force' is read-only"): @@ -268,7 +267,7 @@ def test_slice_exceptions(self): def test_lb_slice_set_get(self): lbf = self.lb_class(**self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf ref_density = 1. + np.arange(np.prod(lbf.shape)).reshape(lbf.shape) lbf[:, :, :].density = ref_density densities = np.copy(lbf[:, :, :].density) @@ -325,7 +324,7 @@ def test_pressure_tensor_observable(self): lbf = self.lb_class(kT=1., seed=1, ext_force_density=[0, 0, 0], **self.params, **self.lb_params) - system.actors.add(lbf) + system.lb = lbf system.thermostat.set_lb(LB_fluid=lbf, seed=1) system.integrator.run(10) pressure_tensor = np.copy( @@ -345,14 +344,14 @@ def test_pressure_tensor_observable(self): self.assertIsInstance( lbf.pressure_tensor, espressomd.utils.array_locked) - system.actors.remove(lbf) + system.lb = None with self.assertRaisesRegex(RuntimeError, 'LB not activated'): obs.calculate() def test_lb_node_set_get(self): lbf = self.lb_class(kT=0.0, ext_force_density=[0, 0, 0], **self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf self.assertAlmostEqual( lbf[0, 0, 0].density, self.params['density'], delta=1e-4) @@ -389,24 +388,27 @@ def test_lb_node_set_get(self): def test_parameter_change_without_seed(self): lbf = self.lb_class(kT=1.0, seed=42, **self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf self.system.thermostat.set_lb(LB_fluid=lbf, seed=23, gamma=2.0) self.system.thermostat.set_lb(LB_fluid=lbf, gamma=3.0) - if espressomd.has_features("ELECTROSTATICS"): - actor = espressomd.electrostatics.DH( - prefactor=1., kappa=1., r_cut=1.) with self.assertRaisesRegex(Exception, "Temperature change not supported by LB"): self.system.thermostat.turn_off() with self.assertRaisesRegex(Exception, "Time step change not supported by LB"): self.system.time_step /= 2. if espressomd.has_features("ELECTROSTATICS"): - with self.assertRaisesRegex(RuntimeError, "LB does not currently support handling changes of the MD cell geometry"): - self.system.actors.add(actor) - self.assertEqual(len(self.system.actors), 1) + self.system.electrostatics.solver = espressomd.electrostatics.DH( + prefactor=1., kappa=1., r_cut=1.) # should not fail + self.system.electrostatics.clear() + with self.assertRaisesRegex(RuntimeError, "MD cell geometry change not supported by LB"): + self.system.box_l = [1., 2., 3.] + np.testing.assert_allclose( + np.copy(self.system.box_l), [1., 2., 3.], atol=1e-7) + with self.assertRaisesRegex(RuntimeError, "MPI topology change not supported by LB"): + self.system.cell_system.node_grid = self.system.cell_system.node_grid def test_grid_index(self): lbf = self.lb_class(**self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf # check ranges and out-of-bounds access shape = lbf.shape for i in range(3): @@ -450,14 +452,14 @@ def test_agrid_rounding(self): system.box_l = [l] * 3 * np.array(system.cell_system.node_grid) lbf = self.lb_class(agrid=l / 31, density=1, kinematic_viscosity=1, kT=0, tau=system.time_step, **self.lb_params) - system.actors.add(lbf) + system.lb = lbf system.integrator.run(steps=1) - system.actors.clear() + system.lb = None system.box_l = old_l def test_bool_operations_on_node(self): lbf = self.lb_class(kT=1.0, seed=42, **self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf # test __eq()__ where a node is equal to itself and not equal to any # other node assert lbf[0, 0, 0] == lbf[0, 0, 0] @@ -476,7 +478,7 @@ def test_bool_operations_on_node(self): @utx.skipIfMissingFeatures("EXTERNAL_FORCES") def test_viscous_coupling(self): lbf = self.lb_class(**self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf self.system.thermostat.set_lb(LB_fluid=lbf, seed=3, gamma=self.gamma) # Random velocities @@ -513,7 +515,7 @@ def test_viscous_coupling(self): def test_viscous_coupling_pairs(self): lbf = self.lb_class(**self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf self.system.thermostat.set_lb(LB_fluid=lbf, seed=3, gamma=self.gamma) # Random velocities @@ -572,7 +574,7 @@ def test_thermalization_force_balance(self): system.part.all().mass = 0.1 + np.random.random(len(system.part)) lbf = self.lb_class(kT=1.5, seed=4, **self.params, **self.lb_params) - system.actors.add(lbf) + system.lb = lbf system.thermostat.set_lb(LB_fluid=lbf, seed=3, gamma=self.gamma) for _ in range(20): @@ -586,7 +588,7 @@ def test_thermalization_force_balance(self): def test_force_interpolation(self): lbf = self.lb_class(**self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf self.system.thermostat.set_lb(LB_fluid=lbf, seed=3, gamma=self.gamma) position = np.array([1., 2., 3.]) @@ -612,7 +614,7 @@ def test_ext_force_density(self): ext_force_density = [2.3, 1.2, 0.1] lbf = self.lb_class(ext_force_density=ext_force_density, **self.params, **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf n_time_steps = 1 self.system.integrator.run(n_time_steps) # ext_force_density is a force density, therefore v = ext_force_density @@ -656,40 +658,40 @@ def params_with_tau(tau): lbf = self.lb_class(**params_with_tau(self.system.time_step), **self.lb_params) sim_time = 100 * self.params['tau'] - self.system.actors.add(lbf) + self.system.lb = lbf self.system.thermostat.set_lb(LB_fluid=lbf, gamma=0.1) self.system.integrator.run( int(round(sim_time / self.system.time_step))) probe_pos = np.array(self.system.box_l) / 2. v1 = np.copy(lbf.get_interpolated_velocity(pos=probe_pos)) f1 = np.copy(p.f) - self.system.actors.clear() + self.system.lb = None # get fresh LBfluid and change time steps with self.assertRaises(Exception): - self.system.actors.add( + self.system.lb = \ self.lb_class(**params_with_tau(0.5 * self.system.time_step), - **self.lb_params)) - self.system.actors.clear() + **self.lb_params) + self.system.lb = None with self.assertRaises(Exception): - self.system.actors.add( + self.system.lb = \ self.lb_class(**params_with_tau(1.1 * self.system.time_step), - **self.lb_params)) - self.system.actors.clear() + **self.lb_params) + self.system.lb = None - self.system.actors.add( + self.system.lb = \ self.lb_class(**params_with_tau(self.system.time_step), - **self.lb_params)) + **self.lb_params) with self.assertRaisesRegex(ValueError, r"LB tau \(0\.0100[0-9]+\) must be >= MD time_step \(0\.0200[0-9]+\)"): self.system.time_step = 2.0 * lbf.get_params()["tau"] with self.assertRaisesRegex(ValueError, r"LB tau \(0\.0100[0-9]+\) must be an integer multiple of the MD time_step \(0\.0080[0-9]+\)"): self.system.time_step = 0.8 * lbf.get_params()["tau"] - self.system.actors.clear() + self.system.lb = None self.system.time_step = 0.5 * self.params['tau'] lbf = self.lb_class(**params_with_tau(self.system.time_step), **self.lb_params) - self.system.actors.add(lbf) + self.system.lb = lbf self.system.integrator.run( int(round(sim_time / self.system.time_step))) v2 = np.copy(lbf.get_interpolated_velocity(pos=probe_pos)) diff --git a/testsuite/python/lb_boundary.py b/testsuite/python/lb_boundary.py index 5f1bf61dd83..e40de2a91c2 100644 --- a/testsuite/python/lb_boundary.py +++ b/testsuite/python/lb_boundary.py @@ -36,10 +36,10 @@ def setUp(self): self.lbf = self.lb_class( kinematic_viscosity=1.0, density=1.0, agrid=0.5, tau=1.0, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf def tearDown(self): - self.system.actors.clear() + self.system.lb = None def check_boundary_flags(self, slip_velocity1, slip_velocity2): def vbb2vel(values): diff --git a/testsuite/python/lb_boundary_velocity.py b/testsuite/python/lb_boundary_velocity.py index c5374abd1eb..7220851eeee 100644 --- a/testsuite/python/lb_boundary_velocity.py +++ b/testsuite/python/lb_boundary_velocity.py @@ -40,11 +40,11 @@ class LBBoundaryVelocityTest(ut.TestCase): system.cell_system.skin = 0.1 def tearDown(self): - self.system.actors.clear() + self.system.lb = None def setUp(self): self.lb_fluid = espressomd.lb.LBFluidWalberla(**self.lb_params) - self.system.actors.add(self.lb_fluid) + self.system.lb = self.lb_fluid def check_wall_slip(self, v_boundary, atol): """ diff --git a/testsuite/python/lb_boundary_volume_force.py b/testsuite/python/lb_boundary_volume_force.py index 26192c73a34..76beda388fe 100644 --- a/testsuite/python/lb_boundary_volume_force.py +++ b/testsuite/python/lb_boundary_volume_force.py @@ -48,10 +48,10 @@ class LBBoundaryForceCommon: def setUp(self): self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf def tearDown(self): - self.system.actors.clear() + self.system.lb = None def test(self): """ diff --git a/testsuite/python/lb_buoyancy_force.py b/testsuite/python/lb_buoyancy_force.py index 6e0f6b3e55e..0b04efb552f 100644 --- a/testsuite/python/lb_buoyancy_force.py +++ b/testsuite/python/lb_buoyancy_force.py @@ -52,10 +52,10 @@ class LBBuoyancy: def setUp(self): self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf def tearDown(self): - self.system.actors.clear() + self.system.lb = None def test(self): # Setup walls diff --git a/testsuite/python/lb_circular_couette.py b/testsuite/python/lb_circular_couette.py index fb90a055a20..71f58366712 100644 --- a/testsuite/python/lb_circular_couette.py +++ b/testsuite/python/lb_circular_couette.py @@ -51,7 +51,7 @@ class LBCircularCouetteCommon: system.periodicity = [False, False, True] def tearDown(self): - self.system.actors.clear() + self.system.lb = None def test_taylor_couette_flow(self): """ @@ -64,7 +64,7 @@ def test_taylor_couette_flow(self): lb_fluid = espressomd.lb.LBFluidWalberla( agrid=AGRID, density=0.5, kinematic_viscosity=3.2, tau=system.time_step) - system.actors.add(lb_fluid) + self.system.lb = lb_fluid # set up two cylinders cyl_center = AGRID * (GRID_SIZE // 2 + 0.5) * [1, 1, 0] diff --git a/testsuite/python/lb_electrohydrodynamics.py b/testsuite/python/lb_electrohydrodynamics.py index 500fb7f4ec6..e8d7808f3c2 100644 --- a/testsuite/python/lb_electrohydrodynamics.py +++ b/testsuite/python/lb_electrohydrodynamics.py @@ -52,11 +52,11 @@ def setUp(self): kT=self.params['temp'] ) - system.actors.add(lbf) + self.system.lb = lbf system.thermostat.set_lb(LB_fluid=lbf, gamma=self.params['friction']) def tearDown(self): - self.system.actors.clear() + self.system.lb = None self.system.part.clear() def test(self): diff --git a/testsuite/python/lb_interpolation.py b/testsuite/python/lb_interpolation.py index cee0fa86601..267886eaafe 100644 --- a/testsuite/python/lb_interpolation.py +++ b/testsuite/python/lb_interpolation.py @@ -58,10 +58,10 @@ class LBInterpolation: def setUp(self): self.lbf = self.lb_class(**LB_PARAMETERS, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf def tearDown(self): - self.system.actors.clear() + self.system.lb = None def set_boundaries(self, velocity): """Place boundaries *not* exactly on a LB node.""" diff --git a/testsuite/python/lb_lees_edwards.py b/testsuite/python/lb_lees_edwards.py index 3e0273497c2..9fd79a63665 100644 --- a/testsuite/python/lb_lees_edwards.py +++ b/testsuite/python/lb_lees_edwards.py @@ -42,12 +42,12 @@ def __init__(self, **kwargs): def __enter__(self): self.lbf = espressomd.lb.LBFluidWalberla( agrid=1., density=1., kinematic_viscosity=1., tau=system.time_step, **self.kwargs) - system.actors.add(self.lbf) + system.lb = self.lbf system.thermostat.set_lb(LB_fluid=self.lbf, gamma=1.0) return self.lbf def __exit__(self, exc_type, exc_value, exc_traceback): - system.actors.remove(self.lbf) + system.lb = None system.thermostat.turn_off() @@ -87,7 +87,7 @@ def setUp(self): protocol=espressomd.lees_edwards.Off()) def tearDown(self): - system.actors.clear() + system.lb = None system.thermostat.turn_off() system.part.clear() @@ -192,6 +192,12 @@ def test_velocity_shift_from_particles(self): for profile in self.sample_lb_velocities(lbf): self.check_profile(profile, stencil_D2Q8, '', 'SNWE', tol) + # order of instantiation doesn't matter + with LBContextManager() as lbf: + with LEContextManager('x', 'y', 0): + for profile in self.sample_lb_velocities(lbf): + self.check_profile(profile, stencil_D2Q8, '', 'SNWE', tol) + le_offset = 6 # North and South are sheared horizontally @@ -248,6 +254,13 @@ def create_impulse(lbf, stencil): for profile in self.sample_lb_velocities(lbf): self.check_profile(profile, stencil_D2Q8, '', 'SNWE', tol) + # order of instantiation doesn't matter + with LBContextManager() as lbf: + with LEContextManager('x', 'y', 0): + create_impulse(lbf, stencil_D2Q8) + for profile in self.sample_lb_velocities(lbf): + self.check_profile(profile, stencil_D2Q8, '', 'SNWE', tol) + le_offset = 6 # North and South are sheared horizontally @@ -276,10 +289,12 @@ def test_lebc_mismatch(self): Check that MD LEbc and LB LEbc always agree. """ err_msg = "MD and LB Lees-Edwards boundary conditions disagree" - # LEbc must be set before instantiating LB + # MD and LB LEbc must match with self.assertRaisesRegex(RuntimeError, err_msg): with LBContextManager() as lbf: LEContextManager('y', 'x', 1.).initialize() + with LBContextManager() as lbf: + LEContextManager('x', 'y', 1.).initialize() # when a LB actor with LEbc is active, the MD LEbc shear directions # are immutable with LEContextManager('x', 'y', 1.): @@ -296,37 +311,37 @@ def test_lebc_mismatch(self): # the MD LEbc must have the same shear directions with self.assertRaisesRegex(Exception, err_msg): with LEContextManager('z', 'y', 1.): - system.actors.add(lbf) - self.assertEqual(len(system.actors), 0) + system.lb = lbf + self.assertIsNone(system.lb) # LB only implements shear_plane_normal="y" err_msg = "Lees-Edwards LB only supports shear_plane_normal=\"y\"" for shear_dir, shear_plane_normal in itertools.product("xyz", "xz"): if shear_dir != shear_plane_normal: with self.assertRaisesRegex(ValueError, err_msg): with LEContextManager(shear_dir, shear_plane_normal, 1.): - system.actors.add(espressomd.lb.LBFluidWalberla( + system.lb = espressomd.lb.LBFluidWalberla( agrid=1., density=1., kinematic_viscosity=1., - tau=system.time_step)) - self.assertEqual(len(system.actors), 0) + tau=system.time_step) + self.assertIsNone(system.lb) # while LB and MD LEbc must agree on the shear directions, # the offset can change with LEContextManager('x', 'y', -1.): - system.actors.add(lbf) - system.actors.clear() + system.lb = lbf + system.lb = None # no thermalization with self.assertRaisesRegex(RuntimeError, "Lees-Edwards LB doesn't support thermalization"): with LEContextManager('x', 'y', 1.): - system.actors.add(espressomd.lb.LBFluidWalberla( + system.lb = espressomd.lb.LBFluidWalberla( agrid=1., density=1., kinematic_viscosity=1., kT=1., seed=42, - tau=system.time_step)) - self.assertEqual(len(system.actors), 0) + tau=system.time_step) + self.assertIsNone(system.lb) with self.assertRaisesRegex(ValueError, "Lees-Edwards sweep is implemented for a ghost layer of thickness 1"): lattice = espressomd.lb.LatticeWalberla(agrid=1., n_ghost_layers=2) with LEContextManager('x', 'y', 1.): - system.actors.add(espressomd.lb.LBFluidWalberla( + system.lb = espressomd.lb.LBFluidWalberla( lattice=lattice, density=1., kinematic_viscosity=1., - tau=system.time_step)) + tau=system.time_step) if __name__ == "__main__": diff --git a/testsuite/python/lb_lees_edwards_particle_coupling.py b/testsuite/python/lb_lees_edwards_particle_coupling.py index 373d2eca250..103197c721f 100644 --- a/testsuite/python/lb_lees_edwards_particle_coupling.py +++ b/testsuite/python/lb_lees_edwards_particle_coupling.py @@ -42,7 +42,7 @@ def test(self): lbf = espressomd.lb.LBFluidWalberla( agrid=1., density=1., kinematic_viscosity=1., tau=system.time_step) - system.actors.add(lbf) + system.lb = lbf system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1) pos = [system.box_l[0] / 2., 0., system.box_l[0] / 2.] diff --git a/testsuite/python/lb_mass_conservation.py b/testsuite/python/lb_mass_conservation.py index 22dcbb6a5a1..cf8ee65c500 100644 --- a/testsuite/python/lb_mass_conservation.py +++ b/testsuite/python/lb_mass_conservation.py @@ -47,11 +47,11 @@ class LBMassCommon: def setUp(self): self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf self.system.thermostat.set_lb(LB_fluid=self.lbf, seed=3, gamma=2.0) def tearDown(self): - self.system.actors.clear() + self.system.lb = None self.system.thermostat.turn_off() def test_mass_conservation(self): diff --git a/testsuite/python/lb_momentum_conservation.py b/testsuite/python/lb_momentum_conservation.py index 0abb38b5030..558185a911b 100644 --- a/testsuite/python/lb_momentum_conservation.py +++ b/testsuite/python/lb_momentum_conservation.py @@ -58,12 +58,12 @@ def setUp(self): self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params) def tearDown(self): - self.system.actors.clear() + self.system.lb = None self.system.thermostat.turn_off() self.system.part.clear() def test(self): - self.system.actors.add(self.lbf) + self.system.lb = self.lbf self.system.thermostat.set_lb(LB_fluid=self.lbf, gamma=GAMMA, seed=1) np.testing.assert_allclose( self.lbf.ext_force_density, diff --git a/testsuite/python/lb_planar_couette.py b/testsuite/python/lb_planar_couette.py index 34ed3b08da1..dfdabd4eb72 100644 --- a/testsuite/python/lb_planar_couette.py +++ b/testsuite/python/lb_planar_couette.py @@ -69,7 +69,7 @@ def setUp(self): self.system.time = 0. def tearDown(self): - self.system.actors.clear() + self.system.lb = None self.system.lees_edwards = espressomd.lees_edwards.LeesEdwards() def check_profile(self, u_getter, **kwargs): @@ -89,7 +89,7 @@ def check_profile(self, u_getter, **kwargs): protocol=protocol, **kwargs) lbf = self.lb_class(**LB_PARAMS, **self.lb_params) - system.actors.add(lbf) + self.system.lb = lbf # warmup system.integrator.run(8) diff --git a/testsuite/python/lb_poiseuille.py b/testsuite/python/lb_poiseuille.py index bb3caf9f728..2176263ed8b 100644 --- a/testsuite/python/lb_poiseuille.py +++ b/testsuite/python/lb_poiseuille.py @@ -68,10 +68,10 @@ class LBPoiseuilleCommon: def setUp(self): self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf def tearDown(self): - self.system.actors.clear() + self.system.lb = None def prepare(self): """ diff --git a/testsuite/python/lb_poiseuille_cylinder.py b/testsuite/python/lb_poiseuille_cylinder.py index c3a95ebab9f..ccc1fa0519c 100644 --- a/testsuite/python/lb_poiseuille_cylinder.py +++ b/testsuite/python/lb_poiseuille_cylinder.py @@ -83,7 +83,7 @@ class LBPoiseuilleCommon: 'orientation': [1, 0, 0]} def tearDown(self): - self.system.actors.clear() + self.system.lb = None def prepare(self): """ @@ -98,7 +98,7 @@ def prepare(self): local_lb_params['ext_force_density'] = np.array( self.params['axis']) * EXT_FORCE self.lbf = self.lb_class(**local_lb_params, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf cylinder_shape = espressomd.shapes.Cylinder( center=self.system.box_l / 2.0, axis=self.params['axis'], diff --git a/testsuite/python/lb_pressure_tensor.py b/testsuite/python/lb_pressure_tensor.py index 6446081867e..21f43797567 100644 --- a/testsuite/python/lb_pressure_tensor.py +++ b/testsuite/python/lb_pressure_tensor.py @@ -44,18 +44,17 @@ class TestLBPressureTensor: system.cell_system.skin = 0 def tearDown(self): - self.system.actors.clear() + self.system.lb = None self.system.thermostat.turn_off() def setUp(self): # Setup - system = self.system self.lbf = self.lb_class(**self.params, **self.lb_params) - system.actors.add(self.lbf) - system.thermostat.set_lb(LB_fluid=self.lbf, seed=42) + self.system.lb = self.lbf + self.system.thermostat.set_lb(LB_fluid=self.lbf, seed=42) # Warmup - system.integrator.run(500) + self.system.integrator.run(500) # Sampling self.p_global = np.zeros((self.steps, 3, 3)) @@ -71,7 +70,7 @@ def setUp(self): self.p_node1[i] = node1.pressure_tensor self.p_global[i] = self.lbf.pressure_tensor - system.integrator.run(2) + self.system.integrator.run(2) def assert_allclose_matrix(self, x, y, atol_diag, atol_offdiag): """Assert that all elements x_ij, y_ij are close with diff --git a/testsuite/python/lb_shear.py b/testsuite/python/lb_shear.py index fcb95eb16c3..fef0838ba6f 100644 --- a/testsuite/python/lb_shear.py +++ b/testsuite/python/lb_shear.py @@ -88,7 +88,7 @@ def setUp(self): self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params) def tearDown(self): - self.system.actors.clear() + self.system.lb = None def check_profile(self, shear_plane_normal, shear_direction): """ @@ -99,7 +99,7 @@ def check_profile(self, shear_plane_normal, shear_direction): self.system.box_l = np.max( ((W, W, W), shear_plane_normal * (H + 2 * AGRID)), 0) self.setUp() - self.system.actors.add(self.lbf) + self.system.lb = self.lbf self.lbf.clear_boundaries() wall_shape1 = espressomd.shapes.Wall( diff --git a/testsuite/python/lb_slice.py b/testsuite/python/lb_slice.py index b5a46367405..840d2501189 100644 --- a/testsuite/python/lb_slice.py +++ b/testsuite/python/lb_slice.py @@ -40,10 +40,10 @@ def setUp(self): self.lb_fluid = self.lb_class( agrid=1., density=1., kinematic_viscosity=1., tau=self.system.time_step, **self.lb_params) - self.system.actors.add(self.lb_fluid) + self.system.lb = self.lb_fluid def tearDown(self): - self.system.actors.clear() + self.system.lb = None def test_slicing(self): lb_fluid = self.lb_fluid diff --git a/testsuite/python/lb_stats.py b/testsuite/python/lb_stats.py index 1687cae6220..1f5044a0eca 100644 --- a/testsuite/python/lb_stats.py +++ b/testsuite/python/lb_stats.py @@ -47,7 +47,7 @@ class TestLB: n_nodes = system.cell_system.get_state()["n_nodes"] def tearDown(self): - self.system.actors.clear() + self.system.lb = None self.system.part.clear() self.system.thermostat.turn_off() @@ -70,7 +70,7 @@ def test_mass_momentum_thermostat(self): agrid=self.params['agrid'], tau=self.system.time_step, ext_force_density=[0, 0, 0], seed=4) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf self.system.thermostat.set_lb( LB_fluid=self.lbf, seed=3, diff --git a/testsuite/python/lb_stokes_sphere.py b/testsuite/python/lb_stokes_sphere.py index dc33a9a7bfe..b4c12c6d219 100644 --- a/testsuite/python/lb_stokes_sphere.py +++ b/testsuite/python/lb_stokes_sphere.py @@ -60,11 +60,11 @@ class Stokes: def setUp(self): self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf self.system.thermostat.set_lb(LB_fluid=self.lbf, gamma=1.0) def tearDown(self): - self.system.actors.clear() + self.system.lb = None self.system.thermostat.turn_off() def test_stokes(self): diff --git a/testsuite/python/lb_streaming.py b/testsuite/python/lb_streaming.py index 3281ffc1768..1a5740cee74 100644 --- a/testsuite/python/lb_streaming.py +++ b/testsuite/python/lb_streaming.py @@ -96,10 +96,10 @@ class LBStreamingCommon: def setUp(self): self.lbf = self.lb_class(**LB_PARAMETERS, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf def tearDown(self): - self.system.actors.clear() + self.system.lb = None def test_population_streaming(self): pop_default = np.zeros(19) + 1e-10 diff --git a/testsuite/python/lb_thermo_virtual.py b/testsuite/python/lb_thermo_virtual.py index 5ce1474bd2a..ab9875f4df4 100644 --- a/testsuite/python/lb_thermo_virtual.py +++ b/testsuite/python/lb_thermo_virtual.py @@ -39,44 +39,42 @@ class LBBoundaryThermoVirtualTest(ut.TestCase): system.cell_system.skin = 0.1 def tearDown(self): + self.system.lb = None self.system.part.clear() - for a in self.system.actors: - self.system.actors.remove(a) - def check_virtual(self, fluid_class): - s = self.system + system = self.system lb_fluid = fluid_class( agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=1.0) - s.actors.add(lb_fluid) + system.lb = lb_fluid - virtual = s.part.add(pos=[0, 0, 0], virtual=True, v=[1, 0, 0]) - physical = s.part.add(pos=[0, 0, 0], virtual=False, v=[1, 0, 0]) + virtual = system.part.add(pos=[0, 0, 0], virtual=True, v=[1, 0, 0]) + physical = system.part.add(pos=[0, 0, 0], virtual=False, v=[1, 0, 0]) - s.thermostat.set_lb( + system.thermostat.set_lb( LB_fluid=lb_fluid, act_on_virtual=False, gamma=1.0) - s.integrator.run(1) + system.integrator.run(1) np.testing.assert_almost_equal(np.copy(virtual.f), [0, 0, 0]) np.testing.assert_almost_equal(np.copy(physical.f), [-1, 0, 0]) - s.thermostat.set_lb(LB_fluid=lb_fluid, act_on_virtual=True) + system.thermostat.set_lb(LB_fluid=lb_fluid, act_on_virtual=True) virtual.v = [1, 0, 0] physical.v = [1, 0, 0] - s.actors.remove(lb_fluid) + system.lb = None lb_fluid = fluid_class( agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=1.0) - s.actors.add(lb_fluid) - s.thermostat.set_lb(LB_fluid=lb_fluid, gamma=1.0) + system.lb = lb_fluid + system.thermostat.set_lb(LB_fluid=lb_fluid, gamma=1.0) virtual.pos = physical.pos virtual.v = [1, 0, 0] physical.v = [1, 0, 0] - s.integrator.run(1) + system.integrator.run(1) # The forces are not exactly -1 because the fluid is not at # rest anymore because of the previous check. diff --git a/testsuite/python/lb_thermostat.py b/testsuite/python/lb_thermostat.py index 516dece5ffe..046b4801db6 100644 --- a/testsuite/python/lb_thermostat.py +++ b/testsuite/python/lb_thermostat.py @@ -67,12 +67,12 @@ class LBThermostatCommon(thermostats_common.ThermostatsCommon): def setUp(self): self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf self.system.thermostat.set_lb( LB_fluid=self.lbf, seed=5, gamma=self.global_gamma) def tearDown(self): - self.system.actors.clear() + self.system.lb = None self.system.thermostat.turn_off() self.system.part.clear() diff --git a/testsuite/python/linear_momentum_lb.py b/testsuite/python/linear_momentum_lb.py index 5994ffc48d8..1bd4beb1142 100644 --- a/testsuite/python/linear_momentum_lb.py +++ b/testsuite/python/linear_momentum_lb.py @@ -53,10 +53,10 @@ class LBLinearMomentum: def setUp(self): self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf def tearDown(self): - self.system.actors.clear() + self.system.lb = None def test(self): """ diff --git a/testsuite/python/long_range_actors.py b/testsuite/python/long_range_actors.py index 7b37911c2bb..d8ac73e4f1d 100644 --- a/testsuite/python/long_range_actors.py +++ b/testsuite/python/long_range_actors.py @@ -128,7 +128,7 @@ def test_magnetostatics_registration(self): self.assertIsNone(dp3m.call_method("unknown")) self.system.magnetostatics.solver = dp3m - self.assertIsNone(self.system.electrostatics.call_method("unknown")) + self.assertIsNone(self.system.magnetostatics.call_method("unknown")) def check_obs_stats(self, key): # check observable statistics have the correct shape diff --git a/testsuite/python/observable_cylindricalLB.py b/testsuite/python/observable_cylindricalLB.py index dd9f31637d5..60423580898 100644 --- a/testsuite/python/observable_cylindricalLB.py +++ b/testsuite/python/observable_cylindricalLB.py @@ -65,10 +65,10 @@ class CylindricalLBObservableCommon: def setUp(self): self.lbf = self.lb_class(**self.lb_params, **self.lb_params_extra) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf def tearDown(self): - self.system.actors.clear() + self.system.lb = None self.system.part.clear() def calc_vel_at_pos(self, positions): diff --git a/testsuite/python/observable_profileLB.py b/testsuite/python/observable_profileLB.py index ac810b24626..c2dbc8f1efa 100644 --- a/testsuite/python/observable_profileLB.py +++ b/testsuite/python/observable_profileLB.py @@ -71,10 +71,10 @@ class ObservableProfileLBCommon: def setUp(self): self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf def tearDown(self): - self.system.actors.clear() + self.system.lb = None def set_fluid_velocities(self): """Set an x dependent fluid velocity.""" @@ -105,7 +105,7 @@ def test_velocity_profile(self): LB_VELOCITY_PROFILE_PARAMS['n_z_bins'] * 3) def test_error_if_no_LB(self): - self.system.actors.clear() + self.system.lb = None obs = espressomd.observables.LBVelocityProfile( **LB_VELOCITY_PROFILE_PARAMS) with self.assertRaises(RuntimeError): diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index 9981fdc8f92..1886a5c862d 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -93,9 +93,9 @@ lattice=lb_lattice, density=1.5, kT=2.0, diffusion=0.2, valency=0.1, advection=False, friction_coupling=False, ext_efield=[0.1, 0.2, 0.3], single_precision=False, tau=system.time_step) - system.ekcontainer.solver = ek_solver - system.ekcontainer.tau = ek_species.tau - system.ekcontainer.add(ek_species) + ekcontainer = espressomd.electrokinetics.EKContainer( + solver=ek_solver, tau=ek_species.tau) + ekcontainer.add(ek_species) ek_species.add_boundary_from_shape( shape=wall1, value=1e-3 * np.array([1., 2., 3.]), boundary_type=espressomd.electrokinetics.FluxBoundary) @@ -358,7 +358,8 @@ "p2nfft_alpha": "0.37"}) if lbf_class: - system.actors.add(lbf) + system.lb = lbf + system.ekcontainer = ekcontainer if 'THERM.LB' in modes: system.thermostat.set_lb(LB_fluid=lbf, seed=23, gamma=2.0) # Create a 3D grid with deterministic values to fill the LB fluid lattice @@ -505,7 +506,7 @@ def test_lb_checkpointing_exceptions(self): lbf.save_checkpoint(str(lbf_cpt_root / "lb_err.cpt"), 2) # deactivate LB actor - system.actors.remove(lbf) + system.lb = None # read the valid LB checkpoint file lbf_cpt_data = lbf_cpt_path.read_bytes() diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 21094ae46d4..d9effd289b5 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -69,24 +69,10 @@ def setUpClass(cls): 'THERM.SDM' in modes or 'INT.SDM' in modes): cls.ref_periodicity = np.array([False, False, False]) - def get_active_actor_of_type(self, actor_type): - for actor in system.actors.active_actors: - if isinstance(actor, actor_type): - return actor - self.fail( - f"system doesn't have an actor of type {actor_type.__name__}") - - def test_get_active_actor_of_type(self): - if system.actors.active_actors: - actor = system.actors.active_actors[0] - self.assertEqual(self.get_active_actor_of_type(type(actor)), actor) - with self.assertRaisesRegex(AssertionError, "system doesn't have an actor of type Wall"): - self.get_active_actor_of_type(espressomd.shapes.Wall) - @utx.skipIfMissingFeatures(["WALBERLA"]) @ut.skipIf(not has_lb_mode, "Skipping test due to missing LB feature.") def test_lb_fluid(self): - lbf = self.get_active_actor_of_type(espressomd.lb.LBFluidWalberla) + lbf = system.lb cpt_mode = 0 if 'LB.ASCII' in modes else 1 cpt_root = pathlib.Path(self.checkpoint.checkpoint_dir) cpt_path = str(cpt_root / "lb") + "{}.cpt" @@ -177,8 +163,6 @@ def test_ek_species(self): self.assertEqual(len(system.ekcontainer), 1) ek_species = system.ekcontainer[0] - self.assertTrue( - system.ekcontainer.call_method("is_poisson_solver_set")) self.assertAlmostEqual(system.ekcontainer.tau, system.time_step, delta=1e-7) self.assertIsInstance(system.ekcontainer.solver, @@ -277,7 +261,7 @@ def generator(value, shape): @utx.skipIfMissingFeatures(["WALBERLA"]) @ut.skipIf(not has_lb_mode, "Skipping test due to missing LB feature.") def test_lb_vtk(self): - lbf = self.get_active_actor_of_type(espressomd.lb.LBFluidWalberla) + lbf = system.lb self.assertEqual(len(lbf.vtk_writers), 2) vtk_suffix = config.test_name key_auto = f"vtk_out/auto_lb_{vtk_suffix}" diff --git a/testsuite/python/virtual_sites_tracers_common.py b/testsuite/python/virtual_sites_tracers_common.py index 869f2b7a1a2..3de7dfa0d8b 100644 --- a/testsuite/python/virtual_sites_tracers_common.py +++ b/testsuite/python/virtual_sites_tracers_common.py @@ -39,16 +39,16 @@ def setUp(self): self.system.box_l = (self.box_lw, self.box_lw, self.box_height) def tearDown(self): - self.system.actors.clear() + self.system.lb = None self.system.part.clear() self.system.thermostat.turn_off() def set_lb(self, ext_force_density=(0, 0, 0), dir_walls=2): - self.system.actors.clear() + self.system.lb = None self.lbf = self.LBClass( kT=0.0, agrid=1., density=1., kinematic_viscosity=1.8, tau=self.system.time_step, ext_force_density=ext_force_density) - self.system.actors.add(self.lbf) + self.system.lb = self.lbf self.system.thermostat.set_lb( LB_fluid=self.lbf, act_on_virtual=False, @@ -160,7 +160,7 @@ def test_advection(self): tracer_dist = p.pos[direction] - pos_initial[direction] self.assertAlmostEqual(tracer_dist / ref_dist, 1., delta=0.01) - system.actors.clear() + system.lb = None def test_zz_exceptions_without_lb(self): """Check behaviour without lb. Ignore non-virtual particles, complain on @@ -170,7 +170,7 @@ def test_zz_exceptions_without_lb(self): self.set_lb() system = self.system system.virtual_sites = espressomd.virtual_sites.VirtualSitesInertialessTracers() - system.actors.clear() + system.lb = None system.part.clear() p = system.part.add(pos=(0, 0, 0)) system.integrator.run(1)