Rewrite LB/EK interface (#4773)

Description of changes: - API changes: - the `espressomd.System` class now has a new member `lb` - the `espressomd.System` class member `ekreactions` is now the `reaction` attribute of `system.ekcontainer` - the `espressomd.System` class no longer has an `actors` member - further API changes will take place in an upcoming PR - remove Boost dependency from LB/EK sources - follow-up to [**Remove remaining boost components** walberla/walberla#605](https://i10git.cs.fau.de/walberla/walberla/-/merge_requests/605) - partial fix for #4723 (comment) - unit tests from the walberla bridge still depend on `Boost::unit_test_framework` - fix regressions in EK - EK now throws an error when the box_l or node_grid changes, when the MD time step and EK tau disagree, when the box size is not an integer multiple of agrid - rewrite the LB and EK code in the core from scratch to enforce SOLID principles - folder `src/grid_based_algorithms` was split into `src/core/lb` and `src/core/ek` - LB and EK global variables (`lattice_switch`, `ek_container`, `ek_reactions`, `lb_walberla_instance`, `lb_walberla_params_instance`) were removed; these objects are now members of the `System` class (partial fix for #2628)
espressomd · Aug 22, 2023 · 08869c6 · 08869c6
2 parents 4aa06e7 + 947ddc8
commit 08869c6
Show file tree

Hide file tree

Showing 183 changed files with 2,987 additions and 2,211 deletions.
diff --git 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
@@ -161,29 +161,34 @@ 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
 ^^^^^^^^^^^^^^^^^
 ::
 
     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,15 +334,15 @@ 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
         vtk_grids = vtk_reader.parse(path_vtk)
         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/electrostatics.rst b/doc/sphinx/electrostatics.rst
@@ -27,8 +27,7 @@ interactions as efficiently as possible, but almost all of them require
 some knowledge to use them properly. Uneducated use can result in
 completely unphysical simulations.
 
-Coulomb interactions have to be added to the list of active actors of the system object to become
-active. This is done by calling the add-method of :attr:`espressomd.system.System.actors`.
+Coulomb interactions have to be attached to the system object to become active.
 Only one electrostatics method can be active at any time.
 
 Note that using the electrostatic interaction also requires assigning charges to
@@ -299,8 +298,8 @@ Usage notes:
     elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=box_l * 0.2, maxPWerror=1e-3)
     system.electrostatics.solver = elc
 
-Although it is technically feasible to remove ``elc`` from the list of actors
-and then to add the ``p3m`` object, it is not recommended because the P3M
+Although it is technically feasible to detach ``elc`` from the system
+and then to attach the ``p3m`` object, it is not recommended because the P3M
 parameters are mutated by *ELC*, e.g. the ``epsilon`` is made metallic.
 It is safer to instantiate a new P3M object instead of recycling one that
 has been adapted by *ELC*.

diff --git 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/io.rst b/doc/sphinx/io.rst
@@ -111,11 +111,10 @@ Be aware of the following limitations:
   for a specific combination of features, please share your findings
   with the |es| community.
 
-* The active actors, i.e., the content of ``system.actors`` resp.
-  ``system.ekcontainers``, are checkpointed. For lattice-based methods like
-  lattice-Boltzmann fluids and advection-diffusion-reaction models, this only
-  includes the parameters such as the lattice constant (``agrid``) and initial
-  densities.
+* All long-range solvers for electrostatics, magnetostatics, lattice-Boltzmann
+  and advection-diffusion-reaction are checkpointed. For lattice-Boltzmann
+  fluids and advection-diffusion-reaction models, this only includes the
+  parameters such as the lattice constant (``agrid``) and initial densities.
   The actual fields have to be saved separately with the lattice-specific
   methods :meth:`espressomd.lb.LBFluidWalberla.save_checkpoint
   <espressomd.detail.walberla.LatticeModel.save_checkpoint>` resp.
@@ -181,10 +180,10 @@ Be aware of the following limitations:
   This is because loading a checkpoint reinitializes the system and enforces
   a recalculation of the forces. However, this computes the forces from the
   velocities at the current time step and not at the previous half time step.
-  Please note that long-range actors can make trajectories non-reproducible.
+  Please note that long-range solvers can make trajectories non-reproducible.
   For example, lattice-Boltzmann introduces errors of the order of 1e-15 with
   binary checkpoint files, or 1e-7 with ASCII checkpoint files. In addition,
-  several electrostatic and magnetostatic actors automatically introduce
+  several electrostatic and magnetostatic solvers automatically introduce
   a deviation of the order of 1e-7, either due to floating-point rounding
   errors (:class:`~espressomd.electrostatics.P3MGPU`), or due to re-tuning
   using the most recent system state (:class:`~espressomd.electrostatics.MMM1D`,

diff --git 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/sphinx/magnetostatics.rst b/doc/sphinx/magnetostatics.rst
@@ -25,7 +25,7 @@ The prefactor :math:`D` is can be set by the user and is given by
 where :math:`\mu_0` and :math:`\mu` are the vacuum permittivity and the
 relative permittivity of the background material, respectively.
 
-Magnetostatic interactions are activated via the actor framework::
+Magnetostatic interactions are activated when attached to the system::
 
     import espressomd
     import espressomd.magnetostatics

diff --git 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/charged_system/NotesForTutor.md b/doc/tutorials/charged_system/NotesForTutor.md
@@ -17,5 +17,5 @@
 
 * Setting up WCA interactions (should be clear from previous tutorial)
 * Using particle properties like ``q``, ``fix`` and ``type``
-* Setting up electrostatic methods, the concept of ``actors``
+* Setting up electrostatic methods, the concept of system solver
 * There can only be one ESPResSo system. For multiple runs, the relevant parts of the system have to be reset manually
diff --git 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
@@ -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",
     "```"
    ]
   },