From d4fd0709ee0fcfb143ac71320251d984d77f8725 Mon Sep 17 00:00:00 2001 From: Patrick Kreissl Date: Tue, 13 Sep 2022 00:52:11 +0200 Subject: [PATCH 1/3] Add draft for sedimentation tutorial. --- .../lattice_boltzmann_sedimentation.ipynb | 396 ++++++++++++++++++ 1 file changed, 396 insertions(+) create mode 100644 doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb diff --git a/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb new file mode 100644 index 00000000000..08b40eb5ee3 --- /dev/null +++ b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb @@ -0,0 +1,396 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sedimentation in a Fluid" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The purpose of this tutorial is to demonstrate how hydrodynamic interactions can have a dramatic impact on the overall dynamics of a molecular dynamics (MD) system.\n", + "\n", + "We will set up a simple semi-two-dimensional system of sedimenting particles and simulate it first with Langevin dynamics.\n", + "Susequently, the system will be coupled to a lattice-Boltzmann fluid instead.\n", + "\n", + "Both scenarios will then be compared by visualizing the data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import espressomd\n", + "import espressomd.lb\n", + "import espressomd.lbboundaries\n", + "import espressomd.shapes\n", + "\n", + "# imports for data handling, plotting, and progress bar\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from tqdm import tqdm" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We intend to use a purely repulsive WCA interaction potential for the particles, with the following parameters:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "lj_sigma = 1\n", + "lj_epsilon = 1\n", + "lj_cutoff = 2**(1. / 6) * lj_sigma" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The intiial particles positions will be chosen to form a two-dimensional hexagonal Bravais lattice structure.\n", + "The spacing we use is large enough for the particles not to interact initially. If particles are positioned closer particles would rearrange to minimize the interaction energy, breaking up the initial structure.\n", + "We setup the simulation box such that the particle structure is periodically repeated in $x$-direction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def hexagonal_lattice(n_rows=5, parts_per_row=5, spacing=1):\n", + " positions = np.array([[x + 0.12 + (0.5 if y % 2 == 1 else 0),\n", + " y * np.sqrt(3) / 2, 0]\n", + " for x in range(parts_per_row)\n", + " for y in range(n_rows)]) * spacing\n", + " return positions" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, we are ready to define the system parameters and initialize the simulation system.\n", + "\n", + "**possible task: set up WCA potential**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# parameters\n", + "\n", + "# lattice spacing and number of lattice rows to use\n", + "spacing = lj_cutoff\n", + "n_rows = 10\n", + "\n", + "# system size in units of lattice spacing\n", + "n_height = 40\n", + "n_width = 20\n", + "n_depth = 2\n", + "\n", + "# resulting box geometry\n", + "box_height = n_height * spacing\n", + "box_width = n_width * spacing\n", + "box_depth = n_depth * spacing\n", + "\n", + "# system setup\n", + "system = espressomd.System(box_l=[box_width, box_height, box_depth])\n", + "system.time_step = 0.01\n", + "system.cell_system.skin = 0.4\n", + "\n", + "# add non-bonded WCA interaction\n", + "system.non_bonded_inter[0, 0].lennard_jones.set_params(\n", + " epsilon=lj_epsilon, sigma=lj_sigma, cutoff=lj_cutoff, shift=\"auto\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We add a wall constraint on bottom and top of the simulation box, respectively.\n", + "\n", + "**possible task: set up WCA wall constaint (wall shapes provided)**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create wall shapes bottom (b) and top (t)\n", + "wall_shape_b = espressomd.shapes.Wall(normal=[0, 1, 0], dist=1)\n", + "wall_shape_t = espressomd.shapes.Wall(\n", + " normal=[0, -1, 0], dist=-(box_height - 1))\n", + "\n", + "# add wall constraints\n", + "for wall_shape in [wall_shape_b, wall_shape_t]:\n", + " system.constraints.add(shape=wall_shape, particle_type=0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# prepare for sampling\n", + "sampling_steps = 700\n", + "\n", + "initial_positions = hexagonal_lattice(n_rows, n_width, spacing)\n", + "# total number of particles\n", + "n_parts = n_rows * n_width\n", + "assert initial_positions.shape == (n_parts, 3)\n", + "\n", + "# we introduce a small imperfection into the initial lattice structure\n", + "initial_positions[-1, 0] -= 0.3 * spacing\n", + "\n", + "# shift initial positions to the top\n", + "y_max = np.amax(initial_positions[:, 1])\n", + "initial_positions += (box_height - 1 - lj_cutoff) - y_max\n", + "\n", + "# data structures to hold particle trajectories\n", + "data_lv = np.full((sampling_steps, n_parts, 3), np.nan)\n", + "data_lb = np.full_like(data_lv, np.nan)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will set an external force acting on all particles. You can take this to model the effect of gravity, for example." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "f_gravity = [0, -3, 0]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, we start with sampling using the Langevin thermostat (with temperature $T=0$).\n", + "\n", + "**possible tasks: set Langevin thermostat, add particles**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Langevin simulation\n", + "system.thermostat.set_langevin(kT=0., gamma=15, seed=12)\n", + "\n", + "parts = system.part.add(pos=initial_positions, ext_force=[f_gravity] * n_parts)\n", + "\n", + "system.integrator.run(0)\n", + "\n", + "for step in tqdm(range(sampling_steps)):\n", + " data_lv[step, :] = parts.pos_folded \n", + " system.integrator.run(25)\n", + "data_lv[-1, :] = parts.pos_folded\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we want to sample the same system, but with a coupled lattice-Boltzmann fluid. We first reset the particles to their initial positions and remove any particle velocities.\n", + "Then we set up the LB fluid. The wall constraints that were previously added have to be also registered as LB boundaries.\n", + "\n", + "**possible tasks: set up LB fluid, add lbboundaries**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# LB simulation cleanup and setup\n", + "parts.pos = initial_positions\n", + "parts.v = [0, 0, 0]\n", + "system.thermostat.turn_off()\n", + "\n", + "lbf = espressomd.lb.LBFluid(agrid=spacing, dens=1.0, visc=1.0, tau=system.time_step, kT=0)\n", + "system.actors.add(lbf)\n", + "system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=15)\n", + "\n", + "# add LB boundaries at walls\n", + "for wall_shape in [wall_shape_b, wall_shape_t]:\n", + " no_slip_wall = espressomd.lbboundaries.LBBoundary(\n", + " shape=wall_shape, velocity=[0, 0, 0])\n", + " system.lbboundaries.add(no_slip_wall)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# numpy array to hold flowfield data\n", + "data_flowfield = np.full((sampling_steps, n_height, n_width, 3), np.nan)\n", + "\n", + "# return n_x x n_x array of lattice velocities, averaged in z-direction \n", + "def mean_flowfield(xs=n_width, ys=n_height, zs=n_depth):\n", + " dataframe = np.array([[[lbf[x, y, z].velocity for z in range(zs)] for x in range(xs)] \n", + " for y in range(ys)])\n", + " return np.mean(dataframe, axis=2)\n", + "\n", + "for step in tqdm(range(sampling_steps)):\n", + " data_lb[step, :] = parts.pos_folded\n", + " data_flowfield[step, :] = mean_flowfield() \n", + " system.integrator.run(25)\n", + "data_lb[-1, :] = parts.pos_folded\n", + "data_flowfield[-1, :] = mean_flowfield()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's visualize the data. First some imports and definitions for inline visualization." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.animation as animation\n", + "import tempfile\n", + "import base64\n", + "\n", + "from matplotlib.quiver import Quiver\n", + "\n", + "VIDEO_TAG = \"\"\"\"\"\"\n", + "\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", + "\n", + "animation.Animation._repr_html_ = anim_to_html\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')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And now the actual visualization code." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# setup figure and prepare axes\n", + "fig = plt.figure(figsize=(2 * 5, 5 / box_width * box_height))\n", + "gs = fig.add_gridspec(1, 2, wspace=0.1)\n", + "(ax1, ax2) = gs.subplots(sharey=True)\n", + "\n", + "ax1.set_title(\"Langevin\")\n", + "ax1.set_xlim((0, box_width))\n", + "ax1.set_ylim((0, box_height))\n", + "\n", + "ax2.set_title(\"LB Fluid\")\n", + "ax2.set_xlim((0, box_width))\n", + "ax2.set_ylim((0, box_height))\n", + "\n", + "# draw walls\n", + "for ax in [ax1, ax2]:\n", + " ax.hlines((1, box_height-1), 0, box_width, color=\"gray\")\n", + "\n", + "# create meshgrid for quiver plot\n", + "xs = np.array([x for x in range(n_width)]) * spacing\n", + "ys = np.array([y for y in range(n_height)]) * spacing\n", + "X, Y = np.meshgrid(xs, ys)\n", + "\n", + "# initialize plot objects\n", + "lb_ff = ax2.quiver(X, Y, data_flowfield[0, :, :, 0],\n", + " data_flowfield[0, :, :, 1])\n", + "lb_particles, = ax2.plot([], [], 'o')\n", + "lv_particles, = ax1.plot([], [], 'o')\n", + "\n", + "def draw_frame(t):\n", + " # manually remove Quivers from ax2\n", + " for artist in ax2.get_children():\n", + " if isinstance(artist, Quiver):\n", + " artist.remove()\n", + " \n", + " # draw new quivers\n", + " lb_ff = ax2.quiver(X, Y, data_flowfield[t, :, :, 0],\n", + " data_flowfield[t, :, :, 1])\n", + " \n", + " # draw particles\n", + " lv_particles.set_data(data_lv[t, :, 0], data_lv[t, :, 1])\n", + " lb_particles.set_data(data_lb[t, :, 0], data_lb[t, :, 1])\n", + "\n", + " return [lv_particles, lb_particles, lb_ff]\n", + "\n", + "animation.FuncAnimation(fig, draw_frame, frames=sampling_steps, blit=True, interval=0, repeat=False)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.12" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} From bdb02c70099f5398cf52a20f016f26d79255591e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Sat, 1 Oct 2022 18:18:53 +0200 Subject: [PATCH 2/3] doc: Finish LB sedimentation tutorial Use normalized quiver plots. Speed-up sampling using observables. Introduce crystalline defects to the sides of the box to observe the drop at the center of the box. Explain LD and LB thermostats. --- doc/sphinx/integration.rst | 1 + doc/tutorials/Readme.md | 3 +- .../lattice_boltzmann/CMakeLists.txt | 5 +- .../lattice_boltzmann/NotesForTutor.md | 16 + .../lattice_boltzmann_sedimentation.ipynb | 443 +++++++++++++----- 5 files changed, 356 insertions(+), 112 deletions(-) diff --git a/doc/sphinx/integration.rst b/doc/sphinx/integration.rst index fac5beef1b0..bb759c377f4 100644 --- a/doc/sphinx/integration.rst +++ b/doc/sphinx/integration.rst @@ -440,6 +440,7 @@ Best explained in an example:: As explained before the temperature is set as thermal energy :math:`k_\mathrm{B} T`. The Langevin thermostat is based on an extension of Newton's equation of motion to +account for drag and collisions with a fluid: .. math:: m_i \dot{v}_i(t) = f_i(\{x_j\},v_i,t) - \gamma v_i(t) + \sqrt{2\gamma k_B T} \eta_i(t). diff --git a/doc/tutorials/Readme.md b/doc/tutorials/Readme.md index 2f4eb711b49..fd3d5fd307b 100644 --- a/doc/tutorials/Readme.md +++ b/doc/tutorials/Readme.md @@ -39,7 +39,8 @@ physical systems. Simulations including hydrodynamic interactions using the Lattice-Boltzmann method. Guide [Part 1](lattice_boltzmann/lattice_boltzmann_theory.ipynb) | - [Part 2](lattice_boltzmann/lattice_boltzmann_poiseuille_flow.ipynb) + [Part 2](lattice_boltzmann/lattice_boltzmann_poiseuille_flow.ipynb) | + [Part 3](lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb) * **Polymers** Modelling polymers with hydrodynamic interactions. [Guide](polymers/polymers.ipynb) diff --git a/doc/tutorials/lattice_boltzmann/CMakeLists.txt b/doc/tutorials/lattice_boltzmann/CMakeLists.txt index 9eea954b7be..b47b538d82e 100644 --- a/doc/tutorials/lattice_boltzmann/CMakeLists.txt +++ b/doc/tutorials/lattice_boltzmann/CMakeLists.txt @@ -1,9 +1,12 @@ configure_tutorial_target( TARGET tutorial_lb DEPENDS lattice_boltzmann_theory.ipynb - lattice_boltzmann_poiseuille_flow.ipynb figures/latticeboltzmann-grid.png + lattice_boltzmann_poiseuille_flow.ipynb lattice_boltzmann_sedimentation.ipynb + figures/latticeboltzmann-grid.png figures/latticeboltzmann-momentumexchange.png) nb_export(TARGET tutorial_lb SUFFIX "1" FILE "lattice_boltzmann_theory.ipynb" HTML_RUN) nb_export(TARGET tutorial_lb SUFFIX "2" FILE "lattice_boltzmann_poiseuille_flow.ipynb" HTML_RUN) +nb_export(TARGET tutorial_lb SUFFIX "3" FILE + "lattice_boltzmann_sedimentation.ipynb" HTML_RUN) diff --git a/doc/tutorials/lattice_boltzmann/NotesForTutor.md b/doc/tutorials/lattice_boltzmann/NotesForTutor.md index fdb93fde28e..18b1a0e6f3a 100644 --- a/doc/tutorials/lattice_boltzmann/NotesForTutor.md +++ b/doc/tutorials/lattice_boltzmann/NotesForTutor.md @@ -26,3 +26,19 @@ In the course of this tutorial, students should learn to: * Set up an LB fluid with LB boundaries * Retrieve data from the LB grid + +# Part 3: Sedimentation + +## Physics learning goals + +After the tutorial, students should be able to: + +* Describe hydrodynamic effects on particle diffusion + +## ESPResSo learning goals + +In the course of this tutorial, students should learn to: + +* Set up external forces on particles (e.g. gravity, constant electric field) +* Set up particle observables and LB profile observables +* Register observables in accumulators diff --git a/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb index 08b40eb5ee3..8c479f64930 100644 --- a/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb +++ b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb @@ -4,19 +4,64 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Sedimentation in a Fluid" + "# Sedimentation in a fluid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The purpose of this tutorial is to demonstrate how hydrodynamic interactions can have a dramatic impact on the overall dynamics of a molecular dynamics (MD) system.\n", - "\n", - "We will set up a simple semi-two-dimensional system of sedimenting particles and simulate it first with Langevin dynamics.\n", - "Susequently, the system will be coupled to a lattice-Boltzmann fluid instead.\n", - "\n", - "Both scenarios will then be compared by visualizing the data." + "The purpose of this tutorial is to demonstrate how hydrodynamic interactions can have a dramatic impact\n", + "on the overall dynamics of a molecular system. We will set up a simple semi-two-dimensional system of\n", + "sedimenting particles and simulate their interaction with a solvent.\n", + "\n", + "In the first scenario, Langevin dynamics (LD) will model the friction from a \"simple\" implicit solvent.\n", + "Although LD in general doesn't fully model an implicit solvent due to the lack of contribution\n", + "for electrostatic screening and the hydrophobic effect, these limitations won't affect this particular\n", + "simulation, since we represent the sediment particles as electrically neutral and fully solvated beads.\n", + "\n", + "LD adds two extra terms to Newton's second law, which account for the friction and\n", + "random collisions from a solvent:\n", + "\n", + "\\begin{equation}\n", + "m_i \\dot{v}_i(t) = f_i(t) - \\gamma v_i(t) + \\sqrt{2\\gamma k_{\\mathrm{B}} T} \\cdot \\eta_i(t)\n", + "\\end{equation}\n", + "\n", + "with $m_i$, $v_i(t)$, $f_i(t)$ the mass, velocity and forces acting on particle $i$ at time $t$,\n", + "$\\gamma$ the dampening constant, $k_{\\mathrm{B}}$ the Boltzmann constant, $T$ the temperature,\n", + "$\\eta_i(t)$ a random uniform number drawn for particle $i$ at time $t$.\n", + "The second term in the right-hand side is responsible for drag due to friction with the solvent,\n", + "while the third term is responsible for heat exchange due to collisions with the solvent.\n", + "This method doesn't capture hydrodynamic effects, since the particles only interact with each\n", + "other via their inter-atomic potentials; when two particles are further away than their\n", + "potential cutoff, they cannot interact.\n", + "\n", + "In the second scenario, we will use the lattice-Boltzmann (LB) method to model the solvent using\n", + "a grid-based solver for the Navier-Stokes equation in the limit of a low Reynolds number.\n", + "Particle–fluid interaction is achieved through momentum exchange via frictional coupling:\n", + "\n", + "\\begin{equation}\n", + "m_i \\dot{v}_i(t) = f_i(t) - \\gamma \\left( v_i(t) - u(x_i(t), t) \\right) + \\sqrt{2\\gamma k_{\\mathrm{B}} T} \\cdot \\eta_i(t)\n", + "\\end{equation}\n", + "\n", + "with $u(x_i(t), t)$ the interpolated fluid velocity at the position of the particle.\n", + "The frictional term is also back-propagated to the nearest nodes of the LB fluid.\n", + "While computationally more demanding, this method recovers hydrodynamic effects, i.e. particles\n", + "in close proximity influence each other due to coupling with the fluid inbetween them, even when\n", + "their separation distance exceeds the inter-atomic potential cutoff. In addition, the fluid can\n", + "develop a macroscopic flow field, in which case all particles will influence each other over\n", + "large distances.\n", + "\n", + "We will compare both scenarios by generating a video of the sedimentation dynamics.\n", + "To make the effect more clearly visible, we will not couple the particles and fluid to\n", + "a thermal bath." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## System setup" ] }, { @@ -29,37 +74,23 @@ "import espressomd.lb\n", "import espressomd.lbboundaries\n", "import espressomd.shapes\n", + "import espressomd.observables\n", + "import espressomd.accumulators\n", "\n", "# imports for data handling, plotting, and progress bar\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", - "from tqdm import tqdm" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We intend to use a purely repulsive WCA interaction potential for the particles, with the following parameters:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "lj_sigma = 1\n", - "lj_epsilon = 1\n", - "lj_cutoff = 2**(1. / 6) * lj_sigma" + "import tqdm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The intiial particles positions will be chosen to form a two-dimensional hexagonal Bravais lattice structure.\n", - "The spacing we use is large enough for the particles not to interact initially. If particles are positioned closer particles would rearrange to minimize the interaction energy, breaking up the initial structure.\n", + "The initial particles positions will be chosen to form a two-dimensional hexagonal Bravais lattice structure.\n", + "The spacing we use is large enough for the particles not to interact initially. If particles were positioned\n", + "closer, they would rearrange to minimize the interaction energy, breaking up the initial structure.\n", + "\n", "We setup the simulation box such that the particle structure is periodically repeated in $x$-direction." ] }, @@ -69,7 +100,7 @@ "metadata": {}, "outputs": [], "source": [ - "def hexagonal_lattice(n_rows=5, parts_per_row=5, spacing=1):\n", + "def hexagonal_lattice(n_rows, parts_per_row, spacing):\n", " positions = np.array([[x + 0.12 + (0.5 if y % 2 == 1 else 0),\n", " y * np.sqrt(3) / 2, 0]\n", " for x in range(parts_per_row)\n", @@ -81,9 +112,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now, we are ready to define the system parameters and initialize the simulation system.\n", - "\n", - "**possible task: set up WCA potential**" + "Now, we are ready to define the system parameters and initialize the simulation system." ] }, { @@ -92,7 +121,10 @@ "metadata": {}, "outputs": [], "source": [ - "# parameters\n", + "# WCA potential parameters\n", + "lj_sigma = 1.\n", + "lj_epsilon = 1.\n", + "lj_cutoff = 2**(1. / 6.) * lj_sigma\n", "\n", "# lattice spacing and number of lattice rows to use\n", "spacing = lj_cutoff\n", @@ -120,19 +152,28 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, "source": [ "We add a wall constraint on bottom and top of the simulation box, respectively.\n", "\n", - "**possible task: set up WCA wall constaint (wall shapes provided)**" + "**Exercise:**\n", + "* set up two wall constraints ``wall_shape_b`` (bottom) and ``wall_shape_t`` (top) in the *xz*-plane\n", + " at a distance of 1 unit from the top and bottom sides of the box\n", + " ([user guide](https://espressomd.github.io/doc/espressomd.html#espressomd.shapes.Wall))\n", + "* add these two shapes to the system constraints using particle type 0\n", + " ([user guide](https://espressomd.github.io/doc/constraints.html#adding-shape-based-constraints-to-the-system))" ] }, { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, "source": [ + "```python\n", "# create wall shapes bottom (b) and top (t)\n", "wall_shape_b = espressomd.shapes.Wall(normal=[0, 1, 0], dist=1)\n", "wall_shape_t = espressomd.shapes.Wall(\n", @@ -140,7 +181,8 @@ "\n", "# add wall constraints\n", "for wall_shape in [wall_shape_b, wall_shape_t]:\n", - " system.constraints.add(shape=wall_shape, particle_type=0)" + " system.constraints.add(shape=wall_shape, particle_type=0)\n", + "```" ] }, { @@ -148,32 +190,72 @@ "execution_count": null, "metadata": {}, "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, "source": [ - "# prepare for sampling\n", + "We will now calculate the particle initial positions and introduce a small crystalline defect to\n", + "help break the symmetry of the system.\n", + "We will also configure an external force acting on all particles, which models the effect of gravity." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# gravitational force\n", + "f_gravity = [0., -3., 0.]\n", + "\n", + "# number of frames in the output trajectory\n", "sampling_steps = 700\n", "\n", + "# lattice positions\n", "initial_positions = hexagonal_lattice(n_rows, n_width, spacing)\n", - "# total number of particles\n", - "n_parts = n_rows * n_width\n", - "assert initial_positions.shape == (n_parts, 3)\n", "\n", - "# we introduce a small imperfection into the initial lattice structure\n", - "initial_positions[-1, 0] -= 0.3 * spacing\n", + "# introduce a small imperfection in the initial lattice structure\n", + "initial_positions[99, 0] += 0.15 * spacing\n", + "initial_positions[109, 0] -= 0.15 * spacing\n", "\n", "# shift initial positions to the top\n", "y_max = np.amax(initial_positions[:, 1])\n", "initial_positions += (box_height - 1 - lj_cutoff) - y_max\n", + "initial_positions = np.remainder(initial_positions, np.copy(system.box_l))\n", "\n", - "# data structures to hold particle trajectories\n", - "data_lv = np.full((sampling_steps, n_parts, 3), np.nan)\n", - "data_lb = np.full_like(data_lv, np.nan)" + "# total number of particles\n", + "n_parts = n_rows * n_width\n", + "assert initial_positions.shape == (n_parts, 3)" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, "source": [ - "We will set an external force acting on all particles. You can take this to model the effect of gravity, for example." + "## Langevin dynamics\n", + "\n", + "In this scenario, we will sample the sedimentation dynamics using the Langevin thermostat.\n", + "\n", + "**Exercise:**\n", + "* set up an unthermalized Langevin thermostat to add a purely frictional term to the equations\n", + " of motion with $\\gamma=15$\n", + " ([user guide](https://espressomd.github.io/doc/integration.html#langevin-thermostat))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "system.thermostat.set_langevin(kT=0., gamma=15., seed=12)\n", + "```" ] }, { @@ -181,17 +263,14 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "f_gravity = [0, -3, 0]" - ] + "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now, we start with sampling using the Langevin thermostat (with temperature $T=0$).\n", - "\n", - "**possible tasks: set Langevin thermostat, add particles**" + "We can now sample the particle positions as a function of time.\n", + "We will use a particle observable and an accumulator to record the trajectory." ] }, { @@ -200,27 +279,24 @@ "metadata": {}, "outputs": [], "source": [ - "# Langevin simulation\n", - "system.thermostat.set_langevin(kT=0., gamma=15, seed=12)\n", - "\n", "parts = system.part.add(pos=initial_positions, ext_force=[f_gravity] * n_parts)\n", + "obs_particle_pos = espressomd.observables.ParticlePositions(ids=list(range(n_parts)))\n", + "acc_particle_pos = espressomd.accumulators.TimeSeries(obs=obs_particle_pos, delta_N=1)\n", "\n", "system.integrator.run(0)\n", "\n", - "for step in tqdm(range(sampling_steps)):\n", - " data_lv[step, :] = parts.pos_folded \n", + "for step in tqdm.tqdm(range(sampling_steps)):\n", + " acc_particle_pos.update()\n", " system.integrator.run(25)\n", - "data_lv[-1, :] = parts.pos_folded\n" + "\n", + "data_ld = np.remainder(np.reshape(acc_particle_pos.time_series(), (-1, n_parts, 3)), system.box_l)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now we want to sample the same system, but with a coupled lattice-Boltzmann fluid. We first reset the particles to their initial positions and remove any particle velocities.\n", - "Then we set up the LB fluid. The wall constraints that were previously added have to be also registered as LB boundaries.\n", - "\n", - "**possible tasks: set up LB fluid, add lbboundaries**" + "We will now disable the thermostat, reset the particles to their initial positions and zero out particle velocities." ] }, { @@ -229,20 +305,82 @@ "metadata": {}, "outputs": [], "source": [ - "# LB simulation cleanup and setup\n", "parts.pos = initial_positions\n", "parts.v = [0, 0, 0]\n", - "system.thermostat.turn_off()\n", - "\n", - "lbf = espressomd.lb.LBFluid(agrid=spacing, dens=1.0, visc=1.0, tau=system.time_step, kT=0)\n", + "system.thermostat.turn_off()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Hydrodynamics" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "In this scenario, we want to sample the same system coupled to a lattice-Boltzmann fluid.\n", + "\n", + "**Exercise:**\n", + "* create an unthermalized lattice-Boltzmann object with viscosity 1, density 1, grid size equal to\n", + " ``spacing`` and LB time step equal to the MD time step and add this object to the system list of actors\n", + " ([user guide](https://espressomd.github.io/doc/lb.html#setting-up-a-lb-fluid))\n", + "* activate particle coupling to the fluid by setting the LB thermostat with $\\gamma=15$\n", + " ([user guide](https://espressomd.github.io/doc/lb.html#coupling-lb-to-a-md-simulation))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "lbf = espressomd.lb.LBFluid(agrid=spacing, dens=1., visc=1., tau=system.time_step, kT=0.)\n", "system.actors.add(lbf)\n", - "system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=15)\n", + "system.thermostat.set_lb(LB_fluid=lbf, gamma=15., seed=0)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "The wall constraints that were previously added now have to be registered as LB boundaries.\n", "\n", - "# add LB boundaries at walls\n", + "**Exercise:**\n", + "* convert the wall shapes to LB boundaries and add them to the system list of LB boundaries\n", + " ([user guide](https://espressomd.github.io/doc/lb.html#using-shapes-as-lattice-boltzmann-boundary))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "# add LB boundaries\n", "for wall_shape in [wall_shape_b, wall_shape_t]:\n", " no_slip_wall = espressomd.lbboundaries.LBBoundary(\n", " shape=wall_shape, velocity=[0, 0, 0])\n", - " system.lbboundaries.add(no_slip_wall)" + " system.lbboundaries.add(no_slip_wall)\n", + "```" ] }, { @@ -250,28 +388,98 @@ "execution_count": null, "metadata": {}, "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, "source": [ - "# numpy array to hold flowfield data\n", - "data_flowfield = np.full((sampling_steps, n_height, n_width, 3), np.nan)\n", - "\n", - "# return n_x x n_x array of lattice velocities, averaged in z-direction \n", - "def mean_flowfield(xs=n_width, ys=n_height, zs=n_depth):\n", - " dataframe = np.array([[[lbf[x, y, z].velocity for z in range(zs)] for x in range(xs)] \n", - " for y in range(ys)])\n", - " return np.mean(dataframe, axis=2)\n", - "\n", - "for step in tqdm(range(sampling_steps)):\n", - " data_lb[step, :] = parts.pos_folded\n", - " data_flowfield[step, :] = mean_flowfield() \n", + "We will plot the fluid flow field in the final video using 2D vectors.\n", + "To this end, we need to record the fluid trajectory with the same frequency as the particle positions.\n", + "\n", + "**Exercise:**\n", + "* create a LB velocity profile in Cartesian coordinates\n", + " ([user guide](https://espressomd.github.io/doc/espressomd.html#espressomd.observables.LBVelocityProfile))\n", + "* register that observable in a time series accumulator named ``acc_lb_vel``\n", + " ([user guide](https://espressomd.github.io/doc/analysis.html#time-series))\n", + "\n", + "**Hints:**\n", + "* the velocity observable takes parameters in MD units, not LB units\n", + "* there is no fluid inside the top an bottom boundaries, therefore the number of bins for the *y*-axis\n", + " is smaller than ``n_height`` by 2 units (parameters ``min_y`` and ``max_y`` are also affected)\n", + "* use ``n_z_bins=1`` to average the velocity along the *z*-direction\n", + "* for ``sampling_delta_*`` and ``sampling_offset_*``, use ``spacing`` and ``0.5 * spacing`` respectively" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "obs_lb_vel = espressomd.observables.LBVelocityProfile(\n", + " n_x_bins=n_width,\n", + " n_y_bins=n_height - 2, # skip data inside the LB boundaries (top and bottom walls)\n", + " n_z_bins=1, # averaged velocity along the z-direction\n", + " min_x=0.0,\n", + " min_y=spacing,\n", + " min_z=0.0,\n", + " max_x=system.box_l[0],\n", + " max_y=system.box_l[1] - spacing,\n", + " max_z=system.box_l[2],\n", + " sampling_delta_x=spacing,\n", + " sampling_delta_y=spacing,\n", + " sampling_delta_z=spacing,\n", + " sampling_offset_x=0.5 * spacing,\n", + " sampling_offset_y=0.5 * spacing,\n", + " sampling_offset_z=0.5 * spacing,\n", + " allow_empty_bins=True)\n", + "acc_lb_vel = espressomd.accumulators.TimeSeries(obs=obs_lb_vel, delta_N=1)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now sample the particle positions and fluid velocity as a function of time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "obs_particle_pos = espressomd.observables.ParticlePositions(ids=list(range(n_parts)))\n", + "acc_particle_pos = espressomd.accumulators.TimeSeries(obs=obs_particle_pos, delta_N=1)\n", + "\n", + "for step in tqdm.tqdm(range(sampling_steps)):\n", + " acc_lb_vel.update()\n", + " acc_particle_pos.update()\n", " system.integrator.run(25)\n", - "data_lb[-1, :] = parts.pos_folded\n", - "data_flowfield[-1, :] = mean_flowfield()" + "\n", + "data_lb = np.remainder(np.reshape(acc_particle_pos.time_series(), (-1, n_parts, 3)), system.box_l)\n", + "data_flowfield = acc_lb_vel.time_series()[:, :, :, 0, 0:2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "## Visualization\n", + "\n", "Now let's visualize the data. First some imports and definitions for inline visualization." ] }, @@ -282,16 +490,18 @@ "outputs": [], "source": [ "import matplotlib.animation as animation\n", + "import matplotlib.quiver\n", "import tempfile\n", "import base64\n", "\n", - "from matplotlib.quiver import Quiver\n", - "\n", "VIDEO_TAG = \"\"\"\"\"\"\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", @@ -303,31 +513,33 @@ " plt.close(anim._fig)\n", " return VIDEO_TAG.format(anim._encoded_video)\n", "\n", - "\n", - "animation.Animation._repr_html_ = anim_to_html\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')" + "animation.Animation._repr_html_ = anim_to_html" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "And now the actual visualization code." + "And now the actual visualization code.\n", + "\n", + "Note: if Jupyter encapsulates the video in a scrollable area, click on the ``Out[]:`` text to the left of\n", + "the output cell to toggle auto-scrolling off. This can also be achieved by hitting the key combination\n", + "Shift+O after highlighting the output cell." ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": false + }, "outputs": [], "source": [ "# setup figure and prepare axes\n", "fig = plt.figure(figsize=(2 * 5, 5 / box_width * box_height))\n", "gs = fig.add_gridspec(1, 2, wspace=0.1)\n", - "(ax1, ax2) = gs.subplots(sharey=True)\n", + "ax1 = plt.subplot(gs[0])\n", + "ax2 = plt.subplot(gs[1], sharey=ax1)\n", "\n", "ax1.set_title(\"Langevin\")\n", "ax1.set_xlim((0, box_width))\n", @@ -343,30 +555,41 @@ "\n", "# create meshgrid for quiver plot\n", "xs = np.array([x for x in range(n_width)]) * spacing\n", - "ys = np.array([y for y in range(n_height)]) * spacing\n", + "ys = np.array([y for y in range(1, n_height-1)]) * spacing\n", "X, Y = np.meshgrid(xs, ys)\n", "\n", + "# create a transposed flow field for quiver plot\n", + "data_flowfield_t = np.transpose(data_flowfield, axes=(0, 2, 1, 3))\n", + "\n", + "# set quiver scale (fraction of the highest velocity in the XY plane)\n", + "lb_vel_max = np.sum(np.square(data_flowfield), axis=-1)\n", + "quiver_scale = np.sqrt(np.max(lb_vel_max))\n", + "\n", + "def plot_lb_vel(ax, X, Y, flowfield, t, scale):\n", + " return ax.quiver(X, Y,\n", + " flowfield[t, :, :, 0],\n", + " flowfield[t, :, :, 1],\n", + " scale_units=\"xy\", scale=scale)\n", + "\n", "# initialize plot objects\n", - "lb_ff = ax2.quiver(X, Y, data_flowfield[0, :, :, 0],\n", - " data_flowfield[0, :, :, 1])\n", + "lb_ff = plot_lb_vel(ax2, X, Y, data_flowfield_t, 0, quiver_scale)\n", "lb_particles, = ax2.plot([], [], 'o')\n", - "lv_particles, = ax1.plot([], [], 'o')\n", + "ld_particles, = ax1.plot([], [], 'o')\n", "\n", "def draw_frame(t):\n", " # manually remove Quivers from ax2\n", " for artist in ax2.get_children():\n", - " if isinstance(artist, Quiver):\n", + " if isinstance(artist, matplotlib.quiver.Quiver):\n", " artist.remove()\n", - " \n", + "\n", " # draw new quivers\n", - " lb_ff = ax2.quiver(X, Y, data_flowfield[t, :, :, 0],\n", - " data_flowfield[t, :, :, 1])\n", - " \n", + " lb_ff = plot_lb_vel(ax2, X, Y, data_flowfield_t, t, quiver_scale)\n", + "\n", " # draw particles\n", - " lv_particles.set_data(data_lv[t, :, 0], data_lv[t, :, 1])\n", + " ld_particles.set_data(data_ld[t, :, 0], data_ld[t, :, 1])\n", " lb_particles.set_data(data_lb[t, :, 0], data_lb[t, :, 1])\n", "\n", - " return [lv_particles, lb_particles, lb_ff]\n", + " return [ld_particles, lb_particles, lb_ff]\n", "\n", "animation.FuncAnimation(fig, draw_frame, frames=sampling_steps, blit=True, interval=0, repeat=False)" ] From c276fe55074e62423d9c07a42e5fce7ac956ac2c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Sat, 1 Oct 2022 18:29:22 +0200 Subject: [PATCH 3/3] tests: Check LB sedimentation tutorial --- testsuite/scripts/tutorials/CMakeLists.txt | 1 + .../test_lattice_boltzmann_sedimentation.py | 71 +++++++++++++++++++ 2 files changed, 72 insertions(+) create mode 100644 testsuite/scripts/tutorials/test_lattice_boltzmann_sedimentation.py diff --git a/testsuite/scripts/tutorials/CMakeLists.txt b/testsuite/scripts/tutorials/CMakeLists.txt index 6ea9dfb57b8..b368a6da371 100644 --- a/testsuite/scripts/tutorials/CMakeLists.txt +++ b/testsuite/scripts/tutorials/CMakeLists.txt @@ -52,6 +52,7 @@ tutorial_test(FILE test_langevin_dynamics.py) tutorial_test(FILE test_polymers.py SUFFIX rouse) tutorial_test(FILE test_polymers.py SUFFIX zimm LABELS "gpu") tutorial_test(FILE test_lattice_boltzmann_poiseuille_flow.py LABELS "gpu") +tutorial_test(FILE test_lattice_boltzmann_sedimentation.py) tutorial_test(FILE test_raspberry_electrophoresis.py LABELS "gpu") tutorial_test(FILE test_active_matter.py LABELS "gpu") tutorial_test(FILE test_electrokinetics.py LABELS "gpu") diff --git a/testsuite/scripts/tutorials/test_lattice_boltzmann_sedimentation.py b/testsuite/scripts/tutorials/test_lattice_boltzmann_sedimentation.py new file mode 100644 index 00000000000..47390a41a76 --- /dev/null +++ b/testsuite/scripts/tutorials/test_lattice_boltzmann_sedimentation.py @@ -0,0 +1,71 @@ +# +# Copyright (C) 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 . +# + +import unittest as ut +import importlib_wrapper +import numpy as np +import scipy.stats + + +tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( + "@TUTORIALS_DIR@/lattice_boltzmann/lattice_boltzmann_sedimentation.py", + sampling_steps=400) + + +def curl2d(flow, spacing): + # curl_z = d(v_y) / dx - d(v_x) / dy + dvy_dx = np.gradient(flow[:, :, 1], spacing[0], axis=0) + dvx_dy = np.gradient(flow[:, :, 0], spacing[1], axis=1) + return dvy_dx - dvx_dy + + +def get_peak_position(mat, kernel): + return np.unravel_index(kernel(mat, axis=None), mat.shape) + + +@skipIfMissingFeatures +class Tutorial(ut.TestCase): + system = tutorial.system + + def test_flow_profile(self): + # slice trajectory to keep only the flow field onset + flow_field = tutorial.data_flowfield[200:400, :, :, 0:2] + vortices = np.zeros((2, flow_field.shape[0], 2)) + for i in range(flow_field.shape[0]): + curl = curl2d(flow_field[i], 2 * [tutorial.spacing]) + vortices[0, i] = get_peak_position(curl, np.argmax) + vortices[1, i] = get_peak_position(curl, np.argmin) + + # check flow field curl + ref_pos_x = [12.5, 7.5] # LB units + ref_pos_y = 2 * [11.] # LB units + for i in range(2): + width = tutorial.n_width # LB units + vortex_avg_x = scipy.stats.circmean(vortices[i, :, 0], high=width) + vortex_std_x = scipy.stats.circstd(vortices[i, :, 0], high=width) + vortex_avg_y = np.mean(vortices[i, :, 1]) + vortex_std_y = np.std(vortices[i, :, 1]) + self.assertAlmostEqual(vortex_avg_x, ref_pos_x[i], delta=2.) + self.assertAlmostEqual(vortex_avg_y, ref_pos_y[i], delta=2.) + self.assertLess(vortex_std_x, 4.) + self.assertLess(vortex_std_y, 6.) + + +if __name__ == "__main__": + ut.main()