Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Tutorial 4 - Lattice-Boltzmann #3893

Merged
merged 19 commits into from
Oct 8, 2020
Merged

Tutorial 4 - Lattice-Boltzmann #3893

Show file tree
Hide file tree
Changes from 13 commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
808e94f
Tutorial 04: Spell check and grammar check
jngrad Sep 11, 2020
bed58d9
python: Document LB parameters
jngrad Sep 11, 2020
c1aed23
tutorial 04: part 1: Update text
jngrad Sep 11, 2020
4dfbd96
tutorial 04: part 2: Simplify code
jngrad Sep 11, 2020
a9a94cf
tutorial 04: part 2: Provide solutions
jngrad Sep 11, 2020
32c410f
docs: Fix documentation
jngrad Sep 19, 2020
216a2fe
tutorial 04: part 3: Move solution to the notebook
jngrad Sep 19, 2020
20882e4
tutorial 04: part 3: Provide solutions
jngrad Sep 20, 2020
4595490
Merge branch 'python' into tutorial_04_summer_school_2020
jngrad Sep 26, 2020
7094cb0
Merge branch 'python' into tutorial_04_summer_school_2020
jngrad Oct 1, 2020
15956f9
tutorial 04: part 3: Proper error analysis
jngrad Oct 4, 2020
7d43fd6
tutorial 04: part 3: Green-Kubo diffusion coefficient
jngrad Oct 6, 2020
02e102f
tutorial 04: Add tutor notes
jngrad Oct 6, 2020
f7cd6ab
tutorial 04: part 4: Switch exercise
jngrad Oct 7, 2020
59984d1
tutorial 04: Create solution cells
jngrad Oct 7, 2020
b71581d
tutorial 04: part 1: Update text, add links
jngrad Oct 7, 2020
c7d42c8
tutorial 04: Polymer diffusion with implicit solvent
jngrad Oct 7, 2020
8469df0
tutorial 04: Remove mentions of LB
jngrad Oct 7, 2020
446fd10
Merge branch 'python' into tutorial_04_summer_school_2020
RudolfWeeber Oct 8, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 14 additions & 15 deletions doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,12 @@
"#### Before you start:\n",
jngrad marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess the equation should read 

\begin{align}
 \rho &= \sum_i f_i \
 \vec{j} = \rho \\
 \vec{u} &= \sum_i f_i \vec{c_i} \\
 \Pi^{\alpha \beta} &= \sum_i f_i \vec{c_i}^{\alpha}\vec{c_i}^{\beta}
 \label{eq:fields}
\end{align}

for better readability?


Reply via ReviewNB

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ReviewNB's support for LaTeX seems incomplete, the original equation looks good in Jupyter, with one equation per line

