diff --git a/doc/tutorials/electrokinetics/CMakeLists.txt b/doc/tutorials/electrokinetics/CMakeLists.txt
index a3c419cb6fc..cf2099bada2 100644
--- a/doc/tutorials/electrokinetics/CMakeLists.txt
+++ b/doc/tutorials/electrokinetics/CMakeLists.txt
@@ -18,6 +18,6 @@
 #
 
 configure_tutorial_target(TARGET tutorial_ek DEPENDS electrokinetics.ipynb
-                          figures/schlitzpore_3d.png scripts/eof_analytical.py)
+                          figures/schlitzpore_3d.png)
 
 nb_export(TARGET tutorial_ek SUFFIX "" FILE "electrokinetics.ipynb" HTML_RUN)
diff --git a/doc/tutorials/electrokinetics/electrokinetics.ipynb b/doc/tutorials/electrokinetics/electrokinetics.ipynb
index cde7cc5e494..28510f0d347 100644
--- a/doc/tutorials/electrokinetics/electrokinetics.ipynb
+++ b/doc/tutorials/electrokinetics/electrokinetics.ipynb
@@ -4,346 +4,200 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Electrokinetics"
+    "# Electrokinetics\n",
+    "### Table of contents\n",
+    "1. [Introduction](#1.-Introduction)\n",
+    "2. [Advection-Diffusion in 2D](#2.-Advection-Diffusion-equation-in-2D)\n",
+    "3. [Electroosmotic flow](#3.-Electroosmotic-flow)\n",
+    "4. [Reaction in turbulent flow](#4.-Reaction-in-turbulent-flow)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Table of Contents\n",
-    "1. [Introduction](#Introduction)\n",
-    "2. [Theoretical Background](#Theoretical-Background)\n",
-    "   1. [The Electrokinetic Equations](#The-Electrokinetic-Equations)\n",
-    "   2. [EOF in the Slit Pore Geometry](#EOF-in-the-Slit-Pore-Geometry)\n",
-    "3. [Simulation using ESPResSo](#Simulation-using-ESPResSo)\n",
-    "   1. [Mapping SI and Simulation Units](#Mapping-SI-and-Simulation-Units)\n",
-    "   2. [Setting up the slit pore system](#Setting-up-the-slit-pore-system)\n",
-    "4. [References](#References)\n",
-    "    "
+    "## 1. Introduction"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Introduction\n",
+    "In this tutorial we're looking at the electrokinetics feature of ESPResSo, which allows us to describe the motion of potentially charged chemical species solvated in a fluid on a continuum level. The govering equations for the solvent are known as the Poisson-Nernst-Planck equations, which is the combination of the electrostatic Poisson equation and the dynamics of the chemical species described by the Nernst-Planck equation. For the advection we solve the incompressible Navier-Stokes equation. The total set of equations is given by\n",
     "\n",
-    "In recent years the lattice-Boltzmann method (LBM) has proven itself to be a viable way to introduce hydrodynamic interactions into coarse-grained MD simulations with moderate computational cost.\n",
-    "ESPResSo features such an algorithm, which can make use of the LBM and extend it to coarse-grain not only the solvent molecules but also ionic solutes. It is called EK and explicitly treats the ionic solutes in a continuum fashion and is valid for a wide range of salt concentrations [1-3].\n",
-    "\n",
-    "### Tutorial Outline\n",
-    "\n",
-    "To make our first steps using ELECTROKINETICS we will work on one of the few systems for which analytic solutions for the electrokinetic equations exist: the slip pore geometry with a counterion-only electrolyte.\n",
-    "The same slit pore system is also treated in the LBM tutorial, but there, the ionic species were modeled as explicit particles.\n",
-    "For this system, the two approaches lead to exactly the same results [4].\n",
-    "Differences became significant for multivalent ions, very high salt concentrations, and very high surface charge, since then the mean-field approach the EK employs, is basically solving the Poisson-Nernst-Planck formalism plus the Navier-Stokes equation on a lattice.\n",
-    "This leads to significantly different results from explicit ion approaches [5-7].\n",
-    "This tutorial is mainly divided into two sections.\n",
-    "* **Theoretical Background** introduces the electrokinetic equations and the analytical solution for the slit pore system.\n",
-    "* **Simulation using ESPResSo** deals exclusively with the simulation. \n",
+    "\\begin{aligned}\n",
+    "\\partial_{t} n_{i} &= - \\vec{\\nabla} \\cdot \\vec{j}_{i} \\\\\n",
+    "\\vec{j}_{i} &= - D_{i} \\vec{\\nabla} n_{i} - \\frac{z_{i} e}{k_{B} T} n_{i} \\vec{\\nabla} \\phi + n_{i} \\vec{u} \\\\\n",
+    "\\Delta \\phi &= \\frac{1}{4 \\pi \\varepsilon_{0} \\varepsilon_{\\mathrm{r}}} \\sum_{i} z_{i} e n_{i} \\\\\n",
+    "\\rho (\\partial_{t} \\vec{u} + (\\vec{u} \\cdot \\vec{\\nabla}) \\vec{u}) &= - \\vec{\\nabla} p + \\eta \\Delta \\vec{u} + \\sum_{i} \\frac{k_{B} T}{D_{i}} \\vec{j}_{i} + \\vec{f}_{\\mathrm{ext}} \\\\\n",
+    "\\vec{\\nabla} \\cdot \\vec{u} &= 0,\n",
+    "\\end{aligned}\n",
     "\n",
-    "If you already know about simple diffusion-migration-advection equations, continuum electrostatics, and Navier-Stokes, then you can skip the first section."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Theoretical Background\n"
+    "where $n_{i}$ denotes the ion density of species $i$, $\\vec{j}_{i}$ the density flux, $D_{i}$ the diffusion coefficient, $z_{i}$ the valency, $e$ the elementary charge, $k_{B}$ the Boltzmann constant, $T$ the temperature, $\\phi$ the electrostatic potential, $\\varepsilon_{0}$ the vacuum permittivity, $\\varepsilon_{\\mathrm{r}}$ the relative permittivity, $\\rho$ the fluid density, $\\vec{u}$ the fluid velocity, $p$ the hydrostatic pressure, $\\eta$ the dynamic viscosity, and $\\vec{f}_{\\mathrm{ext}}$ an external force density."
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### The Electrokinetic Equations\n",
-    "\n",
-    "In the following, we will derive the equations modeling the time evolution of the concentrations of dissolved species as well as the solvent in the standard electrokinetic model.\n",
-    "We do so, neglecting the salt ions' contribution to the overall mass density, which allows us to treat the dynamics of the ions and the fluid separately [7].\n",
-    "The solvent fluid will be modeled using the Navier-Stokes equations while we use a set of diffusion-migration-advection equations for the ionic species.\n"
+    "# 2. Advection-Diffusion equation in 2D"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "#### Ionic Species\n",
-    "\n",
-    "The description starts with the ionic species' concentrations $c_{k}(\\vec{r}, t)$ (number density) and the associated flux densities $\\vec{j}_{k}(\\vec{r}, t)$, for which mass conservation holds\n",
-    "\n",
-    "\\begin{equation}\n",
-    "\\partial_{t} c_{k} = -\\nabla \\cdot\\vec{j}_{k} . \n",
-    "\\end{equation}\n",
-    "\n",
-    "Here $\\vec{r}$ denotes the spatial coordinate and $t$ the time, while $k$ enumerates the ionic species.\n",
-    "The fluxes are caused by diffusion (due to density variations and external forces) and advection.\n",
-    "\n",
-    "The advective contribution to the flux is given by\n",
+    "The first system that is simulated in this tutorial is the simple advection-diffusion of a drop of uncharged chemical species in a constant velocity field. To keep the computation time small, we restrict ourselves to a 2D problem, but the algorithm is also capable of solving the 3D advection-diffusion equation. Furthermore, we can also skip solving the electrostatic Poisson equation, since there are is no charged species present. The equations we solve thus reduce to\n",
     "\n",
     "\\begin{equation}\n",
-    "\\vec{j}_{k}^{\\mathrm{adv.}} = c_{k} \\vec{u} ,\n",
+    "\\partial_{t} n = D \\Delta n - \\vec{\\nabla} \\cdot (\\vec{v} n).\n",
     "\\end{equation}\n",
     "\n",
-    "where $\\vec{u}(\\vec{r}, t)$ denotes the fluid velocity (advective velocity).\n",
-    "This equation models advection as a simple co-movement of the dissolved ions with the surrounding fluid.\n",
-    "All inertial effects of the ions are neglected.\n",
-    "\n",
-    "The diffusive behavior of the ions is best described in a reference frame co-moving with the local advective velocity $\\vec{u}$.\n",
-    "We assume that the species' relative fluxes instantaneously relax to a local equilibrium.\n",
-    "This assumption allows us to derive the diffusive fluxes from a local free-energy density, which we define as\n",
+    "The fundamental solution of this partial diffential equation can be found analytically in the case of a constant velocity field $\\vec{v}$ and a constant diffusion coefficient $D$. For a $d$-dimensional system, the solution of an initally infinitessimaly small droplet at the origin can be written as\n",
     "\n",
     "\\begin{equation}\n",
-    "f \\big( c_{k}(\\vec{r}) \\big) = \\sum_{k} \\underbrace{k_{\\mathrm{B}}T c_{k}(\\vec{r}) \\left[ \\log \\left\\lbrace \\Lambda_{k}^{3} c_{k}(\\vec{r}) \\right\\rbrace - 1 \\right] }_{\\mathrm{ideal~gas~contribution}} + \\underbrace{z_{k} e c_{k}(\\vec{r}) \\Phi(\\vec{r})}_{\\mathrm{electrostatic~contribution}} ,\n",
+    "n(\\vec{x},t) = \\frac{1}{(4 \\pi D t)^{d/2}} \\exp \\left( - \\frac{(\\vec{x} - \\vec{v} t)^2}{4 D t} \\right).\n",
     "\\end{equation}\n",
     "\n",
-    "with the $\\Lambda_{k}$ the species' thermal de Broglie wavelengths, $z_{k}$ their valencies, and $\\Phi(\\vec{r})$ the electrostatic potential.\n",
-    "This free-energy density consists of only an ideal-gas and an electrostatic contribution.\n",
-    "The same assumptions form the basis of Poisson-Boltzmann (PB) theory.\n",
-    "Hence, the limitations of this model are the same as those of PB.\n",
-    "That means this model applies to monovalent ions at low to intermediate densities and surface charges [5,6,10,11].\n",
-    "\n",
-    "The species' chemical potentials $\\mu_{k}$ implied by the free-energy density read\n",
-    "\n",
-    "\\begin{equation}\n",
-    "\\mu_{k}(\\vec{r}) = \\delta_{c_k} f(c_{k}\\big( \\vec r ) \\big) = k_{\\mathrm{B}}T \\log(\\Lambda_{k}^{3} c_{k}(\\vec{r})) + z_{k} e \\Phi(\\vec{r}) .\n",
-    "\\end{equation}\n",
+    "After importing the necessary packages, we start by defining the necessary parameters for the simulation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import espressomd\n",
+    "import espressomd.lb\n",
+    "import espressomd.electrokinetics\n",
+    "import espressomd.shapes\n",
     "\n",
-    "This in turn allows us to formulate the first-order approximation to the thermodynamic driving force as the gradient of the chemical potential, which we use to define an expression for the diffusive flux\n",
+    "espressomd.assert_features([\"WALBERLA\", \"WALBERLA_FFT\"])\n",
     "\n",
-    "\\begin{aligned}\n",
-    "\\vec{j}_{k}^{\\mathrm{diff}} &= \\xi_{k} \\left( -c_{k} \\nabla \\mu_{k} \\right) = -k_{\\mathrm{B}}T \\xi_{k} \\nabla c_{k} - \\xi_{k} z_{k} e c_{k} \\nabla\\Phi \\\\\n",
-    "&= -D_{k} \\nabla c_{k} - \\xi_{k} z_{k} e c_{k} \\nabla \\Phi .\n",
-    "\\end{aligned}\n",
+    "import tqdm\n",
+    "import numpy as np\n",
     "\n",
-    "Here, $\\xi_{k}$ and $D_{k}$ denote the mobility and the diffusion coefficient of species $k$, which are related by the Einstein-Smoluchowski relation $D_{k} / \\xi_{k} = k_{\\mathrm{B}}T$ [11,12].\n",
+    "import scipy.optimize\n",
     "\n",
