Tutorial 4 - Lattice-Boltzmann (#3893)

Description of changes:
- Update part 1 (fixes #3828)
- Write solutions for parts 2 and 3 and check them in CI
- Calculate the diffusion coefficient of a polymer with the MSD and Green-Kubo methods
- Replace LB hydrodynamics by LD hydrodynamics for polymer diffusion

doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb

@@ -9,11 +9,12 @@
     "#### Before you start:\n",
     "\n",
     "With this tutorial you can get started using the lattice-Boltzmann method\n",
-    "for scientific applications. We give a brief introduction about the theory\n",
+    "for hydrodynamics. We give a brief introduction about the theory\n",
     "and how to use the method in **ESPResSo**. We have selected two interesting problems for\n",
     "which LB can be applied and which are well understood. You can start with any of them.\n",
     "\n",
-    "The tutorial is relatively long and working through it carefully is work for at least a full day. You can however get a glimpse of different aspects by starting to work on the tasks.\n",
+    "The tutorial is relatively long and working through it carefully will take at least a full day.\n",
+    "You can however get a glimpse of different aspects by starting to work on the tasks.\n",
     "\n",
     "Note: LB cannot be used as a black box. It is unavoidable to spend time\n",
     "learning the theory and gaining practical experience."
@@ -28,14 +29,9 @@
     "In this tutorial, you will learn the basics of the lattice-Boltzmann method (LBM) with\n",
     "special focus on the application on soft matter simulations or more precisely on how to\n",
     "apply it in combination with molecular dynamics to take into account hydrodynamic\n",
-    "solvent effects without the need to introduce thousands of solvent particles.\n",
-    "\n",
-    "The LBM – its theory as well as its applications – is still a very active field of research.\n",
-    "After almost 20 years of development there are many cases in which the LBM has proven\n",
-    "to be fruitful, in other cases the LBM is considered promising, and in some cases it has\n",
-    "not been of any help. We encourage you to contribute to the scientific discussion of\n",
-    "the LBM because there is still a lot that is unknown or only vaguely known about this\n",
-    "fascinating method.\n",
+    "solvent effects without the need to introduce thousands of solvent particles. The LBM\n",
+    "– its theory as well as its applications – is a very active field of research with\n",
+    "more than 30 years of development.\n",
     "\n",
     "### Tutorial Outline\n",
     "\n",
@@ -46,7 +42,7 @@
     "\n",
     "* **Polymer diffusion.** We analyze the length dependence of the diffusion of polymers.\n",
     "\n",
-    "* **Poiseuille flow.** We reproduce the flow profile between two walls.\n",
+    "* **Poiseuille flow.** We reproduce the flow profile of a fluid moving through a pipe under a homogeneous force density.\n",
     "\n",
     "### Notes on the **ESPResSo** version you will need\n",
     "\n",
@@ -55,6 +51,7 @@
     "\n",
     "For the tutorial you will have to compile in the following features:\n",
     "```c++\n",
+    "#define LB_BOUNDARIES\n",
     "#define LB_BOUNDARIES_GPU\n",
     "#define LENNARD_JONES\n",
     "```\n",
@@ -96,7 +93,7 @@
     "\n",
     "The LBM discretizes the Boltzmann equation not only in real space (the lattice!) and\n",
     "time, but also the velocity space is discretized. A surprisingly small number of velocities,\n",
-    "in 3D usually 19, is sufficient to describe incompressible, viscous flow correctly. Mostly\n",
+    "usually 19 in 3D, is sufficient to describe incompressible, viscous flow correctly. Mostly\n",
     "we will refer to the three-dimensional model with a discrete set of 19 velocities, which is\n",
     "conventionally called D3Q19. These velocities, $\\vec{c_i}$ , are chosen such that they correspond to\n",
     "the movement from one lattice node to another in one time step. A two step scheme is\n",
@@ -149,7 +146,7 @@
     "why the mode space is helpful. A 19 dimensional matrix that\n",
     "conserves the first 4 modes (with the eigenvalue 1) is diagonal in the\n",
     "first four rows and columns.\n",
-    "Some struggling with lattice symmetries shows that four independent\n",
+    "By symmetry consideration one finds that only four independent\n",
     "variables are enough to characterize the linear relaxation\n",
     "process so that all symmetries of the lattice are obeyed. \n",
     "Two of them are closely related to \n",
@@ -171,9 +168,7 @@
     "In other words: each mode is relaxed towards\n",
     "its equilibrium value with a relaxation rate $\\gamma_i$.\n",
     "The conserved modes are not relaxed, i.e. their relaxation rate is 1.\n",
-    "\n",
-    "By symmetry consideration one finds that only four independent\n",
-    "relaxation rates are allowed. We summarize them here.\n",
+    "We summarize them here:\n",
     "\n",
     "\\begin{align*}\n",
     "  m^\\star_i &= \\gamma_i m_i  \\\\\n",
@@ -188,7 +183,8 @@
     "random fluctuations are added to the nonconserved modes $5\\dots 19$ on every LB node so that\n",
     "the LB fluid temperature is well defined and the corresponding\n",
     "fluctuation formula holds, according to the fluctuation dissipation theorem.\n",
-    "An extensive discussion of this topic is found in [3]."
+    "An extensive discussion of this topic is found in [3].\n",
+    "Thermalization of the fluid is optional in **ESPResSo**."
    ]
   },
   {
@@ -197,7 +193,7 @@
    "source": [
     "### Particle coupling\n",
     "\n",
-    "In **ESPResSo** particles are coupled to the LB fluid via the so called the force coupling method:\n",
+    "In **ESPResSo** particles are coupled to the LB fluid via the so-called force coupling method:\n",
     "the fluid velocity $\\vec{u}$ at the position of a particle is calculated \n",
     "via trilinear interpolation and a force is applied on the particle\n",
     "that is proportional to the velocity difference between particle and fluid:\n",
@@ -224,8 +220,12 @@
    "source": [
     "## 3 The LB interface in ESPResSo\n",
     "\n",
-    "**ESPResSo** features two virtually independent implementations of LB. One implementation uses CPUs and one uses a GPU to perform the computational work. For this, we provide two actor classes <tt>LBFluid</tt> and <tt>LBFluidGPU</tt> in the module <tt>espressomd.lb</tt>,\n",
-    "as well as the optional <tt>LBBoundary</tt> class found in <tt>espressomd.lbboundaries</tt>.\n",
+    "**ESPResSo** features two virtually independent implementations of LB. One implementation uses CPUs and one uses a GPU to perform the computational work. For this, we provide two actor classes\n",
+    "[<tt>LBFluid</tt>](http://espressomd.org/html/doc/espressomd.html#espressomd.lb.LBFluid) and\n",
+    "[<tt>LBFluidGPU</tt>](http://espressomd.org/html/doc/espressomd.html#espressomd.lb.LBFluidGPU) in the module\n",
+    "[<tt>espressomd.lb</tt>](http://espressomd.org/html/doc/espressomd.html#module-espressomd.lb), as well as the optional\n",
+    "[<tt>LBBoundary</tt>](http://espressomd.org/html/doc/espressomd.html#espressomd.lbboundaries.LBBoundary) class found in\n",
+    "[<tt>espressomd.lbboundaries</tt>](http://espressomd.org/html/doc/espressomd.html#module-espressomd.lbboundaries).\n",
     "\n",
     "The LB lattice is a cubic lattice, with a lattice constant <tt>agrid</tt> that\n",
     "is the same in all spatial directions. The chosen box length must be an integer multiple\n",
@@ -268,9 +268,11 @@
     "+ <tt>agrid</tt>: The lattice constant of the fluid. It is used to determine the number of LB nodes per direction from <tt>box_l</tt>. *They have to be compatible.*\n",
     "+ <tt>visc</tt>: The kinematic viscosity\n",
     "+ <tt>tau</tt>: The time step of LB. It has to be an integer multiple of <tt>System.time_step</tt>.\n",
+    "+ <tt>kT</tt>: Thermal energy of the simulated heat bath for thermalized fluids, use 0 to deactivate thermalization.\n",
+    "+ <tt>seed</tt>: The random number generator seed, only relevant for thermalized fluids (i.e. <tt>kT</tt> \\> 0).\n",
     "+ <tt>ext_force_density</tt>: An external force density applied to every node. This is given as a list, tuple or array with three components.\n",
     "\n",
-    "Using these arguments, one can initialize an <tt>LBFluid</tt> object. This object then needs to be added to the system's actor list. The code below provides a minimum example.\n",
+    "Using these arguments, one can initialize an <tt>LBFluid</tt> object. This object then needs to be added to the system's actor list. The code below provides a minimal example.\n",
     "\n",
     "```python\n",
     "import espressomd\n",
@@ -281,7 +283,7 @@
     "system.time_step = 0.01\n",
     "system.cell_system.skin = 0.4\n",
     "\n",
-    "# Initialize and LBFluid with the minimum set of valid parameters.\n",
+    "# Initialize an LBFluid with the minimum set of valid parameters.\n",
     "lbf = lb.LBFluidGPU(agrid=1, dens=10, visc=.1, tau=0.01)\n",
     "# Activate the LB by adding it to the System's actor list.\n",
     "system.actors.add(lbf)\n",
@@ -313,15 +315,13 @@
    "source": [
     "### The <tt>LBBoundary</tt> class\n",
     "\n",
-    "The <tt>LBBoundary</tt> class represents a boundary on the <tt>LBFluid</tt>\n",
-    "lattice. It depends on the classes of the module <tt>espressomd.shapes</tt> as it derives its geometry from them. For the initialization, the arguments <tt>shape</tt> and <tt>velocity</tt> are supported. The <tt>shape</tt> argument takes an object from the <tt>shapes</tt> module and the <tt>velocity</tt> argument expects a list, tuple or array containing 3 floats. Setting the <tt>velocity</tt> will results in a slip boundary condition.\n",
+    "The [<tt>LBBoundary</tt>](http://espressomd.org/html/doc/espressomd.html#espressomd.lbboundaries.LBBoundary) class represents a boundary on the\n",
+    "[<tt>LBFluid</tt>](http://espressomd.org/html/doc/espressomd.html#espressomd.lb.LBFluid) lattice.\n",
+    "It depends on the classes of the module <tt>espressomd.shapes</tt> as it derives its geometry from them. For the initialization, the arguments <tt>shape</tt> and <tt>velocity</tt> are supported. The <tt>shape</tt> argument takes an object from the <tt>shapes</tt> module and the <tt>velocity</tt> argument expects a list, tuple or array containing 3 floats. Setting the <tt>velocity</tt> will result in a slip boundary condition.\n",
     "\n",
     "Note that the boundaries are not constructed through the periodic boundary. If, for example, one would set a sphere with its center in one of the corner of the boxes, a sphere fragment will be generated. To avoid this, make sure the sphere, or any other boundary, fits inside the central box.\n",
     "\n",
-    "This part of the LB implementation is still experimental, so please tell us about your experience with it. In general even the simple case of no-slip\n",
-    "boundary is still an important research topic in the LB community and in combination with point particle coupling not much experience exists. This means: do research on that topic, play around with parameters and figure out what happens.\n",
-    "\n",
-    "Boundaries are initialized by passing a shape object to the <tt>LBBoundary</tt> class. One way to initialize a wall is:\n",
+    "Boundaries are instantiated by passing a shape object to the <tt>LBBoundary</tt> class. Here is one way to construct a wall and add it to an existing `system` instance:\n",
     "\n",
     "```python\n",
     "import espressomd.lbboundaries\n",
@@ -332,7 +332,7 @@
     "system.lbboundaries.add(wall)\n",
     "```\n",
     "\n",
-    "Note that all used variables are inherited from previous examples. This will create a wall with a surface normal of $(1, 0, 0)$ at a distance of 1 from the origin of the coordinate system in direction of the normal vector. The wall exhibits a slip boundary condition with a velocity of $(0, 0, 0.01)$. For a no-slip boundary condition, leave out the velocity argument or set it to zero. Please refer to the user guide for a complete list of constraints.\n",
+    "This will create a wall with a surface normal of $(1, 0, 0)$ at a distance of 1 from the origin of the coordinate system in direction of the normal vector. The wall exhibits a slip boundary condition with a velocity of $(0, 0, 0.01)$. For a no-slip boundary condition, leave out the velocity argument or set it to zero. Please refer to the user guide for a complete list of constraints.\n",
     "\n",
     "In **ESPResSo** the so-called *link bounce back* method is implemented, where the effective hydrodynamic boundary is located midway between boundary and fluid node."
    ]