jngrad marked this conversation as resolved.
Show resolved Hide resolved
jngrad marked this conversation as resolved.
Show resolved Hide resolved
jngrad marked this conversation as resolved.
Show resolved Hide resolved
jngrad marked this conversation as resolved.
Show resolved Hide resolved
"\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."
Expand Down Expand Up @@ -46,7 +47,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",
Expand All @@ -55,6 +56,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",
Expand Down Expand Up @@ -96,7 +98,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",
Expand Down Expand Up @@ -149,7 +151,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",
Expand All @@ -171,9 +173,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",
Expand All @@ -197,7 +197,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",
Expand Down Expand Up @@ -268,6 +268,8 @@
"+ <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",
Expand Down Expand Up @@ -314,14 +316,11 @@
"### 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",
"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 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",
Expand All @@ -332,7 +331,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."
]
Expand Down
177 changes: 152 additions & 25 deletions doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,15 +29,15 @@
"and the diffusion coefficient of a long polymer is proportional to its inverse $R_h$.\n",
"\n",
"For shorter polymers there is a transition region. It can be described\n",
"by the Kirkwood-Zimm model:\n",
"by the Kirkwood–Zimm model:\n",
"\n",
"\\begin{equation}\n",
" D=\\frac{D_0}{N} + \\frac{k_B T}{6 \\pi \\eta } \\left\\langle \\frac{1}{R_h} \\right\\rangle\n",
"\\end{equation}\n",
"\n",
"Here $D_0$ is the monomer diffusion coefficient and $\\eta$ the \n",
"viscosity of the fluid. For a finite system size the second part of the\n",
"diffusion is subject of a $1/L$ finite size effect, because\n",
"diffusion is subject to a $1/L$ finite size effect, because\n",
"hydrodynamic interactions are proportional to the inverse\n",
"distance and thus long ranged. It can be taken into account\n",
"by a correction:\n",
Expand All @@ -51,7 +51,7 @@
"too disappointed if you don't manage to do so.\n",
"\n",
"We want to determine the long-time self diffusion coefficient from the mean square\n",
"displacement of the center-of-mass a single polymer. For large $t$ the mean square displacement is\n",
"displacement of the center-of-mass of a single polymer. For large $t$ the mean square displacement is\n",
"proportional to the time and the diffusion coefficient occurs as a \n",
"prefactor:\n",
"\n",
Expand Down Expand Up @@ -98,28 +98,25 @@
"\n",
"logging.basicConfig(level=logging.INFO, stream=sys.stdout)\n",
"\n",
"\n",
"# Constants\n",
"LOOPS = 40000\n",
"STEPS = 10\n",
"STEPS = 400000\n",
"\n",
"# System setup\n",
"system = espressomd.System(box_l=[16] * 3)\n",
"system.time_step = 0.01\n",
"system.cell_system.skin = 0.4\n",
"\n",
"\n",
"lbf = espressomd.lb.LBFluidGPU(kT=1, seed=132, agrid=1, dens=1, visc=5, tau=0.01)\n",
"system.actors.add(lbf)\n",
"\n",
"\n",
"system.part.add(pos=[0, 0, 0])\n",
"\n",
"# Setup observable\n",
"pos = espressomd.observables.ParticlePositions(ids=(0,))\n",
"\n",
"# Run for different friction constants\n",
"# Run for different friction coefficients\n",
"lb_gammas = [1.0, 2.0, 4.0, 10.0]\n",
"tau_results = []\n",
"msd_results = []\n",
"\n",
"for gamma in lb_gammas:\n",
Expand All @@ -131,17 +128,16 @@
" system.integrator.run(1000)\n",
" logging.info(\"Equilibration finished.\")\n",
" # Setup observable correlator\n",
" correlator = espressomd.accumulators.Correlator(obs1=pos, tau_lin=16, tau_max=LOOPS * STEPS, delta_N=1,\n",
" correlator = espressomd.accumulators.Correlator(\n",
" obs1=pos, tau_lin=16, tau_max=STEPS, delta_N=1,\n",
" corr_operation=\"square_distance_componentwise\", compress1=\"discard1\")\n",
" system.auto_update_accumulators.add(correlator)\n",
"\n",
" logging.info(\"Sampling started for gamma = {}.\".format(gamma))\n",
" for i in range(LOOPS):\n",
" system.integrator.run(STEPS)\n",
" system.integrator.run(STEPS)\n",
" correlator.finalize()\n",
" msd_results.append(np.column_stack((correlator.lag_times(),\n",
" correlator.sample_sizes(),\n",
" correlator.result().reshape([-1, 3]))))\n",
" tau_results.append(correlator.lag_times())\n",
" msd_results.append(np.sum(correlator.result().reshape([-1, 3]), axis=1))\n",
"\n",
"logging.info(\"Sampling finished.\")"
]
Expand All @@ -164,16 +160,14 @@
"\n",
"plt.rcParams.update({'font.size': 22})\n",
"\n",
"plt.figure(figsize=(10,10))\n",
"plt.xlabel('tau [$\\Delta t$]')\n",
"plt.ylabel('msd [$\\sigma^2$]')\n",
"for index, r in enumerate(msd_results):\n",
" # adding up componentwise MSDs\n",
"plt.figure(figsize=(10, 10))\n",
"plt.xlabel(r'$\\tau$ [$\\Delta t$]')\n",
"plt.ylabel(r'MSD [$\\sigma^2$]')\n",
"for index, (tau, msd) in enumerate(zip(tau_results, msd_results)):\n",
" # We skip the first entry since it's zero by definition and cannot be displayed\n",
" # in a loglog plot. Furthermore, we only look at the first 100 entries due to\n",
" # lack of good statistics for larger times.\n",
" msd = r[1:100, 2] + r[1:100, 3] + r[1:100, 4]\n",
" plt.loglog(r[1:100, 0], msd, label=r\"$\\gamma=${}\".format(str(lb_gammas[index])))\n",
" # the high variance for larger lag times.\n",
" plt.loglog(tau[1:100], msd[1:100], label=r\"$\\gamma=${:.1f}\".format(lb_gammas[index]))\n",
"plt.legend()\n",
"plt.show()"
]
Expand All @@ -187,7 +181,138 @@
"random forces on the particle and within the LB fluid will cause the particle to move.\n",
"The mean squared displacement is calculated during the simulation via a multiple-tau\n",
"correlator. \n",
"Can you give an explanation for the quadratic time dependency for short times? Use the function <tt>curve_fit</tt> from the module <tt>scipy.optimize</tt> to produce a fit for the linear regime and determine the diffusion coefficients for the different $\\gamma$s. Calculate the diffusion coefficient for all cases and plot them as a function of γ. What relation do you observe?"
"Can you give an explanation for the quadratic time dependency for short times?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The MSD of a Brownian motion can be decomposed in three main regimes [8]:\n",
"* for short lag times $\\tau < \\tau_p$, the particle motion is not\n",
" significantly impeded by solvent collisions: it's in the ballistic mode\n",
" (collision-free regime) where $\\operatorname{MSD}(t) \\sim (k_BT / m) t^2$\n",
"* for long lag times $\\tau > \\tau_f$, the particle motion is determined by\n",
" numerous collisions with the solvent: it's in the diffusive mode where\n",
" $\\operatorname{MSD}(t) \\sim 6t$\n",
"* for lag times between $\\tau_p$ and $\\tau_f$, there is a crossover mode\n",
"\n",
"The values $\\tau_p$ and $\\tau_f$ can be obtained manually through visual\n",
"inspection of the MSD plot, or more accurately by non-linear fitting [9].\n",
"\n",
"The cutoff lag time $\\tau_p$ between the ballistic and crossover modes is proportional\n",
"to the particle mass and inversely proportional to the friction coefficient. \n",
"In the graph below, a parabola is fitted to the data points in the ballistic mode for\n",
"each $\\gamma$ and plotted beyond the crossover region to reveal the deviation from the\n",
"ballistic mode. This deviation is clearly visible in the $\\gamma = 10$ case, because\n",
"the assumption of a collision-free regime quickly breaks down when a particle is\n",
"coupled to its surrounding fluid with a high friction coefficient."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import scipy.optimize\n",
"\n",
"def quadratic(x, a, b, c):\n",
" return a * x**2 + b * x + c\n",
"\n",
"# cutoffs for the ballistic regime (different for each gamma value)\n",
"tau_p_values = [14, 12, 10, 7]\n",
"\n",
"plt.figure(figsize=(10, 10))\n",
"plt.xlabel(r'$\\tau$ [$\\Delta t$]')\n",
"plt.ylabel(r'MSD [$\\sigma^2$]')\n",
"for index, (tau_p, tau, msd) in enumerate(zip(tau_p_values, tau_results, msd_results)):\n",
" (a, b, c), _ = scipy.optimize.curve_fit(quadratic, tau[:tau_p], msd[:tau_p])\n",
" x = np.linspace(tau[0], tau[max(tau_p_values) - 1], 50)\n",
" p = plt.plot(x, quadratic(x, a, b, c), '-')\n",
" plt.plot(tau[:max(tau_p_values)], msd[:max(tau_p_values)], 'o', color=p[0].get_color(),\n",
" label=r'$\\gamma=${:.1f}'.format(lb_gammas[index]))\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the function [<tt>curve_fit()</tt>](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) from the module <tt>scipy.optimize</tt> to produce a fit for the linear regime and determine the diffusion coefficients for the different $\\gamma$s."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For large $t$ the diffusion coefficient can be expressed as:\n",
"\n",
"$$6D = \\lim_{t\\to\\infty} \\frac{\\partial \\operatorname{MSD}(t)}{\\partial t}$$\n",
"\n",
"which is simply the slope of the MSD in the diffusive mode."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def linear(x, a, b):\n",
" return a * x + b\n",
"\n",
"# cutoffs for the diffusive regime (different for each gamma value)\n",
"tau_f_values = [24, 22, 20, 17]\n",
"# cutoff for the data series (larger lag times have larger variance due to undersampling)\n",
"cutoff_limit = 90\n",
"\n",
"diffusion_results = []\n",
"\n",
"plt.figure(figsize=(10, 8))\n",
"plt.xlabel(r'$\\tau$ [$\\Delta t$]')\n",
"plt.ylabel(r'MSD [$\\sigma^2$]')\n",
"for index, (tau_f, tau, msd) in enumerate(zip(tau_f_values, tau_results, msd_results)):\n",
" (a, b), _ = scipy.optimize.curve_fit(linear, tau[tau_f:cutoff_limit], msd[tau_f:cutoff_limit])\n",
" x = np.linspace(tau[tau_f], tau[cutoff_limit - 1], 50)\n",
" p = plt.plot(x, linear(x, a, b), '-')\n",
" plt.plot(tau[tau_f:cutoff_limit], msd[tau_f:cutoff_limit], 'o', color=p[0].get_color(),\n",
" label=r'$\\gamma=${:.1f}'.format(lb_gammas[index]))\n",
" diffusion_results.append(a / 6)\n",
"\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculate the diffusion coefficient for all cases and plot them as a function of $\\gamma$. What relation do you observe?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the diffusive mode, one can derive $D = k_BT / \\gamma$ from the Stokes–Einstein relation [8]."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(10, 8))\n",
"plt.xlabel(r'$\\gamma$')\n",
"plt.ylabel('Diffusion coefficient [$\\sigma^2/t$]')\n",
"gammas = np.linspace(0.9 * min(lb_gammas), 1.1 * max(lb_gammas), 50)\n",
"plt.plot(gammas, lbf.kT / gammas, '-', label=r'$k_BT\\gamma^{-1}$')\n",
"plt.plot(lb_gammas, diffusion_results, 'o', label='D')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
Expand All @@ -202,7 +327,9 @@
"[4] P. G. de Gennes. *Scaling Concepts in Polymer Physics*. Cornell University Press, Ithaca, NY, 1979. \n",
"[5] M. Doi. *Introduction to Polymer Physics.* Clarendon Press, Oxford, 1996. \n",
"[6] Michael Rubinstein and Ralph H. Colby. *Polymer Physics.* Oxford University Press, Oxford, UK, 2003. \n",
"[7] Daan Frenkel and Berend Smit. *Understanding Molecular Simulation.* Academic Press, San Diego, second edition, 2002."
"[7] Daan Frenkel and Berend Smit. *Understanding Molecular Simulation*, section 4.4.1. Academic Press, San Diego, second edition, 2002. \n",
"[8] R. Huang, I. Chavez, K. Taute, et al. Direct observation of the full transition from ballistic to diffusive Brownian motion in a liquid. *Nature Phys.*, 7, 2011. doi:[10.1038/nphys1953](https://doi.org/10.1038/nphys1953) \n",
"[9] M. K. Riahi, I. A. Qattan, J. Hassan, D. Homouz, Identifying short- and long-time modes of the mean-square displacement: An improved nonlinear fitting approach. *AIP Advances*, 9:055112, 2019. doi:[10.1063/1.5098051](https://doi.org/10.1063/1.5098051) "
]
}
],
Expand Down
Loading