-    "Finally, the total number density flux combining effects of diffusion and advection reads\n",
+    "import matplotlib as mpl\n",
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib.animation as animation\n",
+    "import matplotlib.quiver\n",
     "\n",
-    "\\begin{equation}\n",
-    "\\vec{j}_{k} = \\vec{j}_{k}^{\\mathrm{diff}} + \\vec{j}_{k}^{\\mathrm{adv.}} = -D_{k} \\nabla c_{k} - \\xi_{k} z_{k} e c_{k} \\nabla \\Phi + c_{k} \\vec{u} , \n",
-    "\\end{equation}\n",
+    "import tempfile\n",
+    "import base64\n",
     "\n",
-    "where the first term represents Fick's law of diffusion in the absence of an external potential, the second term gives the additional flux due to an external (in this case electrostatic) potential, and the last term introduces the influence of the motion of the underlying solvent."
+    "plt.rcParams.update({'font.size': 16})"
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
+   "outputs": [],
    "source": [
-    "#### Electrostatics\n",
-    "\n",
-    "The dynamics of the charged species in a typical micro- or nanofluidic system are slow compared to the relaxation of the electromagnetic fields.\n",
-    "This allows us to use stationary equations to model electromagnetic effects.\n",
-    "We further assume that the modeled species do not carry permanent magnetic dipoles and that electric currents in the system are small.\n",
-    "Under these conditions, the full set of Maxwell's equations reduces to the Poisson equation\n",
-    "\n",
-    "\\begin{equation}\n",
-    "\\nabla^2 \\Phi = - \\frac{1}{\\varepsilon} \\sum_{k} z_{k} e c_{k} = -4 \\pi l_\\mathrm{B} k_{\\mathrm{B}}T \\sum_{k} z_{k} c_{k} . \n",
-    "\\end{equation}\n",
-    "\n",
-    "Here $\\varepsilon = \\varepsilon_{0} \\varepsilon_r$ denotes the product of the vacuum permittivity $\\varepsilon_{0}$ with the relative permittivity of the solvent $\\varepsilon_r$.\n",
-    "We have also used the Bjerrum-length\n",
-    "\n",
-    "\\begin{equation}\n",
-    "l_\\mathrm{B} = \\frac{e^{2}}{4 \\pi \\varepsilon k_{\\mathrm{B}}T}.\n",
-    "\\end{equation}\n",
-    "\n",
-    "Finally, we have assumed that the permittivity is spatially homogeneous, since this will allow us to use efficient spectral methods to solve this equation.\n"
+    "BOX_L = [80, 80, 1] \n",
+    "AGRID = 1.0\n",
+    "DIFFUSION_COEFFICIENT = 0.06\n",
+    "TAU = 1.0\n",
+    "EXT_FORCE_DENSITY = [0, 0, 0]\n",
+    "FLUID_DENSITY = 1.0\n",
+    "FLUID_VISCOSITY = 1.\n",
+    "FLUID_VELOCITY = [0.04, 0.04, 0.0]\n",
+    "KT = 1.0\n",
+    "\n",
+    "RUN_TIME = 400"
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
+   "outputs": [],
    "source": [
-    "#### Hydrodynamics\n",
-    "\n",
-    "As said before, since the ionic species' contribute at most a few percent to the overall mass, we can safely approximate the overall fluid's mass by the mass of the solvent (typically water) and model the solvents velocity field $\\vec{u}(\\vec{r}, t)$ using the Navier-Stokes equations for an isotropic, incompressible Newtonian fluid\n",
-    "\n",
-    "\\begin{aligned}\n",
-    "\\rho \\big( \\partial_t \\vec{u} + \\left(\\vec{u} \\cdot \\nabla \\right) \\vec{u} \\big) &= -\\nabla p_H + \\eta \\nabla^{2} \\vec{u} + \\vec{f} ,\\\\\n",
-    "\\nabla \\cdot \\vec u &= 0 .\n",
-    "\\end{aligned}\n",
-    "\n",
-    "where $p_H$ denotes the hydrostatic pressure, $\\eta$ the shear viscosity, $\\rho$ the density of the fluid, and $\\vec{f}$ an external body force density.\n",
-    "For the assumption of incompressibility to hold, the Mach number needs to be small – a condition that is fulfilled for micro- and nanofluidic systems with flow velocities on the order of μm/s.\n",
-    "\n",
-    "Earlier we assumed that the ions' velocity relative to the fluid instantaneously relaxes to a stationary state and that this stationary velocity is given by the product of their mobility and the force exerted on them.\n",
-    "For this state to be stationary, all the momentum transferred into the ions by the external force needs to be dissipated into the fluid immediately.\n",
-    "From this we can conclude that the force density acting on the fluid must read\n",
-    "\n",
-    "\\begin{equation}\n",
-    "\\vec{f} = \\sum_{k} \\vec{j}^\\mathrm{diff}_k / \\xi_{k} = -\\sum_{k} (k_\\mathrm{B}T \\nabla c_{k} + z_{k} e c_{k} \\nabla \\Phi) .\n",
-    "\\end{equation}\n",
-    "\n",
-    "Summarizing, the set of electrokinetic equations we solve is given by\n",
-    "\n",
-    "\\begin{aligned}\n",
-    "\\vec{j}_{k} &= -D_{k} \\nabla c_{k} - \\xi_{k} z_{k} e c_{k} \\nabla \\Phi + c_{k} \\vec{u} ,\\\\\n",
-    "\\partial_{t} c_{k} &= -\\nabla \\cdot\\vec{j}_{k} ,\\\\\n",
-    "\\nabla^2 \\Phi &= -4 \\pi l_\\mathrm{B} k_\\mathrm{B}T \\textstyle\\sum_{k} z_{k} c_{k} ,\\\\\n",
-    "\\rho \\big( \\partial_t \\vec{u} + (\\vec{u} \\cdot \\nabla ) \\vec{u} \\big) &= -\\nabla p_H + \\eta \\nabla^{2} \\vec{u} - \\textstyle\\sum_{k} (k_\\mathrm{B}T \\nabla c_{k} + z_{k} e c_{k} \\nabla \\Phi) ,\\\\\n",
-    "\\nabla \\cdot \\vec{u} &= 0 .\n",
-    "\\end{aligned}"
+    "system = espressomd.System(box_l=BOX_L)\n",
+    "system.time_step = TAU\n",
+    "system.cell_system.skin = 0.4"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### EOF in the Slit Pore Geometry\n",
-    "\n",
-    "The slit pore system depicted in Fig. 1 consists of two like charged parallel plates of infinite extent, confining a solution of water and the plates' counterions.\n",
-    "\n",
-    "<figure>\n",
-    "    <img src='figures/schlitzpore_3d.png' alt='missing' style=\"width: 500px;\"/>\n",
-    "    <center>\n",
-    "    <figcaption>Figure 1: Slit pore system and coordinate system used for the analytical calculations.</figcaption>\n",
-    "    </center>\n",
-    "</figure>\n",
-    "\n",
-    "Due to the net neutrality of the system and the translational symmetry in directions parallel to the plates, the potential outside the two plates must be constant.\n",
-    "This means that using periodic or non-periodic boundary conditions makes no difference.\n",
-    "As the system is in equilibrium in the normal direction, the electrokinetic equations for this dimension reduce to the Poisson-Boltzmann equation for the electrostatic potential, which reads\n",
-    "\\begin{equation}\n",
-    "\\partial_x^2 \\Phi(x) = -4 \\pi \\, k_\\mathrm{B}T \\, l_\\mathrm{B} \\, ze \\, c_0 \\cdot \\exp{\\left(-\\frac{ze\\Phi(x)}{k_\\mathrm{B}T}\\right)} \\; ,\n",
-    "\\end{equation}\n",
-    "where $x$ denotes the direction normal to the plates.\n",
-    "The constant $c_0$ has to be chosen such that charge neutrality is fulfilled.\n",
-    "Multiplying by $2 \\partial_x \\Phi(x)$ and applying the inverse chain rule further reduces the equation to first order.\n",
-    "Subsequent separation of variables yields the solution\n",
-    "\\begin{equation}\n",
-    "\\Phi(x) = -\\frac{k_\\mathrm{B}T}{ze} \\cdot \\log \\left[ \\frac{C^2}{8 \\pi \\, k_\\mathrm{B}T \\, l_\\mathrm{B}} \\cdot \\cos^{-2}\\left( \\frac{zeC}{2 k_\\mathrm{B}T} \\cdot x\\right) \\right], \\quad \\left| \\frac{zeC}{2 k_\\mathrm{B}T} \\cdot x \\right| < \\frac \\pi 2\\; .\n",
-    "\\end{equation}\n",
-    "Refer to [4] for details on this calculation.\n",
-    "Knowing that the counterion density $c$ resembles a Boltzmann distribution in the potential $ze \\Phi$ leads to the expression\n",
-    "\\begin{equation}\n",
-    "c(x) = \\frac{C^2}{8 \\pi \\, k_\\mathrm{B}T \\, l_\\mathrm{B}} \\cdot \\cos^{-2} \\left( \\frac{zeC}{2 k_\\mathrm{B}T} \\cdot x \\right) \\; .\n",
-    "\\end{equation}\n",
-    "The constant $C$ is determined by fixing the number of counterions or requiring the E-field to vanish outside the volume contained by the plates.\n",
-    "Both yields\n",
-    "\\begin{equation}\n",
-    "C \\cdot \\tan \\left( \\frac{zed}{4 k_\\mathrm{B}T} \\cdot C \\right) = -4 \\pi \\, k_\\mathrm{B}T \\, l_\\mathrm{B} \\sigma \\; ,\n",
-    "\\end{equation}\n",
-    "where $d$ denotes the distance between the plates and $\\sigma$ their (constant) surface charge density.\n",
-    "\n",
-    "Applying an electric field along one of the directions parallel to the plates does not influence the charge distribution in the normal direction, which allows us to write down the hydrodynamic equations for the parallel direction.\n",
-    "After eliminating all terms from the Navier-Stokes Equations, which vanish due to symmetry, we are left with\n",
-    "\\begin{equation}\n",
-    "\\frac{\\partial_x^2 v_y(x)}{\\partial x^2} = -\\frac{q E C^2}{8 \\, k_\\mathrm{B}T \\, l_\\mathrm{B} \\, \\eta} \\cdot \\cos^{-2} \\left( \\frac{qC}{2 k_\\mathrm{B}T} \\cdot x \\right) \\; ,\n",
-    "\\end{equation}\n",
-    "which yields, by means of simple integration and the application of no-slip boundary conditions\n",
-    "\\begin{equation}\n",
-    "v_y(x) = \\frac{E}{2 \\pi \\, l_\\mathrm{B} \\, \\eta \\, ze} \\cdot \\log \\left[ \\frac{\\cos \\left( \\frac{zeC}{2 k_\\mathrm{B}T} \\cdot x \\right)}{\\cos \\left( \\frac{zeC}{2 k_\\mathrm{B}T} \\cdot \\frac d 2 \\right)} \\right] \\; .\n",
-    "\\end{equation}\n",
-    "\n",
-    "With this tutorial comes a Python script <tt>eof_analytical.py</tt>, which evaluates all these expressions on the same grid as is used in the upcoming simulation."
+    "We use a lattice Boltzmann flow field with constant velocity for advection.\n",
+    "Note that we have to set ``kT=0.0`` here to avoid random fluctuations in the flow velocity."
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
+   "outputs": [],
    "source": [
-    "## Simulation using ESPResSo"
+    "lattice = espressomd.lb.LatticeWalberla(n_ghost_layers=1, agrid=AGRID)\n",
+    "lbf = espressomd.lb.LBFluidWalberla(\n",
+    "                lattice=lattice, density=FLUID_DENSITY, kinematic_viscosity=FLUID_VISCOSITY,\n",
+    "                tau=TAU, ext_force_density=EXT_FORCE_DENSITY, kT=0.0, seed=42)\n",
+    "lbf[:, :, :].velocity = FLUID_VELOCITY\n",
+    "system.lb = lbf"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Mapping SI and Simulation Units\n",
-    "\n",
-    "ESPResSo does not predefine any unit system.\n",
-    "This makes it more flexible but also requires us to spend some thought on the conversion from SI units to simulation units and back.\n",
-    "Since most first time users have trouble with this, we will go through that process in detail here.\n",
-    "\t\n",
-    "Important to note is that ESPResSo's unit system is nothing more than a rescaled variant of the SI system.\n",
-    "All physical formulas you are used to in the SI system remain valid and you can use them to find relations between your units.\n",
-    "Let's start by choosing a unit of length.\n",
-    "Since we are going to deal with Debye layers with extensions of nanometers, a sensible choice is\n",
-    "\n",
-    "\\begin{equation}\n",
-    "[x]=1\\mathrm{nm}.\n",
-    "\\end{equation}\n",
-    "\n",
-    "The involved energies are of the magnitude of $k_{\\mathrm{B}}T$.\n",
-    "We will simulate our system at room temperature ($300\\mathrm{K}$), hence we use as unit of energy\n",
-    "\\begin{equation}\n",
-    "[E]=k_\\mathrm{B}\\cdot 300\\mathrm{K}\\approx 4.14 \\cdot 10^{-21}\\mathrm{J}.\n",
-    "\\end{equation}\n",
-    "\n",
-    "By default ESPResSo has no concept for particle masses (but the feature can be activated).\n",
-    "That means all particle masses are assumed to be $1\\,[\\mathrm{m}]$, which forces us to use the particle mass as mass unit.\n",
-    "For this simulation we use the mass of sodium ions, which is\n",
-    "\\begin{equation}\n",
-    "[m]=23\\mathrm{u}\\approx 3.82\\cdot 10^{-26}\\mathrm{kg}.\n",
-    "\\end{equation}\n",
-    "\n",
-    "For the relation\n",
-    "\\begin{equation}\n",
-    "E=\\frac 1 2 mv^2\n",
-    "\\end{equation}\n",
-    "\n",
-    "to hold, the unit of time has to be defined so that\n",
-    "\\begin{equation}\n",
-    "[E]=[m]\\cdot\\frac{[x]^2}{[t]^2}.\n",
-    "\\end{equation}\n",
-    "\n",
-    "From that we get the missing unit of time\n",
-    "\\begin{equation}\n",
-    "[t]=[x]\\cdot\\sqrt{\\frac{[m]}{[E]}}=1\\mathrm{nm}\\cdot\\sqrt{\\frac{23\\mathrm{u}}{k_B\\cdot 300\\mathrm{K}}}\\approx 3.03760648\\cdot 10^{-12}\\mathrm{s}\\approx 3.04\\mathrm{ps}.\n",
-    "\\end{equation}\n",
-    "\n",
-    "The last unit we need is the one of charge.\n",
-    "We choose it to be the elementary charge\n",
-    "\\begin{equation}\n",
-    "[q]=e\\approx 1.60\\cdot 10^{-19}\\mathrm{C}.\n",
-    "\\end{equation}\n",
-    "\n",
-    "We now have all the units necessary to convert our simulation parameters."
+    "To use the electrokinetics-algorithm in ESPResSo, one needs to create an instance of the `EKContainer`-object and pass it a time step `tau` and Poisson solver `solver`.\n",
+    "Since our species is uncharged, we don't need to solve the electrostatic Poisson equation, so we can use the placeholder-class, which is called `EKNone`."
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
+   "outputs": [],
    "source": [
-    "|parameter |value (SI units) | value (simulation units)|\n",
-    "|:---------|----------------:|------------------------:|\n",
-    "|channel width $d$ | $50\\mathrm{nm}$ | $50\\mathrm{[x]}$|\n",
-    "|thermal energy $k_B T$ | $k_B\\cdot 300\\mathrm{K}$ | $1\\mathrm{[E]}$|\n",
-    "|Bjerrum length $l_B$ | $0.7095\\mathrm{nm}$ | $0.7095\\mathrm{[x]}$|\n",
-    "|counterion charge $q$ | $1e$ | $1\\mathrm{[q]}$|\n",
-    "|counterion diffusion coefficient $D$ | $2.0\\cdot 10^{-9}\\mathrm{m^2/s}$ | $0.006075\\mathrm{[x]^2/[t]}$|\n",
-    "|solvent density $\\rho$ | $1.0\\cdot 10^{3}\\mathrm{kg/m^3}$ | $26.18\\mathrm{[m]/[x]^3}$|\n",
-    "|solvent dynamic viscosity $\\eta$ | $1.0\\cdot 10^{-3}\\mathrm{Pa}\\mathrm{s}$ | $79.53\\mathrm{[m]/([x][t])}$|\n",
-    "|external electric field $E$ | $2.585\\cdot 10^{6}\\mathrm{V/m}$ | $0.1\\mathrm{[E]/([q][x])}$|\n"
+    "eksolver = espressomd.electrokinetics.EKNone(lattice=lattice)\n",
+    "system.ekcontainer = espressomd.electrokinetics.EKContainer(tau=TAU, solver=eksolver)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "ESPResSo determines the strength of the electrostatic interactions via the Bjerrum-length $l_\\mathrm{B}$.\n",
-    "That is the length for which the electrostatic interaction energy of two elementary charges equals the thermal energy\n",
-    "\n",
-    "\\begin{equation}\n",
-    "k_\\mathrm{B} T=\\frac{e^2}{4\\pi\\varepsilon_0\\varepsilon_r}\\cdot\\frac 1 {l_\\mathrm{B}}.\n",
-    "\\end{equation}\n",
-    "\n",
-    "This yields for water at $300K$ with $\\varepsilon_r = 78.54$, a Bjerrum length of $l_\\mathrm{B}\\approx 0.7095\\mathrm{nm}$."
+    "Now, we can add diffusive species to the container to integrate their dynamics."
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "solution2": "hidden",
+    "solution2_first": true
+   },
    "source": [
-    "### Setting up the slit pore system\n",
-    "\n",
-    "The script for this simulation comes with this tutorial and is called <tt>eof_electrokinetics.py</tt>.\n",
-    "All used commands are documented in the User's Guide in the section called **Electrokinetics**.\n",
+    "# Exercise:\n",
+    "- Create an instance of the [`espressomd.electrokinetics.EKSpecies`]() and add it to the system with [`system.ekcontainer.add()`](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKContainer.add). \n",
     "\n",
-    "We first set up a periodic simulation box of the desired dimensions.\n",
-    "Note that the dimensions are, of course, given in simulation units."
+    "# Hint:\n",
+    "- Use the variables `DIFFUSION_COEFFICIENT`, `KT` and `TAU` defined above.\n",
+    "- Enable both `advection` and `friction_coupling`.\n",
+    "- Make sure to initialize the `density` with 0.0, and disable electrostatics by setting `valency` to 0.0 as well."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "solution2": "hidden"
+   },
+   "source": [
+    "```python\n",
+    "species = espressomd.electrokinetics.EKSpecies(\n",
+    "    lattice=lattice, density=0.0, kT=KT,\n",
+    "    diffusion=DIFFUSION_COEFFICIENT, valency=0.0,\n",
+    "    advection=True, friction_coupling=True,\n",
+    "    ext_efield=[0., 0., 0.], tau=TAU)\n",
+    "system.ekcontainer.add(species)\n",
+    "```"
    ]
   },
   {
@@ -351,38 +205,13 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "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",
-    "espressomd.assert_features([\"WALBERLA\", \"WALBERLA_FFT\"])\n",
-    "\n",
-    "import tqdm\n",
-    "import numpy as np\n",
-    "%matplotlib inline\n",
-    "import matplotlib.pyplot as plt\n",
-    "plt.rcParams.update({'font.size': 16})\n",
-    "\n",
-    "box_y = 6\n",
-    "box_z = 6\n",
-    "width = 50\n",
-    "\n",
-    "padding = 1\n",
-    "box_x = width + 2 * padding\n",
-    "\n",
-    "system = espressomd.System(box_l=[box_x, box_y, box_z])"
-   ]
+   "source": []
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We then store all the parameters we calculated earlier.\n",
-    "At this point, these parameters only reside in Python variables.\n",
-    "They will only be used by ESPResSo once they are being passed to the respective initialization functions."
+    "To compare our simulation to the fundamental solution of the advection-diffusion equation, we need to approximate a delta-droplet, which can be achieved by having a non-zero density only at the center of the domain."
    ]
   },
   {
@@ -391,28 +220,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Set the electrokinetic parameters\n",
-    "\n",
-    "agrid = 1.0\n",
-    "dt = 0.5\n",
-    "kT = 4.0\n",
-    "bjerrum_length = 0.7095\n",
-    "permittivity = 1. / (4 * np.pi * bjerrum_length)\n",
-    "D = 0.006075\n",
-    "valency = 1.0\n",
-    "viscosity_dynamic = 79.53\n",
-    "density_water = 26.15\n",
-    "sigma = -0.05\n",
-    "ext_force_density = [0.0, 0.1, 0.0]\n",
-    "\n",
-    "single_precision = False"
+    "system.ekcontainer[0][BOX_L[0] // 2, BOX_L[1] // 2, 0].density = 1.0"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Before we initialize the actual electrokinetics algorithm, we need to set the time step and some other parameters that are not actually used, but would otherwise lead to error messages."
+    "Now everything is set and we can finally run the simulation by running the integrator."
    ]
   },
   {
@@ -421,19 +236,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Set the simulation parameters\n",
-    "\n",
-    "system.time_step = dt\n",
-    "system.cell_system.skin = 0.2\n",
-    "system.thermostat.turn_off()\n",
-    "integration_length = 600"
+    "system.integrator.run(RUN_TIME)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We can now set up the electrokinetics algorithm."
+    "For comparison, we prepare the analytical solution and show the 2D-density as well as a slice through the center of the droplet."
    ]
   },
   {
@@ -442,7 +252,20 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "lattice = espressomd.lb.LatticeWalberla(agrid=agrid, n_ghost_layers=1)"
+    "def calc_gaussian(pos: np.ndarray, time: int, D: float):\n",
+    "    dim = pos.shape[-1]\n",
+    "    return (4 * np.pi * D * time)**(-dim / 2) * np.exp(-np.sum(np.square(pos), axis=-1) / (4 * D * time))\n",
+    "\n",
+    "pos = np.zeros((*BOX_L[:2], 2))\n",
+    "xpos = np.arange(-BOX_L[0] // 2, BOX_L[0] // 2)\n",
+    "ypos = np.arange(-BOX_L[1] // 2, BOX_L[1] // 2)\n",
+    "\n",
+    "pos[..., 1], pos[..., 0] = np.meshgrid(xpos, ypos)\n",
+    "\n",
+    "# add the advection shift\n",
+    "pos -= np.asarray(FLUID_VELOCITY[:2]) * RUN_TIME * TAU\n",
+    "\n",
+    "analytic_density = calc_gaussian(pos=pos, time=RUN_TIME * TAU, D=DIFFUSION_COEFFICIENT)"
    ]
   },
   {
@@ -451,11 +274,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "viscosity_kinematic = viscosity_dynamic / density_water\n",
-    "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"
+    "fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(15, 7))\n",
+    "\n",
+    "ax1.imshow(system.ekcontainer[0][:, :, 0].density, origin=\"lower\", vmin=0, vmax=6e-3)\n",
+    "ax1.set_title(\"simulation\")\n",
+    "\n",
+    "imshow = ax2.imshow(analytic_density, origin=\"lower\", vmin=0, vmax=6e-3)\n",
+    "ax2.set_title(\"analytic\")\n",
+    "fig.colorbar(imshow, ax=[ax1, ax2], shrink=0.8)\n",
+    "plt.show()"
    ]
   },
   {
@@ -464,70 +291,120 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "eksolver = espressomd.electrokinetics.EKFFT(lattice=lattice,\n",
-    "                                            permittivity=permittivity,\n",
-    "                                            single_precision=single_precision)\n",
+    "values_diagonal = np.diagonal(system.ekcontainer[0][:, :, 0].density)\n",
+    "analytic_diagonal = np.diagonal(analytic_density)\n",
+    "positions_diagonal = np.arange(len(values_diagonal)) * np.sqrt(2) * AGRID\n",
+    "\n",
+    "def gaussian_kernel(x, magnitude, mu, sigma):\n",
+    "    return magnitude * np.exp(-(x - mu)**2 / (2 * sigma**2))\n",
+    "\n",
+    "popt, pcov = scipy.optimize.curve_fit(gaussian_kernel, positions_diagonal, analytic_diagonal, p0=[0.05,70,5.])\n",
+    "positions_analytic = np.concatenate([[positions_diagonal[0]],\n",
+    "                                     np.linspace(popt[1] - 5 * popt[2], popt[1] + 5 * popt[2], 120),\n",
+    "                                    [positions_diagonal[-1]]])\n",
+    "values_analytic = gaussian_kernel(positions_analytic, *popt)\n",
+    "\n",
+    "fig = plt.figure(figsize=(8, 5))\n",
+    "ax = fig.gca()\n",
+    "ax.plot(positions_diagonal, values_diagonal, \"o\", mfc=\"none\", label=\"simulation\")\n",
+    "ax.plot(positions_analytic, values_analytic, label=\"analytic\")\n",
+    "\n",
+    "ax.set_xlabel(\"position\")\n",
+    "ax.set_ylabel(\"density\")\n",
     "\n",
-    "system.ekcontainer = espressomd.electrokinetics.EKContainer(solver=eksolver, tau=dt)"
+    "plt.legend()\n",
+    "plt.show()"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "density_counterions = -2.0 * sigma / width\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)"
+    "From the plot one can see that the position of the density-peak matches well. However, one also sees that the droplet in the simulation has spread more than it should. The reason is that the discretization used for the advection term introduces an artifical, additional diffusion to the system. This is a fundamental limitation of the algorithm, which is why it cannot be applied to pure advection problems."
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "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)"
+    "# 3. Electroosmotic flow"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "wall_left = espressomd.shapes.Wall(normal=[1, 0, 0], dist=padding)\n",
-    "wall_right = espressomd.shapes.Wall(normal=[-1, 0, 0], dist=-(padding + width))"
+    "The next system in this tutorial is a simple slit pore, as shown in Figure 1. It consists of an infinite plate capacitor with an electrolytic solution trapped in between the plates. The plates of the capactior carry a constant surface charge and the counterions are solvated in the liquid. \n",
+    "\n",
+    "Charge attraction will cause the ions to accumulate near the surfaces, forming a characteristic ion density profile, which can be calculated analytically using the Poisson-Boltzmann equation. Since the system has translational symmetry in the directions parallel to the plates, the equations for parallel and orthogonal direction decouple. This means that applying an external electric field in a direction parallel to the plates will not change the distribution of the ions along the orthogonal direction. It will however cause motion of the ions and consequently the fluid: The characteristic flow profile of electroosmotic flow.\n",
+    "\n",
+    "<a id=\"fig_1\"></a><figure>\n",
+    "<img src=\"figures/schlitzpore_3d.png\" style=\"float: center; width: 60%\">\n",
+    "<center>\n",
+    "<figcaption>Figure 1: Slit pore system and coordinate system used for the analytical calculations.</figcaption>\n",
+    "</center>\n",
+    "</figure>"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "ekspecies[padding:-padding, :, :].density = density_counterions"
+    "### Analytical solution\n",
+    "\n",
+    "Due to the symmetries of the system, it effectively reduces to a 1D problem along the orthogonal axis. The system can be described by the Poisson-Boltzmann equation:\n",
+    "\n",
+    "$$\n",
+    "\\partial_{x}^2 \\phi(x) = \\frac{1}{\\varepsilon_{0} \\varepsilon_{\\mathrm{r}}} \\sum_{i} z_{i} e n_{i}(x) \\exp \\left( -\\frac{z_{i} e \\phi(x)}{k_{\\mathrm{B}} T} \\right)\n",
+    "$$\n",
+    "\n",
+    "where $x$ is the normal-direction of the plates. Since we will only simulate a single ion species, the counterions, the sum only has a single summand. The solution for the potential is then given by:\n",
+    "\n",
+    "$$\n",
+    "\\phi(x) = -\\frac{k_{B}T}{z e} \\log \\left[ \\frac{C^2 \\varepsilon_{0} \\varepsilon_{\\mathrm{r}}}{2 k_{B}T } \\cos^{-2} \\left( \\frac{z e C}{2 k_{B} T} x \\right) \\right], \\qquad \\text{with } \\left\\| \\frac{z e C}{2 k_{B} T} \\right\\| < \\frac{\\pi}{2},\n",
+    "$$\n",
+    "\n",
+    "where $C$ is an integration constant that is to be determined by the boundary conditions. The ion density follows then from the potential as\n",
+    "\n",
+    "$$\n",
+    "n(x) = \\frac{C^2 \\varepsilon_{0} \\varepsilon_{\\mathrm{r}}}{2 k_{B}T} \\cos^{-2} \\left( \\frac{z e C}{2 k_{B} T} x \\right).\n",
+    "$$\n",
+    "\n",
+    "To find the integration constant we use fact that the total system has to be charge neutral, i.e., the total charge on the plates is counterbalanced by the counterions. This leads to the following equation\n",
+    "\n",
+    "$$\n",
+    "C \\tan \\left( \\frac{z e d}{4 k_{B} T} C \\right) = - \\frac{e^2}{\\varepsilon_{0} \\varepsilon_{\\mathrm{r}}} \\sigma,\n",
+    "$$\n",
+    "\n",
+    "where $\\sigma$ is the surface charge density of the plates. This is a transcendental equation, which must be solved numerically to find $C$. \n",
+    "\n",
+    "The electric field is applied in the $y$-direction, parallel to the plates.\n",
+    "Fluid flow is described by the incompressible Navier-Stokes equation, which due to the symmetries of the system reduces to the one-dimensional problem\n",
+    "\n",
+    "$$\n",
+    "\\frac{\\partial^2 v_{y}(x)}{\\partial x^2} = - \\frac{\\varepsilon_{0} \\varepsilon_{\\mathrm{r}} z e E C^2}{2 k_{B}T \\eta} \\cos^{-2}\\left( \\frac{q C}{2 k_{B} T} x \\right).\n",
+    "$$\n",
+    "\n",
+    "This equation can be solved analytically and the solution is given by\n",
+    "\n",
+    "$$\n",
+    "v_{y}(x) = \\frac{2 \\varepsilon_{0} \\varepsilon_{\\mathrm{r}} k_{B} T E}{\\eta z e} \\log \\left( \\frac{\\cos \\left( \\displaystyle\\frac{z e C}{2 k_{B} T} x \\right)}{\\cos \\left( \\displaystyle\\frac{z e C}{2 k_{B} T} \\frac{d}{2} \\right)} \\right),\n",
+    "$$\n",
+    "\n",
+    "where $d$ denotes the distance between the two plates. Finally, the shear-stress of this problem is given by\n",
+    "\n",
+    "$$\n",
+    "\\sigma(x) = \\mu \\frac{\\partial v_{y}(x)}{\\partial x}\n",
+    "$$"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "ekspecies[:padding, :, :].density = 0.0\n",
-    "ekspecies[-padding:, :, :].density = 0.0\n",
+    "### Numerical solution\n",
     "\n",
-    "ekwallcharge[:padding, :, :].density = -sigma / valency / padding\n",
-    "ekwallcharge[-padding:, :, :].density = -sigma / valency / padding"
+    "We start by resetting the system and defining the necessary parameters."
    ]
   },
   {
@@ -536,12 +413,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "for shape_obj in (wall_left, wall_right):\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.])"
+    "system.ekcontainer = None\n",
+    "system.lb = None"
    ]
   },
   {
@@ -550,101 +423,44 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Integrate the system\n",
-    "for i in tqdm.trange(100):\n",
-    "    system.integrator.run(integration_length)\n",
-    "\n",
-    "# Output\n",
-    "position_list = []\n",
-    "density_list = []\n",
-    "velocity_list = []\n",
-    "pressure_xy_list = []\n",
-    "\n",
-    "for i in range(int(box_x / agrid)):\n",
-    "    if (i * agrid >= padding) and (i * agrid < box_x - padding):\n",
-    "        position = i * agrid - padding - width / 2.0 + agrid / 2.0\n",
-    "        position_list.append(position)\n",
-    "        \n",
-    "        node_idxs = (i, int(box_y / (2 * agrid)), int(box_z / (2 * agrid)))\n",
+    "AGRID = 1.0\n",
+    "TAU = 1.0\n",
+    "KT = 2.0\n",
+    "PERMITTIVITY = 0.28\n",
+    "DIFFUSION_COEFFICIENT = 0.25\n",
+    "VALENCY = 1.0\n",
+    "VISCOSITY_DYNAMIC = 0.5\n",
+    "DENSITY_FLUID = 1.0\n",
+    "SURFACE_CHARGE_DENSITY = -0.05\n",
+    "EXT_FORCE_DENSITY = [0.0, 0.01, 0.0]\n",
+    "\n",
+    "SINGLE_PRECISION = False\n",
     "\n",
-    "        # density\n",
-    "        density_list.append(ekspecies[node_idxs].density)\n",
-    "\n",
-    "        # velocity\n",
-    "        velocity_list.append(lbf[node_idxs].velocity[1])\n",
+    "padding = 1\n",
+    "WIDTH = 126\n",
+    "BOX_L = [(WIDTH + 2 * padding) * AGRID, 1, 1]\n",
     "\n",
-    "        # xz component pressure tensor\n",
-    "        pressure_xy_list.append(lbf[node_idxs].pressure_tensor[0, 1])\n",
+    "system.cell_system.skin = 0.4\n",
+    "system.box_l = BOX_L\n",
+    "system.time_step = TAU\n",
     "\n",
-    "np.savetxt(\"eof_simulation.dat\",\n",
-    "           np.column_stack((position_list,\n",
-    "                            density_list,\n",
-    "                            velocity_list,\n",
-    "                            pressure_xy_list)),\n",
-    "           header=\"#position calculated_density calculated_velocity calculated_pressure_xy\")"
+    "RUN_TIME = 200"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": true
-   },
-   "outputs": [],
+   "cell_type": "markdown",
+   "metadata": {},
    "source": [
-    "from scripts import eof_analytical # executes automatically upon import\n",
-    "\n",
-    "# read analytical solution and simulation data\n",
-    "data_an = np.loadtxt(\"eof_analytical.dat\")\n",
-    "data_ek = np.loadtxt(\"eof_simulation.dat\")\n",
-    "\n",
-    "fig1 = plt.figure(figsize=(16, 4.5))\n",
-    "ax = fig1.add_subplot(131)\n",
-    "ax.plot(data_an[:, 0], data_an[:, 1], label=\"analytical\")\n",
-    "ax.plot(data_ek[:, 0], data_ek[:, 1], \"o\", mfc=\"none\", label=\"simulation\")\n",
-    "ax.set_xlabel(\"x-position\")\n",
-    "ax.set_ylabel(\"Counter-ion density\")\n",
-    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
-    "ax.legend(loc=\"best\")\n",
-    "\n",
-    "ax = fig1.add_subplot(132)\n",
-    "ax.plot(data_an[:, 0], data_an[:, 2], label=\"analytical\")\n",
-    "ax.plot(data_ek[:, 0], data_ek[:, 2], \"o\", mfc=\"none\", label=\"simulation\")\n",
-    "ax.set_xlabel(\"x-position\")\n",
-    "ax.set_ylabel(\"Fluid velocity\")\n",
-    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
-    "ax.legend(loc=\"best\")\n",
-    "\n",
-    "ax = fig1.add_subplot(133)\n",
-    "ax.plot(data_an[:, 0], data_an[:, 3], label=\"analytical\")\n",
-    "ax.plot(data_ek[:, 0], data_ek[:, 3], \"o\", mfc=\"none\", label=\"simulation\")\n",
-    "ax.set_xlabel(\"x-position\")\n",
-    "ax.set_ylabel(\"Fluid shear stress xz\")\n",
-    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
-    "ax.legend(loc=\"best\")\n",
-    "\n",
-    "plt.tight_layout()\n",
-    "plt.show()"
+    "We can now set up the electrokinetics algorithm as in the first part of the tutorial, starting with the LB-method."
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
+   "outputs": [],
    "source": [
-    "## References\n",
-    "\n",
-    "[1] F. Capuani, I. Pagonabarraga and D. Frenkel *Discrete solution of the electrokinetic equations* The Journal of Chemical Physics, 2004  \n",
-    "[2] G. Rempfer *A Lattice based Model for Electrokinetics* Master's thesis, University of Stuttgart, 2013  \n",
-    "[3] G. Rempfer, G. B. Davies, C. Holm and J. de Graaf *Reducing spurious flow in simulations of electrokinetic phenomena* The Journal of Chemical Physics, 2016  \n",
-    "[4] G. Rempfer *Lattice-Boltzmann simulations in complex geometries* Bachelor's thesis, University of Stuttgart, Institute for Computational Physics, 2010  \n",
-    "[5] M. Deserno and C. Holm and S. May, *Fraction of Condensed Counterions around a Charged Rod: Comparison of Poisson-Boltzmann Theory and Computer Simulations* Macromolecules, 2000  \n",
-    "[6] C. Holm, P. K&eacute;kicheff and R. Podgornik *Electrostatic Effects in Soft Matter and Biophysics* Kluwer Academic Publishers, 2001  \n",
-    "[7] M. Deserno and C. Holm *Cell-model and Poisson-Boltzmann-theory: A brief introduction* Electrostatic Effects in Soft Matter and Biophysics, Kluwer Academic Publishers, 2001  \n",
-    "[8] J de Graaf., G. Rempfer and C. Holm *Diffusiophoretic Self-Propulsion for Partially Catalytic Spherical Colloids* IEEE T. Nanobiosci., 2014  \n",
-    "[9] M. Deserno *Counterion condensation for rigid linear polyelectrolytes* Universit&auml;t Mainz, 2000  \n",
-    "[10] J. de Graaf, N Boon, M Dijkstra and R. van Roij *Electrostatic interactions between Janus particles* The Journal of Chemical Physics, 2012   \n",
-    "[11] A. Einstein *&Uuml;ber die von der molekularkinetischen Theorie der W&auml;rme geforderte Bewegung von in ruhenden Fl&uuml;ssigkeiten suspendierten Teilchen* Annalen der Physik, 1905  \n",
-    "[12] M. von Smoluchowski *Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen* Annalen der Physik, 1906  \n"
+    "lattice = espressomd.lb.LatticeWalberla(agrid=AGRID, n_ghost_layers=1)"
    ]
   },
   {
@@ -652,12 +468,816 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "viscosity_kinematic = VISCOSITY_DYNAMIC / DENSITY_FLUID\n",
+    "lbf = espressomd.lb.LBFluidWalberla(lattice=lattice, density=DENSITY_FLUID,\n",
+    "                                    kinematic_viscosity=viscosity_kinematic,\n",
+    "                                    tau=TAU, single_precision=SINGLE_PRECISION)\n",
+    "system.lb = lbf"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Since our species are going to carry a charge now, we need to solve the full electrostatic problem. For that, we have to specify an actual solver."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "solution2": "hidden",
+    "solution2_first": true
+   },
+   "source": [
+    "# Exercise: \n",
+    "- Set up a Poisson solver for the electrostatic interaction and use it to create an instance of the [EKContainer](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKContainer) \n",
+    "- Attach the container to the [`system.ekcontainer`](https://espressomd.github.io/doc/espressomd.html#espressomd.system.System.ekcontainer).\n",
+    "\n",
+    "# Hint:\n",
+    "- Use an [EKFFT](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKFFT)-object as the Poisson-solver."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "solution2": "hidden"
+   },
+   "source": [
+    "```python\n",
+    "eksolver = espressomd.electrokinetics.EKFFT(lattice=lattice, permittivity=PERMITTIVITY,\n",
+    "                                            single_precision=SINGLE_PRECISION)\n",
+    "system.ekcontainer = espressomd.electrokinetics.EKContainer(tau=TAU, solver=eksolver)\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To simulate the system, we will use two different ion species: The counterions are propagated in the fluid. The second species will be used to describe the surface charge on the plates and therefore has to be stationary (i.e. no advection, no diffusion)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ekspecies = espressomd.electrokinetics.EKSpecies(lattice=lattice, density=0.0, kT=KT, diffusion=DIFFUSION_COEFFICIENT, valency=VALENCY, advection=True, friction_coupling=True, ext_efield=EXT_FORCE_DENSITY, single_precision=SINGLE_PRECISION, tau=TAU)\n",
+    "system.ekcontainer.add(ekspecies)\n",
+    "\n",
+    "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=TAU)\n",
+    "system.ekcontainer.add(ekwallcharge)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we set the initial conditions for the ion densities. The counterions will be initialized with a homogeneous distribution, excluding the cells used as boundaries. The surface charge density is homogeneously distributed in the boundary cells."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "density_counterions = -2.0 * SURFACE_CHARGE_DENSITY / VALENCY / WIDTH\n",
+    "ekspecies[padding:-padding, :, :].density = density_counterions\n",
+    "\n",
+    "ekspecies[:padding, :, :].density = 0.0\n",
+    "ekspecies[-padding:, :, :].density = 0.0\n",
+    "\n",
+    "ekwallcharge[:padding, :, :].density = -SURFACE_CHARGE_DENSITY / VALENCY / padding\n",
+    "ekwallcharge[-padding:, :, :].density = -SURFACE_CHARGE_DENSITY / VALENCY / padding"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We now have to specify the boundary conditions. For this, we use ESPResSo's`shapes`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "wall_left = espressomd.shapes.Wall(normal=[1, 0, 0], dist=padding)\n",
+    "wall_right = espressomd.shapes.Wall(normal=[-1, 0, 0], dist=-(padding + WIDTH))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "At both of them we specify no-flux and zero-density boundary conditions for the counterions. Furthermore, we set a no-slip boundary condition for the fluid."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "solution2": "hidden",
+    "solution2_first": true
+   },
+   "source": [
+    "# Exercise\n",
+    "At both walls, set\n",
+    "\n",
+    "- No-flux boundary conditions for the counterions (``ekspecies``)\n",
+    "- Zero-density boundary conditions for the counterions\n",
+    "- No-slip boundary conditions for the fluid (``lbf``)\n",
+    "\n",
+    "# Hints\n",
+    "\n",
+    "- Use the shapes defined above and ``add_boundary_from_shape`` for [EK species](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKSpecies.add_boundary_from_shape) and [LB fluids](https://espressomd.github.io/doc/espressomd.html#espressomd.lb.LBFluidWalberla.add_boundary_from_shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "solution2": "hidden"
+   },
+   "source": [
+    "```python\n",
+    "for wall in (wall_left, wall_right):\n",
+    "    ekspecies.add_boundary_from_shape(shape=wall, value=[0., 0., 0.], boundary_type=espressomd.electrokinetics.FluxBoundary)\n",
+    "    ekspecies.add_boundary_from_shape(shape=wall, value=0.0, boundary_type=espressomd.electrokinetics.DensityBoundary)\n",
+    "    lbf.add_boundary_from_shape(shape=wall, velocity=[0., 0., 0.])\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we can finally integrate the system and extract the ion density profile, the fluid velocity profile as well as the pressure-tensor profile."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "for i in tqdm.trange(80):\n",
+    "    system.integrator.run(RUN_TIME)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mid_y = int(system.box_l[1] / (2 * AGRID))\n",
+    "mid_z = int(system.box_l[2] / (2 * AGRID))\n",
+    "density_eof = ekspecies[padding:-padding, mid_y, mid_z].density\n",
+    "velocity_eof = lbf[padding:-padding, mid_y, mid_z].velocity[:, 1]\n",
+    "pressure_tensor_eof = lbf[padding:-padding, mid_y, mid_z].pressure_tensor[:, 0, 1]\n",
+    "\n",
+    "positions = (np.arange(len(density_eof)) - WIDTH / 2 + 0.5) * AGRID"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For comparison, we calculate the analytic solution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def transcendental_equation(c, distance, kT, sigma, valency, permittivity) -> float:\n",
+    "    elementary_charge = 1.0\n",
+    "    return c * np.tan(valency * elementary_charge * distance / (4 * kT) * c) + sigma / permittivity\n",
+    "\n",
+    "solution = scipy.optimize.fsolve(func=transcendental_equation, x0=0.001, args=(WIDTH, KT, SURFACE_CHARGE_DENSITY, VALENCY, PERMITTIVITY))\n",
+    "\n",
+    "\n",
+    "def eof_density(x, c, permittivity, elementary_charge, valency, kT):\n",
+    "    return c**2 * permittivity / (2 * kT) / (np.cos(valency * elementary_charge * c / (2 * kT) * x))**2\n",
+    "\n",
+    "def eof_velocity(x, c, permittivity, elementary_charge, valency, kT, ext_field, distance, viscosity):\n",
+    "    return 2 * kT * ext_field * permittivity / (viscosity * elementary_charge * valency) * np.log(np.cos(valency * elementary_charge * c / (2 * kT) * x) / np.cos(valency * elementary_charge * c / (2 * kT) * distance / 2))\n",
+    "\n",
+    "def eof_pressure_tensor(x, c, elementary_charge, valency, kT, ext_field, permittivity):\n",
+    "    return permittivity * ext_field * c * np.tan(valency * elementary_charge * c / (2 * kT) * x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "analytic_density_eof = eof_density(x=positions, c=solution, permittivity=PERMITTIVITY, elementary_charge=1.0, valency=VALENCY, kT=KT)\n",
+    "analytic_velocity_eof = eof_velocity(x=positions, c=solution, permittivity=PERMITTIVITY, elementary_charge=1.0, valency=VALENCY, kT=KT, ext_field=EXT_FORCE_DENSITY[1], distance=WIDTH, viscosity=VISCOSITY_DYNAMIC)\n",
+    "analytic_pressure_tensor_eof = eof_pressure_tensor(x=positions, c=solution, elementary_charge=1.0, valency=VALENCY, kT=KT, ext_field=EXT_FORCE_DENSITY[1], permittivity=PERMITTIVITY)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig1 = plt.figure(figsize=(16, 4.5))\n",
+    "fig1.suptitle(\"electroosmotic flow\")\n",
+    "\n",
+    "ax = fig1.add_subplot(131)\n",
+    "ax.plot(positions, density_eof, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n",
+    "ax.plot(positions, analytic_density_eof, label=\"analytic\")\n",
+    "ax.set_xlabel(\"x-position\")\n",
+    "ax.set_ylabel(\"Counter-ion density\")\n",
+    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
+    "ax.legend(loc=\"best\")\n",
+    "\n",
+    "ax = fig1.add_subplot(132)\n",
+    "ax.plot(positions, velocity_eof, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n",
+    "ax.plot(positions, analytic_velocity_eof, label=\"analytic\")\n",
+    "ax.set_xlabel(\"x-position\")\n",
+    "ax.set_ylabel(\"Fluid velocity\")\n",
+    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
+    "ax.legend(loc=\"best\")\n",
+    "\n",
+    "ax = fig1.add_subplot(133)\n",
+    "ax.plot(positions, pressure_tensor_eof, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n",
+    "ax.plot(positions, analytic_pressure_tensor_eof, label=\"analytic\")\n",
+    "ax.set_xlabel(\"x-position\")\n",
+    "ax.set_ylabel(\"Fluid shear stress xz\")\n",
+    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
+    "ax.legend(loc=\"best\")\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In the plots one can see that the analytic solution for the electroosmotic flow matches the simulation very well. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Comparison to pressure-driven flow\n",
+    "We can compare electroosmotic flow to pressure-driven flow. For this, we turn off the external electric field and enable a constant external force density on the fluid."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "EXT_FORCE_DENSITY = [0.0, 0.000004, 0.0]\n",
+    "\n",
+    "ekspecies.ext_efield = [0.0, 0.0, 0.0]\n",
+    "lbf.ext_force_density = EXT_FORCE_DENSITY"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "for i in tqdm.trange(70):\n",
+    "    system.integrator.run(RUN_TIME)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "density_pressure = ekspecies[padding:-padding, mid_y, mid_z].density\n",
+    "velocity_pressure = lbf[padding:-padding, mid_y, mid_z].velocity[:, 1]\n",
+    "pressure_tensor_pressure = lbf[padding:-padding, mid_y, mid_z].pressure_tensor[:, 0, 1]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The analytic solution for pressure-driven flow between two infinite parallel plates is known as the Poiseuille flow."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def pressure_velocity(x, distance, ext_field, viscosity):\n",
+    "    return ext_field / (2 * viscosity) * (distance**2 / 4 - x**2)\n",
+    "\n",
+    "def pressure_pressure_tensor(x, ext_field):\n",
+    "    return ext_field * x"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "analytic_velocity_pressure = pressure_velocity(x=positions, distance=WIDTH, ext_field=EXT_FORCE_DENSITY[1], viscosity=VISCOSITY_DYNAMIC)\n",
+    "analytic_pressure_tensor_pressure = pressure_pressure_tensor(x=positions, ext_field=EXT_FORCE_DENSITY[1])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig1 = plt.figure(figsize=(16, 4.5))\n",
+    "fig1.suptitle(\"pressure-driven flow\")\n",
+    "\n",
+    "ax = fig1.add_subplot(131)\n",
+    "ax.plot(positions, density_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n",
+    "ax.plot(positions, analytic_density_eof, label=\"analytic\")\n",
+    "ax.set_xlabel(\"x-position\")\n",
+    "ax.set_ylabel(\"counter-ion density\")\n",
+    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
+    "ax.legend(loc=\"best\")\n",
+    "\n",
+    "ax = fig1.add_subplot(132)\n",
+    "ax.plot(positions, velocity_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n",
+    "ax.plot(positions, analytic_velocity_pressure, label=\"analytic\")\n",
+    "ax.set_xlabel(\"x-position\")\n",
+    "ax.set_ylabel(\"fluid velocity\")\n",
+    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
+    "ax.legend(loc=\"best\")\n",
+    "\n",
+    "ax = fig1.add_subplot(133)\n",
+    "ax.plot(positions, pressure_tensor_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n",
+    "ax.plot(positions, analytic_pressure_tensor_pressure, label=\"analytic\")\n",
+    "ax.set_xlabel(\"x-position\")\n",
+    "ax.set_ylabel(\"fluid shear stress xz\")\n",
+    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
+    "ax.legend(loc=\"best\")\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As one can again see, the body force on the fluid did non alter the ion-density profile.\n",
+    "However, because the force now applies homogeneously on the whole fluid, the flow profile looks parabolic."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To see the difference between the two types of flows, we plot the simulation data together in one plot."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "fig1 = plt.figure(figsize=(16, 4.5))\n",
+    "fig1.suptitle(\"electroosmotic vs. pressure-driven flow comparison\")\n",
+    "\n",
+    "ax = fig1.add_subplot(131)\n",
+    "ax.plot(positions, density_eof, \"o\", mfc=\"none\", ms=4, markevery=0.015, label=\"eof\")\n",
+    "ax.plot(positions, density_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"pressure\")\n",
+    "ax.set_xlabel(\"x-position\")\n",
+    "ax.set_ylabel(\"counter-ion density\")\n",
+    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
+    "ax.legend(loc=\"best\")\n",
+    "\n",
+    "ax = fig1.add_subplot(132)\n",
+    "ax.plot(positions, velocity_eof, \"o\", mfc=\"none\", markevery=0.015, label=\"eof\")\n",
+    "ax.plot(positions, velocity_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"pressure\")\n",
+    "ax.set_xlabel(\"x-position\")\n",
+    "ax.set_ylabel(\"fluid velocity\")\n",
+    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
+    "ax.legend(loc=\"best\")\n",
+    "\n",
+    "ax = fig1.add_subplot(133)\n",
+    "ax.plot(positions, pressure_tensor_eof, \"o\", mfc=\"none\", markevery=0.015, label=\"eof\")\n",
+    "ax.plot(positions, pressure_tensor_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"pressure\")\n",
+    "ax.set_xlabel(\"x-position\")\n",
+    "ax.set_ylabel(\"fluid shear stress xz\")\n",
+    "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n",
+    "ax.legend(loc=\"best\")\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Looking at the fluid velocity plot, one can see that the electroosmotic flow profile flattens significantly faster towards the center of the channel when compared to the pressure driven flow. The reason for this is the accumulation of the counterion-density towards the oppositely charged plates. Here, the driving electric field causes the highest force on the fluid, which decays towards the center of the channel. In contrast, the Poiseuille-flow is driven by a constant, uniform driving force."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 4. Reaction in turbulent flow"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To showcase the reaction feature of our electrokinetics algorithm, we simulate a simple reaction in complex flow.\n",
+    "For this, we choose a geometry of rigid cylinders.\n",
+    "At large flow velocities, a [Kármán vortex street](https://en.wikipedia.org/wiki/K%C3%A1rm%C3%A1n_vortex_street), i.e., a repeating pattern of swirling vorticies behind the obstacle, develops.\n",
+    " \n",
+    "To this flow, we will add several species undergoing advection-diffusion, which is dominated by the downstream fluid flow in the channel.\n",
+    "The reaction will be included as a bulk-reaction, which means that the reaction can happen anywhere, the only requirement is that both species are present in the same lattice cell. When the reaction occurs, parts of the reactant species will turn into the product. How much of the species will transform within each timestep is determined by the respective reaction rate and the overall structure of the reaction.\n",
+    "\n",
+    "We start again by resetting the system and defining parameters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "system.ekcontainer = None\n",
+    "system.lb = None"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "BOX_L = [80, 32, 1] \n",
+    "AGRID = 1.0\n",
+    "DIFFUSION_COEFFICIENT = 0.01\n",
+    "TAU = 0.03\n",
+    "EXT_FORCE_DENSITY = [0.6, 0, 0]\n",
+    "OBSTACLE_RADIUS = 6\n",
+    "\n",
+    "DENSITY_FLUID = 0.5\n",
+    "VISCOSITY_KINEMATIC = 2.0\n",
+    "KT = 1.0\n",
+    "\n",
+    "TOTAL_FRAMES = 50\n",
+    "\n",
+    "system.time_step = TAU\n",
+    "system.cell_system.skin = 0.4\n",
+    "system.box_l = BOX_L"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lattice = espressomd.lb.LatticeWalberla(n_ghost_layers=1, agrid=AGRID)\n",
+    "lbf = espressomd.lb.LBFluidWalberla(\n",
+    "                lattice=lattice, density=DENSITY_FLUID, kinematic_viscosity=VISCOSITY_KINEMATIC,\n",
+    "                tau=TAU, ext_force_density=EXT_FORCE_DENSITY, kT=KT, seed=42)\n",
+    "system.lb = lbf\n",
+    "system.thermostat.set_lb(LB_fluid=lbf, seed=42)",
+    "\n",
+    "eksolver = espressomd.electrokinetics.EKNone(lattice=lattice)\n",
+    "system.ekcontainer = espressomd.electrokinetics.EKContainer(tau=TAU, solver=eksolver)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we can focus on the reactions. In this tutorial we choose the simple case of $A + B \\rightarrow C$, which means that equal parts of the educt species $A$ and $B$ can turn into the product species $C$.\n",
+    "ESPResSo distinguishes between educts and products by the sign of their respective stoechiometric coefficients, where educts have negative coefficients and products positive coefficients. Intuitively this can be understood that when a reaction happens, the density of the educts will decrease, hence the stoechiometric coefficient is negative. \n",
+    "\n",
+    "The reaction rate constant $r$ is the rate at which the reaction happens. The order $O_i$ for a species $i$ specifies to which order the reaction depends on the density of that species. Positive orders mean that the reaction is faster the more density of this species is present, for negative orders the reaction slows down with higher density.\n",
+    "In general, this process can be written as $\\Gamma = r \\left[ A \\right]^{O_A} \\left[ B \\right]^{O_B} \\left[ C \\right]^{O_C}$, where $\\Gamma$ is known as the reaction rate. This is sometimes also called the [rate equation](https://en.wikipedia.org/wiki/Rate_equation).\n",
+    "\n",
+    "For our specific simulation this means that all stoechiometric coefficients are $-1$ for the educts and $+1$ for the product. We choose the order of the educts as $1$ and the order of the product as $0$. This means that the more amount of both educts is present, the more will react and the amount of product present won't have an influence."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "REACTION_RATE_CONSTANT = 2.5\n",
+    "EDUCT_COEFFS = [-1, -1]\n",
+    "EDUCT_ORDERS = [1,1]\n",
+    "PRODUCT_COEFFS = [1]\n",
+    "PRODUCT_ORDERS = [0]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We create each involved species and directly specify their boundary-conditions for the domain-boundaries. We set the initial density of the species to 0 and also add Dirichlet boundary conditions of zero density at both the inlet and the outlet of the system."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "educt_species = []\n",
+    "product_species = []\n",
+    "reactants = []\n",
+    "for coeff, order in zip(EDUCT_COEFFS, EDUCT_ORDERS):\n",
+    "    species = espressomd.electrokinetics.EKSpecies(\n",
+    "        lattice=lattice, density=0.0, kT=KT,\n",
+    "        diffusion=DIFFUSION_COEFFICIENT, valency=0.0,\n",
+    "        advection=True, friction_coupling=True,\n",
+    "        ext_efield=[0., 0., 0.], tau=TAU)\n",
+    "    system.ekcontainer.add(species)\n",
+    "    reactants.append(\n",
+    "        espressomd.electrokinetics.EKReactant(\n",
+    "            ekspecies=species,\n",
+    "            stoech_coeff=coeff,\n",
+    "            order=order))\n",
+    "    educt_species.append(species)\n",
+    "    species[0,:,:].density_boundary = espressomd.electrokinetics.DensityBoundary(0.0)\n",
+    "    species[-1,:,:].density_boundary = espressomd.electrokinetics.DensityBoundary(0.0)\n",
+    "\n",
+    "for coeff, order in zip(PRODUCT_COEFFS, PRODUCT_ORDERS):\n",
+    "    species = espressomd.electrokinetics.EKSpecies(\n",
+    "        lattice=lattice, density=0.0, diffusion=DIFFUSION_COEFFICIENT,\n",
+    "        kT=KT, valency=0.0, advection=True, friction_coupling=True,\n",
+    "        ext_efield=[0., 0., 0.], tau=TAU)\n",
+    "    system.ekcontainer.add(species)\n",
+    "    reactants.append(\n",
+    "        espressomd.electrokinetics.EKReactant(\n",
+    "            ekspecies=species,\n",
+    "            stoech_coeff=coeff,\n",
+    "            order=order))\n",
+    "    product_species.append(species)\n",
+    "    species[0,:,:].density_boundary = espressomd.electrokinetics.DensityBoundary(0.0)\n",
+    "    species[-1,:,:].density_boundary = espressomd.electrokinetics.DensityBoundary(0.0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "solution2": "hidden",
+    "solution2_first": true
+   },
+   "source": [
+    "# Exercise:\n",
+    "- Create an instance of [`EKBulkReaction`](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKBulkReaction) using the previously created `reactants` and activate the reaction by adding it to [`system.ekcontainer.reactions`](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKContainer.reactions).\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "solution2": "hidden"
+   },
+   "source": [
+    "```python\n",
+    "reaction = espressomd.electrokinetics.EKBulkReaction(\n",
+    "    reactants=reactants, coefficient=REACTION_RATE_CONSTANT, lattice=lattice, tau=TAU)\n",
+    "\n",
+    "system.ekcontainer.reactions.add(reaction)\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The next thing to add to the system is the cylindrical obstacles, which act as the boundaries for the Kármán vortices to form. These are placed close to the inlet of the system and also act as impenetrable boundaries for the species.\n",
+    "Since ESPResSo uses periodic boundary conditions, we need to add a total of three cylinders to the system, which will form two complete cylinders in the periodic system."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cylinder_centers = [\n",
+    "    [BOX_L[0] // 10, 0, 1],\n",
+    "    [BOX_L[0] // 10, BOX_L[1] // 2, 1],\n",
+    "    [BOX_L[0] // 10, BOX_L[1], 1],    \n",
+    "]\n",
+    "\n",
+    "shape_cylinder = []\n",
+    "for cylinder_center in cylinder_centers:\n",
+    "    shape_cylinder.append(espressomd.shapes.Cylinder(\n",
+    "        center=cylinder_center,\n",
+    "        axis=[0, 0, 1],\n",
+    "        length=BOX_L[2],\n",
+    "        radius=OBSTACLE_RADIUS,\n",
+    "    ))\n",
+    "\n",
+    "for shape in shape_cylinder:\n",
+    "    lbf.add_boundary_from_shape(shape)\n",
+    "    for spec in (*educt_species, *product_species):\n",
+    "        spec.add_boundary_from_shape(shape, value=[0,0,0], boundary_type=espressomd.electrokinetics.FluxBoundary)\n",
+    "        spec.add_boundary_from_shape(shape, value=0., boundary_type=espressomd.electrokinetics.DensityBoundary)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Up to this point there is no species present anywhere in the system and also no way for it to enter the system. Since the reaction is irreversible in our setup, we need to introduce some density of both the educt species to the system.\n",
+    "For that we set two additional Dirichlet boundary conditions (sources) in the domain, where we fix the species' density to a constant, non-zero value. The sources are placed some distance apart such that the reaction happens further downstream when diffusion mixes the two species."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "source_boundary = espressomd.electrokinetics.DensityBoundary(10.0)\n",
+    "source_x_pos = BOX_L[1] // 20\n",
+    "\n",
+    "educt_species[0][source_x_pos,BOX_L[1]//4-1:BOX_L[1]//4+1,:].density_boundary = source_boundary\n",
+    "educt_species[1][source_x_pos,3*(BOX_L[1]//4)-1:3*(BOX_L[1]//4)+1,:].density_boundary = source_boundary"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "With this, the system is now finally complete and we can start the integration. To see the system evolve, we will render a movie from the timeseries of the system. For that we have to setup some helper functions for the plotting, which are beyond the scope of this tutorial."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "VIDEO_TAG = \"\"\"<video controls>\n",
+    " <source src=\"data:video/x-m4v;base64,{0}\" type=\"video/mp4\">\n",
+    " Your browser does not support the video tag.\n",
+    "</video>\"\"\"\n",
+    "\n",
+    "# set ignore 'divide' and 'invalid' errors\n",
+    "# these occur when plotting the flowfield containing a zero velocity\n",
+    "np.seterr(divide='ignore', invalid='ignore')\n",
+    "\n",
+    "def anim_to_html(anim):\n",
+    "    if not hasattr(anim, '_encoded_video'):\n",
+    "        with tempfile.NamedTemporaryFile(suffix='.mp4') as f:\n",
+    "            anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264'])\n",
+    "            with open(f.name, \"rb\") as g:\n",
+    "                video = g.read()\n",
+    "        anim._encoded_video = base64.b64encode(video).decode('ascii')\n",
+    "        plt.close(anim._fig)\n",
+    "    return VIDEO_TAG.format(anim._encoded_video)\n",
+    "\n",
+    "animation.Animation._repr_html_ = anim_to_html"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "get_colormap = mpl.colormaps.get_cmap if hasattr(mpl.colormaps, \"get_cmap\") else mpl.cm.get_cmap\n",
+    "box_width = lattice.shape[1]\n",
+    "box_height = lattice.shape[0]\n",
+    "\n",
+    "boundary_mask = lbf[:, :, 0].boundary != None\n",
+    "\n",
+    "cmap = get_colormap(\"viridis\").copy()\n",
+    "cmap.set_bad(color=\"gray\")\n",
+    "cmap_quiver = get_colormap(\"binary\").copy()\n",
+    "cmap_quiver.set_bad(color=\"gray\")\n",
+    "\n",
+    "# setup figure and prepare axes\n",
+    "fig = plt.figure(figsize=(9.8, 5.5))\n",
+    "imshow_kwargs = {\"origin\": \"upper\", \"extent\": (0, BOX_L[1], BOX_L[0], 0)}\n",
+    "gs = fig.add_gridspec(1, 4, wspace=0.1)\n",
+    "ax1 = plt.subplot(gs[0])\n",
+    "ax2 = plt.subplot(gs[1], sharey=ax1)\n",
+    "ax3 = plt.subplot(gs[2], sharey=ax1)\n",
+    "ax4 = plt.subplot(gs[3], sharey=ax1)\n",
+    "ax1.set_yticks(np.arange(0, 80 + 1, 16))\n",
+    "ax1.set_xticks(np.arange(0, 32 + 1, 16))\n",
+    "ax2.set_xticks(np.arange(0, 32 + 1, 16))\n",
+    "ax3.set_xticks(np.arange(0, 32 + 1, 16))\n",
+    "ax4.set_xticks(np.arange(0, 32 + 1, 16))\n",
+    "\n",
+    "# set the background color for the quiver plot\n",
+    "bg_colors = np.copy(boundary_mask).astype(float)\n",
+    "bg_colors[boundary_mask] = np.NaN\n",
+    "ax4.imshow(bg_colors, cmap=cmap_quiver, **imshow_kwargs)\n",
+    "\n",
+    "for ax, title in zip(\n",
+    "    [ax1, ax2, ax3, ax4],\n",
+    "    [\"educt 1\", \"educt 2\", \"product\", \"fluid velocity\"]\n",
+    "):\n",
+    "    ax.set_title(title)\n",
+    "    ax.set_xlim((0, box_width))\n",
+    "    ax.set_ylim((0, box_height))\n",
+    "\n",
+    "# create meshgrid for quiver plot subsampled by a factor 2\n",
+    "xs = np.arange(0, box_width, 2)\n",
+    "ys = np.arange(0, box_height, 2)\n",
+    "X, Y = np.meshgrid(xs, ys)\n",
+    "\n",
+    "flowfield = lbf[:, :, 0].velocity[::2, ::2, :]\n",
+    "quiver = ax4.quiver(X + 1., Y + 1., flowfield[..., 1], flowfield[..., 0], scale=100)\n",
+    "fig.subplots_adjust(left=0.025, bottom=0.075, right=0.975, top=0.925, wspace=0.0125)\n",
+    "\n",
+    "progress_bar = tqdm.tqdm(total=TOTAL_FRAMES)\n",
+    "\n",
+    "def draw_frame(frame):\n",
+    "    system.integrator.run(50)\n",
+    "    \n",
+    "    flowfield = np.copy(lbf[:, :, 0].velocity)\n",
+    "    e1 = np.copy(educt_species[0][:, :, 0].density)\n",
+    "    e2 = np.copy(educt_species[1][:, :, 0].density)\n",
+    "    p = np.copy(product_species[0][:, :, 0].density)\n",
+    "    \n",
+    "    # apply the mask for the boundary\n",
+    "    e1[boundary_mask] = np.NaN\n",
+    "    e2[boundary_mask] = np.NaN\n",
+    "    p[boundary_mask] = np.NaN\n",
+    "    flowfield[boundary_mask] = np.NaN\n",
+    "\n",
+    "    ax1.imshow(e1, cmap=cmap, vmin=0., vmax=source_boundary.density, **imshow_kwargs)\n",
+    "    ax2.imshow(e2, cmap=cmap, vmin=0., vmax=source_boundary.density, **imshow_kwargs)\n",
+    "    ax3.imshow(p, cmap=cmap, vmin=0., vmax=source_boundary.density, **imshow_kwargs)\n",
+    "    quiver.set_UVC((flowfield[::2, ::2, 1] + flowfield[1::2, 1::2, 1]) / 2.,\n",
+    "                   (flowfield[::2, ::2, 0] + flowfield[1::2, 1::2, 0]) / 2.)\n",
+    "    \n",
+    "    progress_bar.update()\n",
+    "\n",
+    "\n",
+    "animation.FuncAnimation(fig, draw_frame, frames=range(TOTAL_FRAMES), interval=300)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Looking at the movie of the species densities one can see that the fluid flow advects the educt species from their source locations past the cylinders into the system. Here, they start to mix and react, such that the product forms.\n",
+    "The vortex street behind the obstacles enhances mixing through fluid turbulence. The density of the product then increases towards the outflow-location of the channel, where it is deleted because of our zero-density boundary condition."
+   ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -671,7 +1291,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.6"
+   "version": "3.10.12"
   }
  },
  "nbformat": 4,
