From b7617a301ef0e3c84ed29a2739d1a695887e61ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 2 Oct 2023 14:44:57 +0200 Subject: [PATCH] Rework electrokinetics tutorial Write tests to check for analytical solutions and turbulent flow. Improve matplotlib figures, ffmpeg videos and png figure captions. --- .../electrokinetics/electrokinetics.ipynb | 219 ++++++++++-------- src/python/espressomd/electrokinetics.py | 28 +-- .../scripts/tutorials/test_electrokinetics.py | 93 +++++++- 3 files changed, 219 insertions(+), 121 deletions(-) diff --git a/doc/tutorials/electrokinetics/electrokinetics.ipynb b/doc/tutorials/electrokinetics/electrokinetics.ipynb index 3c34856287..f59d994f9d 100644 --- a/doc/tutorials/electrokinetics/electrokinetics.ipynb +++ b/doc/tutorials/electrokinetics/electrokinetics.ipynb @@ -28,12 +28,12 @@ "\\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 \\epsilon_{0} \\epsilon_{\\mathrm{r}}} \\sum_{i} z_{i} e n_{i} \\\\\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", - "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, $\\epsilon_{0}$ the vacuum permittivity, $\\epsilon_{\\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." + "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." ] }, { @@ -107,7 +107,7 @@ "FLUID_VELOCITY = [0.04, 0.04, 0.0]\n", "KT = 1.0\n", "\n", - "RUN_TIME = 200" + "RUN_TIME = 400" ] }, { @@ -135,8 +135,7 @@ "metadata": {}, "outputs": [], "source": [ - "lattice = espressomd.lb.LatticeWalberla(\n", - " n_ghost_layers=1, agrid=AGRID)\n", + "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", @@ -212,7 +211,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "To compare our simulation to the fundamental solution of the Advection-Diffusion equation, we need to approximate a delta-droplet, which can be done by having only a non-zero density at the center of the domain." + "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." ] }, { @@ -244,7 +243,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "For comparison, we prepare the analytical solution and show the 2d-density as well as a slice through the center of the droplet." + "For comparison, we prepare the analytical solution and show the 2D-density as well as a slice through the center of the droplet." ] }, { @@ -264,9 +263,9 @@ "pos[..., 1], pos[..., 0] = np.meshgrid(xpos, ypos)\n", "\n", "# add the advection shift\n", - "pos -= np.asarray(FLUID_VELOCITY[:2]) * TAU * RUN_TIME\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)" + "analytic_density = calc_gaussian(pos=pos, time=RUN_TIME * TAU, D=DIFFUSION_COEFFICIENT)" ] }, { @@ -282,7 +281,8 @@ "\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)" + "fig.colorbar(imshow, ax=[ax1, ax2], shrink=0.8)\n", + "plt.show()" ] }, { @@ -295,16 +295,25 @@ "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", - "\n", - "ax.plot(positions_diagonal, values_diagonal, \"x\", label=\"simulation\")\n", - "ax.plot(positions_diagonal, analytic_diagonal, label=\"analytic\")\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", - "plt.legend()" + "plt.legend()\n", + "plt.show()" ] }, { @@ -329,11 +338,11 @@ "\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", - "
\n", - " missing\n", - "
\n", - "
Figure 1: Slit pore system and coordinate system used for the analytical calculations.
\n", - "
\n", + "
\n", + "\n", + "
\n", + "
Figure 1: Slit pore system and coordinate system used for the analytical calculations.
\n", + "
\n", "
" ] }, @@ -346,47 +355,47 @@ "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}{\\epsilon_{0} \\epsilon_{\\mathrm{r}}} \\sum_{i} z_{i} e n_{i}(x) \\exp \\left( -\\frac{z_{i} e \\phi(x)}{k_{\\mathrm{B}} T} \\right)\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", - "\\begin{equation}\n", - "\\phi(x) = -\\frac{k_{B}T}{z e} \\log \\left[ \\frac{C^2 \\epsilon_{0} \\epsilon_{\\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", - "\\end{equation}\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", - "\\begin{equation}\n", - "n(x) = \\frac{C^2 \\epsilon_{0} \\epsilon_{\\mathrm{r}}}{2 k_{B}T} \\cos^{-2} \\left( \\frac{z e C}{2 k_{B} T} x \\right).\n", - "\\end{equation}\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", - "\\begin{equation}\n", - "C \\tan \\left( \\frac{z e d}{4 k_{B} T} C \\right) = - \\frac{e^2}{\\epsilon_{0} \\epsilon_{\\mathrm{r}}} \\sigma,\n", - "\\end{equation}\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", - "\\begin{equation}\n", - "\\frac{\\partial^2 v_{y}(x)}{\\partial x^2} = - \\frac{\\epsilon_{0} \\epsilon_{\\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", - "\\end{equation}\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", - "\\begin{equation}\n", - "v_{y}(x) = \\frac{2 \\epsilon_{0} \\epsilon_{\\mathrm{r}} k_{B} T E}{\\eta z e} \\log \\left( \\frac{\\cos \\left( \\frac{z e C}{2 k_{B} T} x \\right)}{\\cos \\left( \\frac{z e C}{2 k_{B} T} \\frac{d}{2} \\right)} \\right),\n", - "\\end{equation}\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", - "\\begin{equation}\n", + "$$\n", "\\sigma(x) = \\mu \\frac{\\partial v_{y}(x)}{\\partial x}\n", - "\\end{equation}" + "$$" ] }, { @@ -428,14 +437,14 @@ "SINGLE_PRECISION = False\n", "\n", "padding = 1\n", - "WIDTH = 200\n", - "BOX_L = [WIDTH + 2 * padding, 1, 1]\n", + "WIDTH = 126\n", + "BOX_L = [(WIDTH + 2 * padding) * AGRID, 1, 1]\n", "\n", "system.cell_system.skin = 0.4\n", "system.box_l = BOX_L\n", "system.time_step = TAU\n", "\n", - "RUN_TIME = 600" + "RUN_TIME = 200" ] }, { @@ -461,7 +470,9 @@ "outputs": [], "source": [ "viscosity_kinematic = VISCOSITY_DYNAMIC / DENSITY_FLUID\n", - "lbf = espressomd.lb.LBFluidWalberla(lattice=lattice, density=DENSITY_FLUID, kinematic_viscosity=viscosity_kinematic, tau=TAU, single_precision=SINGLE_PRECISION)\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" ] }, @@ -494,8 +505,8 @@ }, "source": [ "```python\n", - "eksolver = espressomd.electrokinetics.EKFFT(lattice=lattice, permittivity=PERMITTIVITY, single_precision=SINGLE_PRECISION)\n", - "\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", "```" ] @@ -600,10 +611,10 @@ }, "source": [ "```python\n", - "for shape_obj in (wall_left, wall_right):\n", - " ekspecies.add_boundary_from_shape(shape=shape_obj, value=[0., 0., 0.], boundary_type=espressomd.electrokinetics.FluxBoundary)\n", - " ekspecies.add_boundary_from_shape(shape=shape_obj, value=0.0, boundary_type=espressomd.electrokinetics.DensityBoundary)\n", - " lbf.add_boundary_from_shape(shape=shape_obj, velocity=[0., 0., 0.])\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", "```" ] }, @@ -629,7 +640,7 @@ }, "outputs": [], "source": [ - "for i in tqdm.trange(100):\n", + "for i in tqdm.trange(80):\n", " system.integrator.run(RUN_TIME)" ] }, @@ -699,7 +710,7 @@ "fig1.suptitle(\"electroosmotic flow\")\n", "\n", "ax = fig1.add_subplot(131)\n", - "ax.plot(positions, density_eof, \"o\", mfc=\"none\", label=\"simulation\")\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", @@ -707,7 +718,7 @@ "ax.legend(loc=\"best\")\n", "\n", "ax = fig1.add_subplot(132)\n", - "ax.plot(positions, velocity_eof, \"o\", mfc=\"none\", label=\"simulation\")\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", @@ -715,7 +726,7 @@ "ax.legend(loc=\"best\")\n", "\n", "ax = fig1.add_subplot(133)\n", - "ax.plot(positions, pressure_tensor_eof, \"o\", mfc=\"none\", label=\"simulation\")\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", @@ -737,7 +748,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Comparison to pressure driven flow\n", + "### 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." ] }, @@ -761,7 +772,7 @@ }, "outputs": [], "source": [ - "for i in tqdm.trange(50):\n", + "for i in tqdm.trange(70):\n", " system.integrator.run(RUN_TIME)" ] }, @@ -780,7 +791,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The anlaytic solution for pressure driven flow between two infinite parallel plates is known as the Poiseuille flow." + "The analytic solution for pressure-driven flow between two infinite parallel plates is known as the Poiseuille flow." ] }, { @@ -813,10 +824,10 @@ "outputs": [], "source": [ "fig1 = plt.figure(figsize=(16, 4.5))\n", - "fig1.suptitle(\"pressure driven flow\")\n", + "fig1.suptitle(\"pressure-driven flow\")\n", "\n", "ax = fig1.add_subplot(131)\n", - "ax.plot(positions, density_pressure, \"o\", mfc=\"none\", label=\"simulation\")\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", @@ -824,7 +835,7 @@ "ax.legend(loc=\"best\")\n", "\n", "ax = fig1.add_subplot(132)\n", - "ax.plot(positions, velocity_pressure, \"o\", mfc=\"none\", label=\"simulation\")\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", @@ -832,7 +843,7 @@ "ax.legend(loc=\"best\")\n", "\n", "ax = fig1.add_subplot(133)\n", - "ax.plot(positions, pressure_tensor_pressure, \"o\", mfc=\"none\", label=\"simulation\")\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", @@ -867,27 +878,27 @@ "outputs": [], "source": [ "fig1 = plt.figure(figsize=(16, 4.5))\n", - "fig1.suptitle(\"electroosmotic - pressure driven flow comparison\")\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\", label=\"eof\")\n", - "ax.plot(positions, density_pressure, \"o\", mfc=\"none\", label=\"pressure\")\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\", label=\"eof\")\n", - "ax.plot(positions, velocity_pressure, \"o\", mfc=\"none\", label=\"pressure\")\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\", label=\"eof\")\n", - "ax.plot(positions, pressure_tensor_pressure, \"o\", mfc=\"none\", label=\"pressure\")\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", @@ -941,18 +952,18 @@ "metadata": {}, "outputs": [], "source": [ - "BOX_L = [80, 30, 1] \n", + "BOX_L = [80, 32, 1] \n", "AGRID = 1.0\n", "DIFFUSION_COEFFICIENT = 0.01\n", - "TAU = 0.01\n", + "TAU = 0.03\n", "EXT_FORCE_DENSITY = [0.6, 0, 0]\n", - "OBSTACLE_RADIUS = 5\n", + "OBSTACLE_RADIUS = 6\n", "\n", "DENSITY_FLUID = 0.5\n", "VISCOSITY_KINEMATIC = 2.0\n", "KT = 1.0\n", "\n", - "TOTAL_FRAMES = 200\n", + "TOTAL_FRAMES = 50\n", "\n", "system.time_step = TAU\n", "system.cell_system.skin = 0.4\n", @@ -965,15 +976,12 @@ "metadata": {}, "outputs": [], "source": [ - "lattice = espressomd.lb.LatticeWalberla(\n", - " n_ghost_layers=1, agrid=AGRID)\n", - "\n", + "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", "\n", - "\n", "eksolver = espressomd.electrokinetics.EKNone(lattice=lattice)\n", "system.ekcontainer = espressomd.electrokinetics.EKContainer(tau=TAU, solver=eksolver)" ] @@ -988,7 +996,7 @@ "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." + "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." ] }, { @@ -997,8 +1005,7 @@ "metadata": {}, "outputs": [], "source": [ - "REACTION_RATE_CONSTANT = 0.1\n", - "\n", + "REACTION_RATE_CONSTANT = 2.5\n", "EDUCT_COEFFS = [-1, -1]\n", "EDUCT_ORDERS = [1,1]\n", "PRODUCT_COEFFS = [1]\n", @@ -1089,8 +1096,8 @@ "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 vorticies 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." + "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." ] }, { @@ -1115,10 +1122,10 @@ " ))\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)" + " 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)" ] }, { @@ -1180,35 +1187,38 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": true - }, + "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", - "progress_bar = tqdm.tqdm(total=TOTAL_FRAMES)\n", - "\n", - "cmap = mpl.cm.get_cmap(\"viridis\").copy()\n", + "cmap = get_colormap(\"viridis\").copy()\n", "cmap.set_bad(color=\"gray\")\n", - "cmap_quiver = mpl.cm.get_cmap(\"binary\").copy()\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=(20, 10))\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)\n", + "ax4.imshow(bg_colors, cmap=cmap_quiver, **imshow_kwargs)\n", "\n", "for ax, title in zip(\n", " [ax1, ax2, ax3, ax4],\n", @@ -1218,16 +1228,19 @@ " ax.set_xlim((0, box_width))\n", " ax.set_ylim((0, box_height))\n", "\n", - "# create meshgrid for quiver plot\n", - "xs = np.arange(box_width)\n", - "ys = np.arange(box_height)\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\n", - "quiver = ax4.quiver(X, Y, flowfield[..., 1], flowfield[..., 0], scale=100)\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(100)\n", + " system.integrator.run(50)\n", " \n", " flowfield = np.copy(lbf[:, :, 0].velocity)\n", " e1 = np.copy(educt_species[0][:, :, 0].density)\n", @@ -1240,14 +1253,14 @@ " p[boundary_mask] = np.NaN\n", " flowfield[boundary_mask] = np.NaN\n", "\n", - "\n", - " ax1.imshow(e1, cmap=cmap)\n", - " ax2.imshow(e2, cmap=cmap)\n", - " ax3.imshow(p, cmap=cmap)\n", - " quiver.set_UVC(flowfield[..., 1], flowfield[..., 0])\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", "\n", "animation.FuncAnimation(fig, draw_frame, frames=range(TOTAL_FRAMES), interval=300)" ] diff --git a/src/python/espressomd/electrokinetics.py b/src/python/espressomd/electrokinetics.py index 2ba4138a69..336cea4780 100644 --- a/src/python/espressomd/electrokinetics.py +++ b/src/python/espressomd/electrokinetics.py @@ -596,10 +596,11 @@ class EKReactant(ScriptInterfaceHelper): Parameters ---------- - ekspecies : :obj:`espressomd.electrokinetics.EKSpecies ` + ekspecies : :obj:`~espressomd.electrokinetics.EKSpecies` EK species to react stoech_coeff: :obj:`float` - Stoechiometric coefficient of this species in the reaction. Products have positive coeffiecients whereas Educts have negative ones. + 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. @@ -620,7 +621,7 @@ class EKBulkReaction(ScriptInterfaceHelper): 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: array_like of :obj:`~espressomd.electrokinetics.EKReactant` Reactants that participate this reaction. """ @@ -641,7 +642,7 @@ class EKIndexedReaction(ScriptInterfaceHelper): 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: array_like of :obj:`~espressomd.electrokinetics.EKReactant` Reactants that participate this reaction. """ @@ -726,7 +727,7 @@ def add(self, reaction): Parameters ---------- - reaction : :obj:`espressomd.electrokinetics.EKBulkReaction ` or :obj:`espressomd.electrokinetics.EKIndexedReaction ` + reaction : :obj:`~espressomd.electrokinetics.EKBulkReaction` or :obj:`~espressomd.electrokinetics.EKIndexedReaction` Reaction to be added. """ @@ -740,7 +741,7 @@ def remove(self, reaction): Parameters ---------- - reaction : :obj:`espressomd.electrokinetics.EKBulkReaction ` or :obj:`espressomd.electrokinetics.EKIndexedReaction ` + reaction : :obj:`~espressomd.electrokinetics.EKBulkReaction` or :obj:`~espressomd.electrokinetics.EKIndexedReaction` Reaction to be removed. """ @@ -750,13 +751,13 @@ def remove(self, reaction): @script_interface_register class EKContainer(ScriptObjectList): """ - Container object holding the :obj:`espressomd.electrokinetics.EKSpecies `. + 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 : :obj:`~espressomd.electrokinetics.EKNone` or :obj:`~espressomd.electrokinetics.EKFFT` Solver defining the treatment of the electrostatic Poisson-equation. Methods @@ -782,19 +783,18 @@ def reactions(self): """ Returns ------- - :obj:`espressomd.electrokinetics.EKReactions ` - Reactions-container of the current reactions. + 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. + Add an :obj:`~espressomd.electrokinetics.EKSpecies` to the container. Parameters ---------- - ekspecies : :obj:`espressomd.electrokinetics.EKSpecies ` + ekspecies : :obj:`~espressomd.electrokinetics.EKSpecies` Species to be added. """ @@ -802,11 +802,11 @@ def add(self, ekspecies): def remove(self, ekspecies): """ - Remove an :obj:`espressomd.electrokinetics.EKSpecies ` from the container. + Remove an :obj:`~espressomd.electrokinetics.EKSpecies` from the container. Parameters ---------- - ekspecies : :obj:`espressomd.electrokinetics.EKSpecies ` + ekspecies : :obj:`~espressomd.electrokinetics.EKSpecies` Species to be removed. """ diff --git a/testsuite/scripts/tutorials/test_electrokinetics.py b/testsuite/scripts/tutorials/test_electrokinetics.py index 4cd7264779..e138fe209a 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,25 +15,109 @@ # # 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 as iw import numpy as np +import scipy.optimize tutorial, skipIfMissingFeatures = iw.configure_and_import( - "@TUTORIALS_DIR@/electrokinetics/electrokinetics.py", RUN_TIME=10, TOTAL_FRAMES=10) + "@TUTORIALS_DIR@/electrokinetics/electrokinetics.py", TOTAL_FRAMES=0) @skipIfMissingFeatures class Tutorial(ut.TestCase): system = tutorial.system - def test_simulation_fields_finite(self): + 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__":