diff --git a/doc/tutorials/electrokinetics/scripts/eof_analytical.py b/doc/tutorials/electrokinetics/scripts/eof_analytical.py
deleted file mode 100644
index 83122d225a8..00000000000
--- a/doc/tutorials/electrokinetics/scripts/eof_analytical.py
+++ /dev/null
@@ -1,146 +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 <http://www.gnu.org/licenses/>.
-# Set the slit pore geometry the width is the non-periodic part of the geometry
-# the padding is used to ensure that there is no field outside the slit
-import numpy as np
-
-box_x = 6
-box_y = 6
-width = 50
-
-padding = 6
-box_z = width + 2 * padding
-
-# Set the electrokinetic parameters
-
-agrid = 1.0
-temperature = 1.0
-bjerrum_length = 0.7095
-valency = 1.0
-viscosity_dynamic = 79.53
-density_water = 26.15
-sigma = -0.05
-force = 0.1
-
-viscosity_kinematic = viscosity_dynamic / density_water
-density_counterions = -2.0 * sigma / width
-
-# Calculate the inverse length xi, which is a combination of various
-# constants (xi = zeC/2kBT), with C a constant that needs to be
-# solved for, or equivalently, xi needs to be solved for
-
-# root finding function
-
-
-def solve(xi=None, d=None, bjerrum_length=None, sigma=None, valency=None):
-    el_char = 1.0
-    return xi * np.tan(xi * d / 2.0) + 2.0 * np.pi * \
-        bjerrum_length * sigma / (valency * el_char)
-
-
-size = np.pi / (2.0 * width)
-
-pnt0 = 0.0
-pntm = pnt0 + size
-pnt1 = pnt0 + 1.9 * size
-
-# the bisection scheme
-tol = 1e-8
-
-while size > tol:
-    val0 = solve(xi=pnt0, d=width,
-                 bjerrum_length=bjerrum_length, sigma=sigma, valency=valency)
-    val1 = solve(xi=pnt1, d=width,
-                 bjerrum_length=bjerrum_length, sigma=sigma, valency=valency)
-    valm = solve(xi=pntm, d=width,
-                 bjerrum_length=bjerrum_length, sigma=sigma, valency=valency)
-
-    if (val0 < 0.0) and (val1 > 0.0):
-        if valm < 0.0:
-            pnt0 = pntm
-            size = size / 2.0
-            pntm = pnt0 + size
-        else:
-            pnt1 = pntm
-            size = size / 2.0
-            pntm = pnt1 - size
-    elif (val0 > 0.0) and (val1 < 0.0):
-        if valm < 0.0:
-            pnt1 = pntm
-            size = size / 2.0
-            pntm = pnt0 - size
-        else:
-            pnt0 = pntm
-            size = size / 2.0
-            pntm = pnt0 - size
-    else:
-        raise Exception(
-            "Bisection method fails:\nTuning of domain boundaries may be required.")
-
-
-# obtain the desired xi value
-xi = pntm
-
-# function to calculate the density
-
-
-def density(x=None, xi=None, bjerrum_length=None):
-    return (xi**2) / (2.0 * np.pi * bjerrum_length * np.cos(xi * x)**2)
-
-# function to calculate the velocity
-
-
-def velocity(x=None, xi=None, d=None, bjerrum_length=None, force=None,
-             viscosity_kinematic=None, density_water=None):
-    return force * np.log(np.cos(xi * x) / np.cos(xi * d / 2.0)) / \
-        (2.0 * np.pi * bjerrum_length * viscosity_kinematic * density_water)
-
-# function to calculate the nonzero component of the pressure tensor
-
-
-def pressure_tensor_offdiagonal(x=None, xi=None, bjerrum_length=None,
-                                force=None):
-    return force * xi * np.tan(xi * x) / (2.0 * np.pi * bjerrum_length)
-
-
-position_list = []
-density_list = []
-velocity_list = []
-pressure_xy_list = []
-
-for i in range(int(box_z / agrid)):
-    if (i * agrid >= padding) and (i * agrid < box_z - padding):
-        position = i * agrid - padding - width / 2.0 + agrid / 2.0
-        position_list.append(position)
-
-        # density
-        density_list.append(density(x=position, xi=xi,
-                                    bjerrum_length=bjerrum_length))
-
-        # velocity
-        velocity_list.append(velocity(x=position, xi=xi,
-                                      d=width, bjerrum_length=bjerrum_length,
-                                      force=force, viscosity_kinematic=viscosity_kinematic,
-                                      density_water=density_water))
-        # xz component pressure tensor
-        pressure_xy_list.append(pressure_tensor_offdiagonal(x=position, xi=xi,
-                                                            bjerrum_length=bjerrum_length, force=force))
-
-np.savetxt(
-    "eof_analytical.dat", np.column_stack(
-        (position_list, density_list, velocity_list, pressure_xy_list)),
-    header="#position calculated_density calculated_velocity calculated_pressure_xy")
diff --git a/src/python/espressomd/electrokinetics.py b/src/python/espressomd/electrokinetics.py
index 1fb2a59ae9a..336cea4780d 100644
--- a/src/python/espressomd/electrokinetics.py
+++ b/src/python/espressomd/electrokinetics.py
@@ -31,7 +31,16 @@
 class EKFFT(ScriptInterfaceHelper):
     """
     A FFT-based Poisson solver.
+    Intrinsically assumes periodic boundary conditions.
 
+    Parameters
+    ----------
+    lattice : :obj:`espressomd.lb.LatticeWalberla <espressomd.detail.walberla.LatticeWalberla>`
+        Lattice object.
+    permittivity : :obj:`float`
+        permittivity of the fluid :math:`\\epsilon_0 \\epsilon_{\\mathrm{r}}`.
+    single_precision : :obj:`bool`, optional
+        Use single-precision floating-point arithmetic.
     """
     _so_name = "walberla::EKFFT"
     _so_features = ("WALBERLA_FFT",)
@@ -44,6 +53,12 @@ class EKNone(ScriptInterfaceHelper):
     The default Poisson solver.
     Imposes a null electrostatic potential everywhere.
 
+    Parameters
+    ----------
+    lattice : :obj:`espressomd.lb.LatticeWalberla <espressomd.detail.walberla.LatticeWalberla>`
+        Lattice object.
+    single_precision : :obj:`bool`, optional
+        Use single-precision floating-point arithmetic.
     """
     _so_name = "walberla::EKNone"
     _so_features = ("WALBERLA",)
@@ -576,16 +591,61 @@ def __repr__(self):
 
 @script_interface_register
 class EKReactant(ScriptInterfaceHelper):
+    """
+    Reactant-object which specifies the contribution of a species to a reaction.
+
+    Parameters
+    ----------
+    ekspecies : :obj:`~espressomd.electrokinetics.EKSpecies`
+        EK species to react
+    stoech_coeff: :obj:`float`
+        Stoechiometric coefficient of this species in the reaction.
+        Products have positive coefficients whereas educts have negative ones.
+    order: :obj:`float`
+        Partial-order of this species in the reaction.
+
+    """
     _so_name = "walberla::EKReactant"
     _so_creation_policy = "GLOBAL"
 
 
 class EKBulkReaction(ScriptInterfaceHelper):
+    """
+    Reaction type that is applied everywhere in the domain.
+
+    Parameters
+    ----------
+    lattice : :obj:`espressomd.electrokinetics.LatticeWalberla <espressomd.detail.walberla.LatticeWalberla>`
+        Lattice object.
+    tau : :obj:`float`
+        EK time step, must be an integer multiple of the MD time step.
+    coefficient : :obj:`float`
+        Reaction rate constant of the reaction.
+    reactants: array_like of :obj:`~espressomd.electrokinetics.EKReactant`
+        Reactants that participate this reaction.
+
+    """
     _so_name = "walberla::EKBulkReaction"
     _so_creation_policy = "GLOBAL"
 
 
 class EKIndexedReaction(ScriptInterfaceHelper):
+    """
+    Reaction type that is applied only on specific cells in the domain.
+    Can be used to model surface-reactions.
+
+    Parameters
+    ----------
+    lattice : :obj:`espressomd.electrokinetics.LatticeWalberla <espressomd.detail.walberla.LatticeWalberla>`
+        Lattice object.
+    tau : :obj:`float`
+        EK time step, must be an integer multiple of the MD time step.
+    coefficient : :obj:`float`
+        Reaction rate constant of the reaction.
+    reactants: array_like of :obj:`~espressomd.electrokinetics.EKReactant`
+        Reactants that participate this reaction.
+
+    """
     _so_name = "walberla::EKIndexedReaction"
     _so_creation_policy = "GLOBAL"
 
@@ -648,21 +708,64 @@ def __setitem__(self, key, values):
 
 @script_interface_register
 class EKReactions(ScriptObjectList):
+    """
+    Container object holding all EK-reactions that are considered.
+
+    Methods
+    -------
+    clear()
+        Remove all reactions.
+
+    """
     _so_name = "walberla::EKReactions"
     _so_creation_policy = "GLOBAL"
     _so_bind_methods = ("clear",)
 
     def add(self, reaction):
+        """
+        Add a reaction to the container.
+
+        Parameters
+        ----------
+        reaction : :obj:`~espressomd.electrokinetics.EKBulkReaction` or :obj:`~espressomd.electrokinetics.EKIndexedReaction`
+            Reaction to be added.
+
+        """
         if not isinstance(reaction, (EKBulkReaction, EKIndexedReaction)):
             raise TypeError("reaction object is not of correct type.")
         self.call_method("add", object=reaction)
 
     def remove(self, reaction):
+        """
+        Remove a reaction from the container.
+
+        Parameters
+        ----------
+        reaction : :obj:`~espressomd.electrokinetics.EKBulkReaction` or :obj:`~espressomd.electrokinetics.EKIndexedReaction`
+            Reaction to be removed.
+
+        """
         self.call_method("remove", object=reaction)
 
 
 @script_interface_register
 class EKContainer(ScriptObjectList):
+    """
+    Container object holding the :obj:`~espressomd.electrokinetics.EKSpecies`.
+
+    Parameters
+    ----------
+    tau : :obj:`float`
+        EK time step, must be an integer multiple of the MD time step.
+    solver : :obj:`~espressomd.electrokinetics.EKNone` or :obj:`~espressomd.electrokinetics.EKFFT`
+        Solver defining the treatment of the electrostatic Poisson-equation.
+
+    Methods
+    -------
+    clear()
+        Removes all species.
+
+    """
     _so_name = "walberla::EKContainer"
     _so_creation_policy = "GLOBAL"
     _so_features = ("WALBERLA",)
@@ -675,8 +778,36 @@ def __init__(self, *args, **kwargs):
         else:
             super().__init__(**kwargs)
 
+    @property
+    def reactions(self):
+        """
+        Returns
+        -------
+        Reactions-container of the current reactions (:obj:`~espressomd.electrokinetics.EKReactions`).
+
+        """
+        return self._getter("reactions")
+
     def add(self, ekspecies):
+        """
+        Add an :obj:`~espressomd.electrokinetics.EKSpecies` to the container.
+
+        Parameters
+        ----------
+        ekspecies : :obj:`~espressomd.electrokinetics.EKSpecies`
+            Species to be added.
+
+        """
         self.call_method("add", object=ekspecies)
 
     def remove(self, ekspecies):
+        """
+        Remove an :obj:`~espressomd.electrokinetics.EKSpecies` from the container.
+
+        Parameters
+        ----------
+        ekspecies : :obj:`~espressomd.electrokinetics.EKSpecies`
+            Species to be removed.
+
+        """
         self.call_method("remove", object=ekspecies)
diff --git a/testsuite/scripts/tutorials/test_electrokinetics.py b/testsuite/scripts/tutorials/test_electrokinetics.py
index 611f8dafceb..e138fe209ad 100644
--- a/testsuite/scripts/tutorials/test_electrokinetics.py
+++ b/testsuite/scripts/tutorials/test_electrokinetics.py
@@ -1,4 +1,5 @@
-# Copyright (C) 2019-2022 The ESPResSo project
+#
+# Copyright (C) 2019-2023 The ESPResSo project
 #
 # This file is part of ESPResSo.
 #
@@ -14,35 +15,109 @@
 #
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
 
 import unittest as ut
 import importlib_wrapper as iw
 import numpy as np
+import scipy.optimize
 
 tutorial, skipIfMissingFeatures = iw.configure_and_import(
-    "@TUTORIALS_DIR@/electrokinetics/electrokinetics.py", integration_length=400)
+    "@TUTORIALS_DIR@/electrokinetics/electrokinetics.py", TOTAL_FRAMES=0)
 
 
 @skipIfMissingFeatures
 class Tutorial(ut.TestCase):
     system = tutorial.system
 
-    def normalize_two_datasets(self, a, b):
-        offset = min(np.min(a), np.min(b))
-        a -= offset
-        b -= offset
-        scale = max(np.max(a), np.max(b))
-        a /= scale
-        b /= scale
-
-    def test_simulation(self):
-        for varname, tol in zip(["density", "velocity"], [2, 5]):
-            sim = np.array(tutorial.__dict__[varname + "_list"])
-            ana = np.array(tutorial.eof_analytical.__dict__[varname + "_list"])
-            self.normalize_two_datasets(sim, ana)
-            accuracy = np.max(np.abs(sim - ana))
-            # expecting at most a few percents deviation
-            self.assertLess(accuracy, tol / 100.)
+    def test_analytic_solutions(self):
+        # check advection-diffusion
+        mu_simulated = tutorial.positions_diagonal[np.argmax(
+            tutorial.values_diagonal)]
+        mu_analytic = tutorial.positions_analytic[np.argmax(
+            tutorial.values_analytic)]
+        tol = 2. * tutorial.AGRID
+        self.assertAlmostEqual(mu_simulated, mu_analytic, delta=tol)
+        # check electroosmotic flow
+        np.testing.assert_allclose(
+            np.copy(tutorial.density_eof),
+            tutorial.analytic_density_eof,
+            rtol=0.005)
+        np.testing.assert_allclose(
+            np.copy(tutorial.velocity_eof),
+            tutorial.analytic_velocity_eof,
+            rtol=0.05)
+        np.testing.assert_allclose(
+            np.copy(tutorial.pressure_tensor_eof),
+            tutorial.analytic_pressure_tensor_eof,
+            rtol=0.05)
+        # check pressure-driven flow
+        np.testing.assert_allclose(
+            np.copy(tutorial.density_pressure),
+            tutorial.analytic_density_eof,
+            rtol=0.005)
+        np.testing.assert_allclose(
+            np.copy(tutorial.velocity_pressure),
+            tutorial.analytic_velocity_pressure,
+            rtol=0.05)
+        np.testing.assert_allclose(
+            np.copy(tutorial.pressure_tensor_pressure),
+            tutorial.analytic_pressure_tensor_pressure,
+            rtol=0.05)
+
+    def test_karman_vortex_street(self):
+        """
+        Check for the formation of a Kármán vortex street. Due to turbulence,
+        a wavy pattern emerges.
+        """
+        def get_frequency_species(species):
+            """
+            Compute the principal wavevector of a chemical species.
+            """
+            fdata = np.zeros(16)
+            for i in range(32, 80):
+                rdata = species[i, :, 0].density
+                fdata += np.abs(np.fft.fft(rdata - np.mean(rdata)))[:16]
+            return np.argmax(fdata)
+
+        def get_phase_karman(species):
+            """
+            Compute the time-dependent phase of a turbulent flow profile.
+            """
+            phase = []
+            k = 2  # wavevector for product species
+            for i in range(36, 68):
+                rdata = species[i, :, 0].density
+                phase.append(np.angle(np.fft.fft(rdata - np.mean(rdata))[k]))
+            return np.array(phase)
+
+        def cosine_kernel(x, magnitude, freq, phase):
+            return magnitude * np.cos(x * freq + phase)
+
+        tutorial.system.integrator.run(2000)
+        # check for finite values
+        for species in (*tutorial.educt_species, *tutorial.product_species):
+            assert np.all(np.isfinite(species[:, :, :].density))
+            assert np.all(species[:, :, :].density >= 0)
+        assert np.all(np.isfinite(tutorial.lbf[:, :, :].velocity))
+        # there is only one inlet per educt, thus wavelength == box width
+        self.assertEqual(get_frequency_species(tutorial.educt_species[0]), 1)
+        self.assertEqual(get_frequency_species(tutorial.educt_species[1]), 1)
+        # reaction happens in the mixing region, thus the frequency doubles
+        self.assertEqual(get_frequency_species(tutorial.product_species[0]), 2)
+        # check for turbulence onset
+        ref_params = np.array([2., 0.12, 1. / 2. * np.pi])
+        sim_phase = get_phase_karman(tutorial.product_species[0])
+        xdata = np.arange(sim_phase.shape[0])
+        popt, _ = scipy.optimize.curve_fit(
+            cosine_kernel, xdata, sim_phase, p0=ref_params,
+            bounds=([-4., 0.08, 0.], [4., 0.24, 2. * np.pi]))
+        fit_phase = cosine_kernel(xdata, *popt)
+        rmsd = np.sqrt(np.mean(np.square(sim_phase - fit_phase)))
+        self.assertAlmostEqual(popt[2], ref_params[2], delta=0.20)
+        self.assertAlmostEqual(popt[1], ref_params[1], delta=0.02)
+        self.assertAlmostEqual(popt[0], ref_params[0], delta=0.80)
+        self.assertLess(rmsd / abs(popt[0]), 0.2)
 
 
 if __name__ == "__main__":