diff --git a/doc/notebooks/07_tutorial_shanchen_twophase.ipynb b/doc/notebooks/07_tutorial_shanchen_twophase.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..b6dbc8cae666eb2d98eaaa5d80c7eb16144fdeca
--- /dev/null
+++ b/doc/notebooks/07_tutorial_shanchen_twophase.ipynb
@@ -0,0 +1,334 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Shan-Chen Two-Phase Single-Component Lattice Boltzmann"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from lbmpy.session import *\n",
+    "from lbmpy.updatekernels import create_stream_pull_with_output_kernel\n",
+    "from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This is based on section 9.3.2 of Krüger et al.'s \"The Lattice Boltzmann Method\", Springer 2017 (http://www.lbmbook.com).\n",
+    "Sample code is available at [https://github.com/lbm-principles-practice/code/](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "N = 64\n",
+    "omega_a = 1.\n",
+    "g_aa = -4.7\n",
+    "rho0 = 1."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Data structures"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dh = ps.create_data_handling((N, N), periodicity=True, default_target='cpu')\n",
+    "\n",
+    "method_a = create_lb_method(relaxation_rate=omega_a, compressible=True)\n",
+    "\n",
+    "src_a = dh.add_array('src_a', values_per_cell=len(method_a.stencil))\n",
+    "dst_a = dh.add_array_like('dst_a', 'src_a')\n",
+    "\n",
+    "ρ_a = dh.add_array('rho_a')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Force & combined velocity"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The force on the fluid is\n",
+    "$\\vec{F}_A(\\vec{x})=-\\psi(\\rho_A(\\vec{x}))g_{AA}\\sum\\limits_{i=1}^{19}w_i\\psi(\\rho_A(\\vec{x}+\\vec{c}_i))\\vec{c}_i$\n",
+    "with \n",
+    "$\\psi(\\rho)=\\rho_0\\left[1-\\exp(-\\rho/\\rho_0)\\right]$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def psi(dens):\n",
+    "    return rho0 * (1. - sp.exp(-dens / rho0));"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "zero_vec = sp.Matrix([0] * dh.dim) \n",
+    "\n",
+    "stencil, weights = method_a.stencil, method_a.weights\n",
+    "force_a = sum((psi(ρ_a[d]) * w_d * sp.Matrix(d)\n",
+    "                for d, w_d in zip(stencil, weights)), \n",
+    "               zero_vec) * psi(ρ_a.center) * -1 * g_aa"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Kernels"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "stream_a = create_stream_pull_with_output_kernel(method_a, src_a, dst_a, {'density': ρ_a})\n",
+    "\n",
+    "# TODO use method above\n",
+    "collision_a = create_lb_update_rule(relaxation_rate=omega_a, \n",
+    "                                    compressible=True,\n",
+    "                                    force_model='guo', \n",
+    "                                    force=force_a,\n",
+    "                                    kernel_type='collide_only',\n",
+    "                                    optimization={'symbolic_field': src_a})\n",
+    "\n",
+    "opts = {'cpu_openmp': False}\n",
+    "stream_a_kernel = ps.create_kernel(stream_a, **opts).compile()\n",
+    "collision_a_kernel = ps.create_kernel(collision_a, **opts).compile()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Initialization"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "init_a = macroscopic_values_setter(method_a, velocity=(0, 0), \n",
+    "                                   pdfs=src_a.center_vector, density=ρ_a.center)\n",
+    "init_a_kernel = ps.create_kernel(init_a, ghost_layers=0).compile()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def init():\n",
+    "    for x in range(N):\n",
+    "        for y in range(N):\n",
+    "            if (x-N/2)**2 + (y-N/2)**2 <= 15**2:\n",
+    "                dh.fill(ρ_a.name, 2.1, slice_obj=[x,y])\n",
+    "            else:\n",
+    "                dh.fill(ρ_a.name, 0.15, slice_obj=[x,y])\n",
+    "\n",
+    "    dh.run_kernel(init_a_kernel)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Timeloop"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sync_pdfs = dh.synchronization_function([src_a.name])\n",
+    "sync_ρs = dh.synchronization_function([ρ_a.name])\n",
+    "\n",
+    "def time_loop(steps):\n",
+    "    dh.all_to_gpu()\n",
+    "    for i in range(steps):\n",
+    "        sync_ρs()\n",
+    "        dh.run_kernel(collision_a_kernel)\n",
+    "        \n",
+    "        sync_pdfs()\n",
+    "        dh.run_kernel(stream_a_kernel)\n",
+    "        \n",
+    "        dh.swap(src_a.name, dst_a.name)\n",
+    "    dh.all_to_cpu()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_ρs():\n",
+    "    plt.title(\"$\\\\rho_A$\")\n",
+    "    plt.scalar_field(dh.gather_array(ρ_a.name), vmin=0, vmax=2.5)\n",
+    "    plt.colorbar()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Run the simulation\n",
+    "### Initial state"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAF0CAYAAAAKF1nQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAeTklEQVR4nO3df6zld1kn8PdDKaKAC9ihzPYHxWRWLYZfOymwbKRaWUtt2v4Bu2Wj22XZNLqwC0YjBaPGZDfBaIgaBHYCSFkRlgDSiSlgrQiaSGVayo92YFtRYWTstAX5IQjpzLN/nANehztz7/nee+45535fr+Sbc74/zvc8uZ/b9j59Pp/nW90dAACAMXjQogMAAADYKRIgAABgNCRAAADAaEiAAACA0ZAAAQAAoyEBAgAARkMCBAAAjIYECAAAGA0JEMBIVNUjqupAVX2hqo5V1c8sOiYA2GkSIIDxeHeSv0zy2CRXJ/n1qnrsYkMCgJ0lAQIYgaq6PEm6+1e7++vd/cdJ/jbJv5qef0JVHa+qcxcZJwDMmwQIYByuSHLDN3eq6kFJ/kWSe6aHXpbk/yT5gZ0PDQB2jgQIYByeluT+Nfs/kuS+7v5UVT0xydEk74sECIBdTgIEsMtV1ZlJ9iV5blU9tKqekOQ1mVR9kuRnkvxqkjsjAQJgl3vwogMAYO5+IMlfJ/lEJlPejiX5n939jqp6cpJnJvm9JGdMNwDYtSRAALvfE5Mc7u5fTPKLJ527LsnTuvsLSVJVf7HTwQHATjIFDmD3e1KSwycfrKp/neRr30x+pv6xqr5nxyIDgB0mAQLY/Z6Y5JMnH+zuW7v7BScd+6Huvv/kawFgp1XVeVX1/qo6XFV3VNVL1rnm4qr6YlXdPt1+acP7dvd8IgYAABioqvYm2dvdt1XVI5LcmuSq7r5zzTUXJ/m57r58s/dVAQIAAJZOdx/t7tum77+cyXTuc7Z6XwkQAACw1KrqgiRPSXLLOqefUVUfrar3TB/1cFo72gXurLPO6gsuuGAnvxKWzl0f/cyiQwBYOvuedP6iQ4CFuvXWW+/r7j2LjmMWP/bDD+v7P3988Odv/djX70jyj2sOHejuAydfV1UPT/LOJC/t7i+ddPq2JI/r7q9U1WVJ3p3Js+9OaUcToAsuuCCHDh3aya+EpXPp2f9t0SEALJ33HnrNokOAhaqqv1l0DLO67/PHc8v7zh38+TP3/uU/dvf+010zfZj3O5O8pbvfdfL5tQlRd99YVa+pqrO6+75T3dMUOAAAYOlUVSV5QybPsnvVKa557PS6VNVFmeQ3p+1m6kGoAADAAJ3jfWKeX/DMJD+Z5ONVdfv02CuSnJ8k3f26JM9N8tNV9UCSryW5ujdocy0BAgAAZtZJTmR+j9Tp7j9LUhtc8+okr57lvhIgAABgkBOZawVoLqwBAgAARkMFCAAAmFmnc/z0y22WkgQIAAAYZJ5rgOZFAgQAAMyskxyXAAEAAGOxihUgTRAAAIDRUAECAABm1okmCAAAwHis3lOANjkFrqoeWVXvqKpPVtXhqnpGVT26qm6qqrumr4+ad7AAAMBy6HSOb2FblM2uAfrNJO/t7u9P8qQkh5Ncl+Tm7t6X5ObpPgAAMAadHN/CtigbJkBV9d1JfijJG5Kku7/R3X+f5Mok108vuz7JVfMKEgAAYDtspgL0vUnuTfI7VfWRqnp9VT0sydndfTRJpq+PWe/DVXVtVR2qqkP33nvvtgUOAAAsTmeyBmjotiibSYAenOSpSV7b3U9J8g+ZYbpbdx/o7v3dvX/Pnj0DwwQAAJZL5fgWtkXZTAJ0JMmR7r5luv+OTBKie6pqb5JMX4/NJ0QAAGDZdJITPXxblA0ToO7+uySfrarvmx66JMmdSQ4muWZ67JokN8wlQgAAgG2y2ecA/fckb6mqhyT5dJIXZJI8vb2qXpjkM0meN58QAQCAZbTIqWxDbSoB6u7bk+xf59Ql2xsOAACwCjq7OAECAAA42YmWAAEAACOwqhWgzXSBAwAA2BVUgAAAgJl1KsdXsJ4iAQIAAAaxBggAABiFVV0DJAECAAAGqBzv1ZsCt3oRAwAADKQCBAAAzKyTnFjBeooECAAAGMQaIAAAYBS6rQECAABYaipAAADAICdMgQMAAMZg8hyg1ZtQJgECAAAGWM01QBIgAABgZqvaBnv1IgYAABhIBQgAABjkeGuCAAAAjECnNEEAAADG44QmCAAAwBisahvs1YsYAABgIBUgAABgZp3SBAEAABiPVXwOkAQIAACYWXdyfAWbIKxexAAAAAOpAAEAAANUTsQaIAAAYAQ6qzkFTgIEAAAMsorPAZIAAQAAM+tUTqxgG+zVS9kAAAAGUgECAAAGMQUOAAAYhU5yQhMEAABgHCrHtcEGAADGYFUrQKsXMQAAwEAqQAAAwCCmwAEAAKPQXabAAQAA43G8HzR420hVnVdV76+qw1V1R1W9ZJ1rqqp+q6rurqqPVdVTN7qvChAAALCMHkjys919W1U9IsmtVXVTd9+55prnJNk33Z6W5LXT11NSAQIAAGbWSU6kBm8b3r/7aHffNn3/5SSHk5xz0mVXJnlzT3woySOrau/p7qsCBAAADFCbmsq2Ld9UdUGSpyS55aRT5yT57Jr9I9NjR091r00lQFX110m+nOR4kge6e39VPTrJ/01yQZK/TvLvu/sLm7kfAACw2ibPAdpSF7izqurQmv0D3X3g5Iuq6uFJ3pnkpd39pZNPnyK0U5qlAvTD3X3fmv3rktzc3a+squum+y+b4X4AAMAKO761FTX3dff+011QVWdmkvy8pbvftc4lR5Kct2b/3CSfO909txLxlUmun76/PslVW7gXAADAt1RVJXlDksPd/apTXHYwyX+adoN7epIvdvcpp78lm68AdZI/rKpO8r+npamzv3nz7j5aVY85ReDXJrk2Sc4///xNfh0AALDMOrXVKXAbeWaSn0zy8aq6fXrsFUnOT5Lufl2SG5NcluTuJF9N8oKNbrrZBOiZ3f25aZJzU1V9crNRT5OlA0myf//+087HAwAAVseJOTaV7u4/y/prfNZe00leNMt9N5UAdffnpq/Hqur3k1yU5J6q2jut/uxNcmyWLwYAAFZXd3J8vhWgudgwZauqh00fPJSqeliSf5fkE5nMt7tmetk1SW6YV5AAAMDyOdE1eFuUzVSAzk7y+5M1SHlwkt/r7vdW1YeTvL2qXpjkM0meN78wAQAAtm7DBKi7P53kSescvz/JJfMICgAAWG6TJgg78yDU7TTLc4AAAAC+5fjpexQsJQkQAAAws04WupZnqNWrWQEAAAykAgQAAAxgDRAAADAiJ6wBAgAAxmBVH4QqAQIAAAYxBQ6AlXHFB+5c9/jBZ124w5EAwM6RAAEAADObPAjVFDgAAGAkNEEAAABGwYNQAQAAlpwKEMASOlWDgt3y3RotAOwOusABAADj0JogAAAAI9HRBAEAABiRVawArd6kPQAAgIFUgAB20CKbGyyTzf4cNEsAWF6r2gZbAgQAAAwiAQIAAEahowscAAAwIqvYBU4TBAAAYDRUgAC2geYG8zHLz1XDBIAd1tYAAQAAI6ELHAAAMCqrmABZAwQAAIyGChAAADAzbbABRkCzg+W13thojAAwXy0BAgAAxmIVnwMkAQIAAGbWK9oGWxMEAABgNFSAAACAQawBAthFNDxYfRojAMyTLnAAAMCIqAABAACj0NEEAQAAYKmpAAEAALPrSSvsVSMBAoiGB2OiMQLA9vEgVAAAYBQ6q9kEwRogAABgNFSAAACAATwHCAAAGJFVbIKw6SlwVXVGVX2kqv5guv/oqrqpqu6avj5qfmECAADLprsGb4syyxqglyQ5vGb/uiQ3d/e+JDdP9wEAgBHo3sUJUFWdm+THk7x+zeErk1w/fX99kqu2NzQAAIDttdk1QL+R5OeTPGLNsbO7+2iSdPfRqnrMeh+sqmuTXJsk559//hZCBQAAlskqNkHYsAJUVZcnOdbdtw75gu4+0N37u3v/nj17htwCAABYQpNpcMO2RdlMBeiZSa6oqsuSPDTJd1fV7ya5p6r2Tqs/e5Mcm2egAADAclnFB6FumAB198uTvDxJquriJD/X3T9RVb+W5Jokr5y+3jDHOAG2zRUfuHPRIbBk1vudOPisCxcQCcDq6Cy2mcFQs3SBO9krkzy7qu5K8uzpPgAAwNKa6UGo3f0nSf5k+v7+JJdsf0gAAMAqWMHnoM6WAAEAACRJepeuAQIAAFjXCpaAtrIGCAAAYC6q6o1VdayqPnGK8xdX1Rer6vbp9kubua8KEAAAMMicp8C9Kcmrk7z5NNf8aXdfPstNJUAAAMAg83ygaXd/sKou2O77mgIHAADMrDOpAA3dkpxVVYfWbNcOCOMZVfXRqnpPVT1hMx9QAQIAAGbXSbY2Be6+7t6/hc/fluRx3f2VqrosybuT7NvoQxIgYFe74gN3LjoEVtR6vzsHn3XhAiIBYD3d/aU172+sqtdU1Vndfd/pPicBAgAABpnnGqCNVNVjk9zT3V1VF2WyvOf+jT4nAQIAAIaZYwJUVW9NcnEma4WOJPnlJGcmSXe/Lslzk/x0VT2Q5GtJru7eOCWTAAEAAAPUXNtgd/fzNzj/6kzaZM9EAgQAAAyzwClwQ2mDDQAAjIYKEAAAMLvOXKfAzYsECAAAGGYFp8BJgAAAgIFWrwJkDRAAADAaKkAAAMAwpsABAACjIQECAABGoZPoAgcAAIxFr2AFSBMEAABgNFSAAACAYVawAiQBAgAAhrEGCAAAGItSAQIAAEahs5JT4DRBAAAARkMFCAAAGKCsAQIAAEZkBafASYAAAIBhVjABsgYIAAAYDRUgAABgmBWsAEmAgF3jig/cuegQ2OXW+x07+KwLFxAJwBLoaIIAAACMhwehAgAA47GCCZAmCAAAwGhIgAAAgNEwBQ4AABjEGiCABVqvG5fOcGwnHd8ATqILHAAAMAodTRAAAACWmQoQAAAwzG6sAFXVQ6vqL6rqo1V1R1X9yvT4o6vqpqq6a/r6qPmHCwAALIvq4duibGYK3NeT/Eh3PynJk5NcWlVPT3Jdkpu7e1+Sm6f7AADAWPQWtgXZMAHqia9Md8+cbp3kyiTXT49fn+SquUQIAACwTTbVBKGqzqiq25McS3JTd9+S5OzuPpok09fHnOKz11bVoao6dO+9925X3AAAwKLtxgpQknT38e5+cpJzk1xUVT+42S/o7gPdvb+79+/Zs2donAAAwBLZyvqfZV8D9C3d/fdJ/iTJpUnuqaq9STJ9Pbbt0QEAAMura/i2IJvpArenqh45ff+dSX40ySeTHExyzfSya5LcMK8gAQCAJbSCU+A28xygvUmur6ozMkmY3t7df1BVf57k7VX1wiSfSfK8OcYJAACwZRsmQN39sSRPWef4/UkumUdQAADA8lvkWp6hNlMBAgAA+HYSIAAAYBQW3M1tqJm6wAEAAKwyFSAAAGCYFawASYAAAIBhJEAAAMBYWAMEAACwxCRAAADAaJgCBwAADLOCU+AkQAAAwOxW9DlAEiAAAGAYCRAAADAaEiCA5XLwWRd+27ErPnDnAiJh1az3uwPA6pMAAQAAM6tYAwQAAIyJBAgAABiFFe0C50GoAADA0qmqN1bVsar6xCnOV1X9VlXdXVUfq6qnbua+EiAAAGCY3sK2sTclufQ055+TZN90uzbJazdzUwkQAAAwzBwToO7+YJLPn+aSK5O8uSc+lOSRVbV3o/taAwQAAAyyxTVAZ1XVoTX7B7r7wAyfPyfJZ9fsH5keO3q6D0mAAACAYbaWAN3X3fu38Pla59iGEZkCBwAArKIjSc5bs39uks9t9CEVIGB0Dj7rwm87dsUH7lxAJCyL9X4nANjA5psZzMvBJC+uqrcleVqSL3b3aae/JRIgAABgoHk+B6iq3prk4kzWCh1J8stJzkyS7n5dkhuTXJbk7iRfTfKCzdxXAgQAAAwzxwSou5+/wflO8qJZ7ysBAgAABplnBWheNEEAAABGQwUIAAAYZgUrQBIgAABgdovvAjeIBAgAAJhZZf0nkS47a4AAAIDRUAECAACGMQUOYDUdfNaF33bsig/cuYBImLf1xhqAYVaxDbYECAAAGEYCBAAAjMYKJkCaIAAAAKOhAgQAAMyurQECAADGRAIEsHvoDLf6dHwDmC8VIAAAYDxWMAHSBAEAABgNFSAAAGCQVZwCt2EFqKrOq6r3V9Xhqrqjql4yPf7oqrqpqu6avj5q/uECAABLobe4LchmKkAPJPnZ7r6tqh6R5NaquinJf05yc3e/sqquS3JdkpfNL1SAxTvVonrNERZPwwOABdiNFaDuPtrdt03ffznJ4STnJLkyyfXTy65PctW8ggQAANgOM60BqqoLkjwlyS1Jzu7uo8kkSaqqx5ziM9cmuTZJzj///K3ECgAALInKLl0D9E1V9fAk70zy0u7+0mY/190Hunt/d+/fs2fPkBgBAIBltEvXAKWqzswk+XlLd79revieqto7rf7sTXJsXkECAADLp3r1SkAbJkBVVUnekORwd79qzamDSa5J8srp6w1ziRBgBWx2Ab5mCbPR2ABgiS24kjPUZipAz0zyk0k+XlW3T4+9IpPE5+1V9cIkn0nyvPmECAAAsD02TIC6+88yWeO0nku2NxwAAGBVrGIThJm6wAEAAHyLBAgAABgLFSAATkuzhAnNDQB2iRVMgDb9HCAAAIBVpwIEAADMrk2BAwAAxkQCBAAAjEFFBQiAbbITTQJO1WhBgwIAdjMJEAAAMEyvXglIAgQAAAxiChwAADAOHU0QAACA8agTi45gdhIggJHS7ACAMZIAAQAAw5gCBwAAjIUmCAAAwDh0tMEGAADGYxUrQA9adAAAAAA7RQUIAAAYZgUrQBIgAABgZpXVnAInAQIAAGbXvZJNEKwBAgAARkMFCAAAGMQUOAAAYDwkQAAAwFioAAEAAOPQSU6sXgakCQIAADAaKkAAAMAwq1cAkgABAADDWAMEAACMhwehAgAAY1E9fNvU/asurapPVdXdVXXdOucvrqovVtXt0+2XNrqnChAAALB0quqMJL+d5NlJjiT5cFUd7O47T7r0T7v78s3eVwUIAACYXW9x29hFSe7u7k939zeSvC3JlVsNWwIEAADMrJJU9+BtE85J8tk1+0emx072jKr6aFW9p6qesNFNTYEDAACGObGlT59VVYfW7B/o7gNr9mudz5ycOd2W5HHd/ZWquizJu5PsO92XSoAAAIBFuK+795/m/JEk563ZPzfJ59Ze0N1fWvP+xqp6TVWd1d33neqmpsABAACDzHkK3IeT7Kuqx1fVQ5JcneTgP/v+qsdWVU3fX5RJfnP/6W6qAgQAAMxu880Mht2++4GqenGS9yU5I8kbu/uOqvqp6fnXJXlukp+uqgeSfC3J1d2nz64kQAAAwAA99wehdveNSW486djr1rx/dZJXz3JPCRAAADDIZh9oukysAQIAAEZDBQgAABhmzlPg5mHDClBVvbGqjlXVJ9Yce3RV3VRVd01fHzXfMAEAgKXSSZ0Yvi3KZqbAvSnJpScduy7Jzd29L8nN030AAGBMuodvC7JhAtTdH0zy+ZMOX5nk+un765Nctc1xAQAAbLuhTRDO7u6jSTJ9fcypLqyqa6vqUFUduvfeewd+HQAAsHR6C9uCzL0LXHcf6O793b1/z5498/46AABgh1T34G1RhiZA91TV3iSZvh7bvpAAAICVsBvXAJ3CwSTXTN9fk+SG7QkHAABYCZ3kxBa2BdlMG+y3JvnzJN9XVUeq6oVJXpnk2VV1V5JnT/cBAACW2oYPQu3u55/i1CXbHAsAALAiKotdyzPUhgkQAADAuiRAAADAaEiAAACAUfhmE4QVM/fnAAEAACwLFSAAAGAQTRAAAIDxkAABAADj0CuZAFkDBAAAjIYKEAAAMLvOSlaAJEAAAMAwK9gGWwIEAAAMogscAAAwHiuYAGmCAAAAjIYKEAAAMLtOcmL1KkASIAAAYIDVfA6QBAgAABhGAgQAAIzGCiZAmiAAAACjoQIEAADMThMEAABgPDrpE4sOYmYSIAAAYBhrgAAAAJaXChAAADA7a4AAAIBRWcEpcBIgAABgGAkQAAAwDr2SCZAmCAAAwGioAAEAALPrJCc8BwgAABiLFZwCJwECAACGkQABAADj0Cv5HCBNEAAAgNFQAQIAAGbXSbcmCAAAwFis4BQ4CRAAADDMCjZBsAYIAAAYDRUgAABgdt0ehAoAAIzICk6BkwABAACDtAoQAAAwDr2SFSBNEAAAgNFQAQIAAGbXWcnnAG2pAlRVl1bVp6rq7qq6bruCAgAAVkCfGL4tyOAKUFWdkeS3kzw7yZEkH66qg91953YFBwAALKdO0iOrAF2U5O7u/nR3fyPJ25JcuT1hAQAAS6177hWgjWac1cRvTc9/rKqeutE9t5IAnZPks2v2j0yPnRzUtVV1qKoO3XvvvVv4OgAAYCzWzDh7TpILkzy/qi486bLnJNk33a5N8tqN7ruVBKjWOfZtNbDuPtDd+7t7/549e7bwdQAAwDLpEz1424TNzDi7Msmbe+JDSR5ZVXtPd9OtJEBHkpy3Zv/cJJ/bwv0AAIBVMt8pcJuZcbapWWlrbaUN9oeT7Kuqxyf52yRXJ/mPp/vArbfeel9V/c0WvpPhzkpy36KD4NsYl+VjTJaTcVlO2zYuVRvOWmHz/POynDYal8ftVCDb5cv5wvv+qN9x1hZu8dCqOrRm/0B3H1izv5kZZ5ualbbW4ASoux+oqhcneV+SM5K8sbvv2OAz5sAtSFUd6u79i46Df864LB9jspyMy3IyLsvJuCyn3Tgu3X3pnL9iMzPOZp6VtqUHoXb3jUlu3Mo9AAAA1rGZGWcHk7y4qt6W5GlJvtjdR0930y0lQAAAAPNwqhlnVfVT0/Ovy6QYc1mSu5N8NckLNrqvBGg8Dmx8CQtgXJaPMVlOxmU5GZflZFyWk3EZYL0ZZ9PE55vvO8mLZrlnTT4DAACw+22lDTYAAMBKkQDtYlX1a1X1yar6WFX9flU9cs25l1fV3VX1qar6sUXGOUZVden0Z393VV236HjGqqrOq6r3V9Xhqrqjql4yPf7oqrqpqu6avj5q0bGOTVWdUVUfqao/mO4bkyVQVY+sqndM/9tyuKqeYWwWq6p+Zvrvr09U1Vur6qHGZDGq6o1VdayqPrHm2CnHwt9iiyMB2t1uSvKD3f3EJP8vycuTpKouzKSLxhOSXJrkNVV1xsKiHJnpz/q3kzwnyYVJnj8dE3beA0l+trt/IMnTk7xoOhbXJbm5u/cluXm6z856SZLDa/aNyXL4zSTv7e7vT/KkTMbI2CxIVZ2T5H8k2d/dP5jJIvGrY0wW5U2Z/F211rpj4W+xxZIA7WLd/Yfd/cB090OZ9EVPkiuTvK27v97df5VJ14yLFhHjSF2U5O7u/nR3fyPJ2zIZE3ZYdx/t7tum77+cyR9z52QyHtdPL7s+yVWLiXCcqurcJD+e5PVrDhuTBauq707yQ0nekCTd/Y3u/vsYm0V7cJLvrKoHJ/muTJ5/YkwWoLs/mOTzJx0+1Vj4W2yBJEDj8V+SvGf6/pwkn11z7sj0GDvDz38JVdUFSZ6S5JYkZ3/zGQLT18csLrJR+o0kP5/kxJpjxmTxvjfJvUl+Zzo98fVV9bAYm4Xp7r9N8utJPpPkaCbPP/nDGJNlcqqx8LfAAkmAVlxV/dF03u/J25VrrvmFTKb6vOWbh9a5lXaAO8fPf8lU1cOTvDPJS7v7S4uOZ8yq6vIkx7r71kXHwrd5cJKnJnltdz8lyT/E1KqFmq4nuTLJ45P8yyQPq6qfWGxUbJK/BRbIc4BWXHf/6OnOV9U1SS5Pckn/U8/zI0nOW3PZuZmUzNkZfv5LpKrOzCT5eUt3v2t6+J6q2tvdR6tqb5Jji4twdJ6Z5IqquizJQ5N8d1X9bozJMjiS5Eh33zLdf0cmCZCxWZwfTfJX3X1vklTVu5L8mxiTZXKqsfC3wAKpAO1iVXVpkpcluaK7v7rm1MEkV1fVd1TV45PsS/IXi4hxpD6cZF9VPb6qHpLJIsiDC45plKqqMlnPcLi7X7Xm1MEk10zfX5Pkhp2Obay6++XdfW53X5DJPxt/3N0/EWOycN39d0k+W1XfNz10SZI7Y2wW6TNJnl5V3zX999klmaxlNCbL41Rj4W+xBfIg1F2squ5O8h1J7p8e+lB3/9T03C9ksi7ogUym/bxn/bswD9P/u/0bmXTseWN3/68FhzRKVfVvk/xpko/nn9abvCKTdUBvT3J+Jn9gPK+7T17YypxV1cVJfq67L6+q74kxWbiqenImzSkekuTTSV6Qyf9MNTYLUlW/kuQ/ZPLf848k+a9JHh5jsuOq6q1JLk5yVpJ7kvxyknfnFGPhb7HFkQABAACjYQocAAAwGhIgAABgNCRAAADAaEiAAACA0ZAAAQAAoyEBAgAARkMCBAAAjIYECAAAGI3/D3miqZb0Upr2AAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 1152x432 with 2 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "init()\n",
+    "plot_ρs()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Check the first time step against reference data\n",
+    "\n",
+    "The reference data was obtained with the [sample code](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp) after making the following changes:\n",
+    "```c++\n",
+    "const int nsteps = 1000;\n",
+    "const int noutput = 1;\n",
+    "```\n",
+    "\n",
+    "Remove the next cell if you changed the parameters at the beginning of this notebook."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "init()\n",
+    "time_loop(1)\n",
+    "ref_a = np.array([0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.136756, 0.220324, 1.2382, 2.26247, 2.26183, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.26183, 2.26247, 1.2382, 0.220324, 0.136756, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15])\n",
+    "assert np.allclose(dh.gather_array(ρ_a.name)[N//2], ref_a)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Run the simulation until converged"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAF0CAYAAAAKF1nQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df6ycV33n8c9n5v7yT4hjx7hJTKjkbQsIGtYKZbOq0qZsQxrhaAW7YdVulmUVtYIuVEXFULVVq10pVSvUdlPIWpBithQW8StRFUrTlIoilRQ7hEBi2Li0Tdy4/pUQ27n2vZ6Z7/4xE3rxOeP7zDN37sxzn/dLGt07Z555nnPnOXfuPXPO+TyOCAEAAABAHTTGXQEAAAAAWC10gAAAAADUBh0gAAAAALVBBwgAAABAbdABAgAAAFAbdIAAAAAA1AYdIAAAAAC1QQcIAAAAQG3QAQKAmrC9yfY+28/aPm77l8ZdJwAAVhsdIACoj89J+jtJL5F0m6Tftf2S8VYJAIDVRQcIAGrA9i2SFBG/HRELEfGXkv5J0r/qPf4K223bV42zngAAjBodIACohzdKuveFO7Ybkl4k6Viv6D2S/o+kH1n9qgEAsHroAAFAPbxW0qkl939S0smI+LbtV0k6KukLogMEAFjj6AABwBpne1rSLklvsj1n+xWSPqDuqI8k/ZKk35b0uOgAAQDWuKlxVwAAMHI/IukfJH1T3SlvxyX9j4j4lO0flXS9pD+R1OzdAABYs+gAAcDa9ypJhyLi1yT92kWP7ZX02oh4VpJs/+1qVw4AgNXEFDgAWPteLenQxYW2/7Wkcy90fnrO27581WoGAMAqowMEAGvfqyR96+LCiDgYEW+9qOzHI+LUxdsCALDabF9t+4u2D9l+zPY7M9vcYPs524/0br++7H4jYjQ1BgAAAICSbO+QtCMiHra9SdJBSbdGxONLtrlB0rsj4pai+2UECAAAAMDEiYijEfFw7/sz6k7nvnLY/dIBAgAAADDRbF8j6VpJD2Uefp3tr9v+fO9SD5e0qilwW7dujWuuuWY1DwlMnCcefXLcVcDE8bgrMCGYkl1nu161c9xVAMbq4MGDJyNi27jrMYif/okNceqZdunnH3x04TFJ55cU7YuIfRdvZ3ujpE9LeldEnL7o4YclvTQiztq+WdLn1L32XV+r2gG65pprdODAgdU8JDBx3nDlL467Cpg0DQbjJUmdzrhrgDH6/IH/Ne4qAGNl+x/HXYdBnXymrYe+cFXp50/v+LvzEbH7Utv0Lub9aUkfi4jPXPz40g5RRNxv+wO2t0bEyX775K8uAAAAgIlj25I+rO617N7fZ5uX9LaT7evU7d9cMs2UC6ECAAAAKCHUjpGO3l8v6eckfcP2I72y90naKUkRcbekN0n6BdstSeck3RbLxFzTAQKAcctN/Vrr0+KY7gYAlReSOiNcvxkRX9YyC2Uj4i5Jdw2yXzpAAAAAAErpqHofaK3xjxgBAAAA4F8wAgQAAABgYKFQ+9LLbSYSHSAAAAAApYxyDdCo0AECAAAAMLCQ1KYDBAAAAKAuqjgCRAgCAAAAgNpgBAgAAADAwEIiBAEAAABAfVTvKkAFp8DZfrHtT9n+lu1Dtl9ne4vtB2w/0ft62agrCwAAAGAyhELtIW7jUnQN0O9L+rOI+GFJr5Z0SNJeSQ9GxC5JD/buAwAAAKiDkNpD3MZl2Q6Q7c2SflzShyUpIhYj4ruS9kja39tsv6RbR1VJAAAAAFgJRUaAflDSCUl/ZPtrtj9ke4Ok7RFxVJJ6X6/IPdn2HbYP2D5w4sSJFas4AAAAgPEJddcAlb2NS5EO0JSk10j6YERcK+l5DTDdLSL2RcTuiNi9bdu2ktUEAAAAMFms9hC3cSnSAToi6UhEPNS7/yl1O0THbO+QpN7X46OpIgAAAIBJE5I6Uf42Lst2gCLinyU9ZfuHekU3Snpc0n2Sbu+V3S7p3pHUEAAAAABWSNHrAP2ipI/ZnpH0HUlvVbfz9Enbb5P0pKQ3j6aKAAAAACbROKeylVWoAxQRj0janXnoxpWtDgAAAIAqCK3hDhAAAAAAXKwTdIAAAAAA1AAjQACAldPJXCGhUSS4cwLlfhYAAMaEDhAAAACAgYWsdqGr6kwWOkAAAAAASmENEAAAAIBaYA0QAAAAgBqx2sEUOADAqFQhGIHAAwDAhKMDBAAAAGBgIalDCAIAAACAumANEAAAAIBaiKjmGqDq1RgAAAAASmIECAAAAEApHabAAQAAAKiD7nWAqjehjA4QAAAAgBKquQaIDhAAAACAgVU1Brt6NQYAAACAkhgBAoCqaFTgM6tcHTud1a8HAGBVtIMQBAAAAAA1EDIhCAAAAADqo0MIAgAAAIA6qGoMdvVqDAAAAAAlMQIEAKupCkEGK22Yn5kABQCYWCETggAAAACgPqp4HSA6QAAAAAAGFiG1KxiCUL0aAwAAAEBJjAABAAAAKMHqiDVAAFBPdQw3WA2DvK4EJgDAqgpVcwocHSAAAAAApVTxOkB0gAAAAAAMLGR1KhiDXb0uGwAAAACUxAgQAAAAgFKYAgcAa92khR14wqYeRIzv2LlzQzACAIxMSOoQggAAAACgHqw2MdgAAAAA6qCqI0DVqzEAAAAAlMQIEAAAAIBSmAIHAFU1znCDXJBBrqxR8I9M0f0NIhduUDTwoFPwuaMIUCh6XglLAICBRZgpcAAAAADqox2N0rfl2L7a9hdtH7L9mO13Zrax7T+wfdj2o7Zfs9x+GQECAAAAMIlakn45Ih62vUnSQdsPRMTjS7Z5g6RdvdtrJX2w97UvRoAAAAAADCwkdeTSt2X3H3E0Ih7ufX9G0iFJV1602R5JH42ur0h6se0dl9ovI0AAAAAASnChqWwrciT7GknXSnroooeulPTUkvtHemVH++2rUAfI9j9IOiOpLakVEbttb5H0fyVdI+kfJP2HiHi2yP4AYKxGEXgwTJBBpj7RzNSx2UzLptKymMo8N1OX6BOM4IIBBW5lggNa7bSsnZa5nXluLoigaIDCpcqLyLUJghEA4JK61wEaKmRnq+0DS+7vi4h9F29ke6OkT0t6V0ScvvjhPlXra5ARoJ+IiJNL7u+V9GBE3Gl7b+/+ewbYHwAAAIAKaw+3ouZkROy+1Aa2p9Xt/HwsIj6T2eSIpKuX3L9K0tOX2ucwNd4jaX/v+/2Sbh1iXwAAAADwPbYt6cOSDkXE+/tsdp+k/9xLg/sxSc9FRN/pb1LxEaCQ9Oe2Q9L/7g1NbX9h5xFx1PYVfSp+h6Q7JGnnzp0FDwcAAABgkoU87BS45Vwv6eckfcP2I72y90naKUkRcbek+yXdLOmwpHlJb11up0U7QNdHxNO9Ts4Dtr9VtNa9ztI+Sdq9e/cIrnIHAAAAYBw6IwyVjogvK7/GZ+k2Ientg+y3UAcoIp7ufT1u+7OSrpN0zPaO3ujPDknHBzkwAAAAgOqKkNqjHQEaiWU7QLY3SGpExJne9/9O0m+pO9/udkl39r7eO8qKAsCyRpHudrE+yWkqmNoWs9OFytrr0rLOXLq/VqasM5PWsTOVS6lLi7oVSosarbSwsZiWTZ1PE98ambLmuQtpdRaKleVS5brlmdS2lU6GyyEtDkCNjXgK3EgUGQHaLumz3TVImpL0JxHxZ7a/KumTtt8m6UlJbx5dNQEAAABgeMt2gCLiO5JenSk/JenGUVQKAAAAwGTrhiCszoVQV9Ig1wECAAAAgO9pXzqjYCLRAQIAAAAwsNDaXQMEAJNnpQMPcuEGubKpNHRAkmJuJinrrE/LWptmk7LFF6VvxYub0p9vcWNan9b6TOBBelh1ctXu9xJm1vQ3MrkDjcW0bGo+/VlmzqZBBDNn0tdh5rlWur8zC+lx5zMHluTzmfJWpuK5YISVDksgGAEAJhYdIAAAAAAlsAYIAAAAQI10WAMEAAAAoA7W7IVQAQAAACCHKXAAMAqrEXjQTI8RM9Np2fp08b4ktTbPJWXnt6ZpBOcuT4+zsCWtz+KmdFF+e31a1plLF9vHVGYBfjOzyL/fh3a5PIB2urFb6c/SOJ9u15xPy2bOpH9+Zp9JkxrWnUrPwdzJfAjC1OnzaR3n0xAFL15In9zOvGYEIwDAmkQHCAAAAMDAuhdCZQocAAAAgJogBAEAAABALVT1QqjVW7UEAAAAACUxAgRgbcsFHkyli+1zgQedzeuSsoWtaZkkPb89fTs9tz099vkt6cL61uZWusN17aSoOZuWTTfThfWNRnoMu1iZJEXm07xcWaeTlrXb6edqiwvp6714Li07vyUtW8iUrcsETkjShmOZYIWT6fMbp88lZdlghFb6eg8VjAAAaxApcAAAAADqIQhBAAAAAFATIUIQAAAAANRIFUeAqjdpDwAAAABKYgQIwGRpDPG5TNHAg9mZpKx92fqk7Nz2dLH92R9I9ydJ8zvSssXLM4voN6aL7adn0xCE6elMCEIm8KA5QLhB0e1ygQdFt2vnymbTc3phXfo6tjakf5Ke35SGU1zYmG8jrfXp8zfOpud13VT6/Oaz80mZtZg5yBDBCP3adic9rwBQBVWNwaYDBAAAAKAUOkAAAAAAaiFEChwAAACAGqliChwhCAAAAABqgxEgAOMxTNiBlA88aKb7jJl0EX0u8GB+x7qk7PTOdKH+/I78gvfW5WmQwdSGNPBgZjYtm8qFGzRyZcOFGxRX7Pm5EITcH5V25mfJ/cytTPDD4nS63ULmnEpSeyY9eq4sGum5TluE1HwmEzCRCzxoZ0IMigYjSPnfBYIRAFRBsAYIAAAAQE2QAgcAAACgVqrYAWINEAAAAIDaYAQIAAAAwMCIwQaA1ZQJQcgFHnQ2pwvez22fS8qygQdXpgvZW1vTEANJmt24kJTNzGSCEQqGGwwfZPD9GkPuL/cHrmgdp5rpdrmfORv8kAlLWJxKwxIkaaE5m5TNN3KBCem5didtE+tb6bEbz2XO1fnF9BCDhCAAQIXlAnEmHR0gAAAAAKVU8TpAdIAAAAAADCwqGoNNCAIAAACA2mAECAAAAEAprAECgJzcle4HkQk80FS6kD3Wp4vgF7amIQhnfyATeLCjWOBBLuxAkmZzgQfNdLF+I/OjFA0TGDbIYBjDHLtogEIz00zs9DX0TJ8DbUyLcmdrvpMGIzQX0zbRXEjbztyFTH1amVCGC5nXa5BghNzvTCcNZQCA8SIFDgAAAECNMAIEAAAAoBZChCAAAAAAwERjBAgAAADA4KKa132mAwRgZY0i8CCzOj7m0pXwrc1zSdnz29O3ufkd6SFal6chBrnAg1zYgZQPPGg2Jj/cYDUU/flyS/yb2ZkVmdABScqFI+SCEdrpTucX03YyNZ8pez5tY9OLaZtwLrCg3SfEoOh/DwQjAJhAXAgVAAAAQC2EqhmCwBogAAAAALXBCBAAAACAErgOEAAAAIAaqWIIQuEpcLabtr9m+09797fYfsD2E72vl42umgAAAAAmTYRL38ZlkBGgd0o6JGlz7/5eSQ9GxJ229/buv2eF6wegbrIpcM2kqLM+jfw6vzUtO7c93d/i5WmK2NSGC0nZTCbxLZf2JpH4thJyr01uakX/1zo9N5FJhmtnzvXi5enngefOpu1u9nS6w+Z8pmwhPYY6fepdxY9PAUDdt681G4Jg+ypJPyPpQ0uK90ja3/t+v6RbV7ZqAAAAALCyio4A/Z6kX5G0aUnZ9og4KkkRcdT2Fbkn2r5D0h2StHPnziGqCgAAAGCSVDEEYdkRINu3SDoeEQfLHCAi9kXE7ojYvW3btjK7AAAAADCButPgyt3GpcgI0PWS3mj7Zklzkjbb/mNJx2zv6I3+7JB0fJQVBQAAADBZqrgGaNkOUES8V9J7Jcn2DZLeHRE/a/t3JN0u6c7e13tHWE8Ak6gxgmspN9I30pidTspam2aTsnOZheznt2Q+YtqYCTyYTcumGp0i1euLwIPhFQ1G6G6bluXOYe5cn9uYazuZYIRMG5t9Nm2LjfnFpMztfICG0ioWl/sd7AyzQwAoLjTeNLeyhvnv5U5Jr7f9hKTX9+4DAAAAwMQa6EKoEfFXkv6q9/0pSTeufJUAAAAAVEEV5zoM1AECAAAAAElSRa8DRAcIAAAAQDkVHAKiAwSstqLBAWt9IbP7fGKUeX1yIQiLL0rfvha2pPtsbW4lZdOzadlUM329m430Xd19gg0mKfBgta7JMK6fud9xc78xzcyvW+5cT2XaxIXNuTaWluXa4vR30zbr82kwQveBTM3HmQ+7GkYRoAJgzbF9j6QXLsnzyszjN6gbxPb3vaLPRMRvLbdfOkAAAAAAShnxFLiPSLpL0kcvsc1fR8Qtg+yUDhAAAACAUkY5YB0RX7J9zUrvlzFoAAAAAAMLdUeAyt4kbbV9YMntjhLVeJ3tr9v+vO1XFHkCI0AAAAAABheShpsCdzIidg/x/IclvTQiztq+WdLnJO1a7kl0gIBJNYpFwpMUrNAnBCEyq9bb6zIhCJvS7RY3Zcbh17WTounptKzZSF+bfoEH47Ja4QZFFa3POAMicucwd65zbeJCpu0sbmpmytK2OJdpsz6b/512K/M6TlIIAoEFACZURJxe8v39tj9ge2tEnLzU8+gAAQAAAChlnJ/X2H6JpGMREbavU3d5z6nlnkcHCAAAAEA5I+wA2f64pBvUXSt0RNJvSJqWpIi4W9KbJP2C7Zakc5Jui1i+S0YHCAAAAEAJHmkMdkS8ZZnH71I3JnsgdIAAAAAAlDNBSxaLogME1EluMXPRYISVXgjd6POJUTNdZN6Zyyw835g+v70+s+B9NhN40Ex/5maj2Dv4ai3on7TAg2HkfpZRvI65feaOnTvX2TaRaTu5NpZri7k2m2vbkqRGKy0bJq9kkn7PAWAC0QECAAAAMLjQSKfAjQodIAAAAADlMAUOAAAAQH1UbwSIyb4AAAAAaoMRIKDuVmPRswf4dGgqXSjeyiwob63PLTxPF3pP5xa3r1KQAaoj1yZywQgX5tLtWuvT36Fcm53JtO2+cr8zw1xtkHADAKNSwT+pdIAAAAAAlEMHCAAAAEAthCRS4AAAAADUxTCzc8eFScEAAAAAaoMRIADj0ScYIabSz2U6M5nAg5ncc9NF641G+tGUMwvec2WrpVPB6QPDyv3MjVU6B0XPf67t5NpYZyYNN8i12Vzb7h67fucfwBpSwREgOkAAAAAAyqngh3h0gAAAAACUUsUrS9ABAgAAADC4UCWnwBGCAAAAAKA2GAECMB79Fn5nyjtTmbJ03bnUXNlwg9ValI/RyZ3DoqET2baTaWO5tphrs4O0eQCoBrMGCAAAAECNVPCzQjpAAAAAAMqpYAeINUAAAAAAaoMRIAAAAADlVHAEiA4QgIkSuQXhufWVufHr7LrzlQ1GwNpUuJ0M0RazbRsAqixECAIAAACA+qjiZ4p0gAAAAACUU8EOECEIAAAAAGqDDhAAAACA2mAKHAAAAIBSWAMEAENyZN5Jc2+unUxZ7qmZdJpcWSUnMWPFFG4nQ7TFbNsGgKojBQ4AAABALYQq+fkha4AAAAAA1AYjQAAAAADKWYsjQLbnbP+t7a/bfsz2b/bKt9h+wPYTva+Xjb66AAAAACaFo/xtXIqMAC1I+smIOGt7WtKXbX9e0r+X9GBE3Gl7r6S9kt4zwroCWEv6LQjPlDdambJ25rntooEHxXQyz21UMe6mxnLnsKhs28m0sVxbzLXZQdo8AFRGBd/Clh0Biq6zvbvTvVtI2iNpf698v6RbR1JDAAAAAFghhUIQbDdtPyLpuKQHIuIhSdsj4qgk9b5e0ee5d9g+YPvAiRMnVqreAAAAAMYthriNSaEOUES0I+JHJV0l6Trbryx6gIjYFxG7I2L3tm3bytYTAAAAwAQZZv3POGeUDxSDHRHflfRXkm6SdMz2DknqfT2+4rUDAAAAMLnC5W9jsmwIgu1tki5ExHdtr5P0U5J+W9J9km6XdGfv672jrCiANabPwm+3OklZYzETgrCYe276mU6nUywYIVfmVfp4KhesMMzi/SoYZ5hE0fOfazu5NpZri7k2m2vbvYPnywGgCir4FlYkBW6HpP22m+qOGH0yIv7U9t9I+qTtt0l6UtKbR1hPAAAAABjash2giHhU0rWZ8lOSbhxFpQAAAABMvipeHaLICBAAAAAApOgAAQAAAKiFMae5lUUHCKi7Tp+F2RdrDBQa+f0GWeTdaidFU+czZfPp21fjfLpovd1O693OLHjnzbDecm0i13ZybWxqPm3fuTaba9t9rXQwwmr8ngNARfA3HwAAAEA5jAABAAAAqA06QAAAAADqooprgJjsCwAAAKA2GAEC6qToQuiizx1mwXSnz0dG7XSheCOzoHzmbPr85ny6QH1xoZkeYjYTjNBIf76pZnqMTmaxvCQ1VvgjsNz++h170q30a9NP0den3SkYlpFpOzOZNpZri7k2m2vbkvr/LpQ1Sb/nADCB6AABAAAAKKeCU+DoAAEAAAAYHNcBAgAAAFArdIAAAAAA1AYdIAArZpiFzFXQ50r3bqc/d/PchaRs5sxspix9S1s8ly5kv7AuLZtqZo7bSOvoMY71Fw0TWK2whNUKNxhGZF6Ldidd1H/hQtomlGk7M2cyIQhn0nCDXJvNte1eJfPlk2IU70UEKwAYIzpAAAAAAAZmsQYIAAAAQJ3QAQIAAABQCxVNgWMSLgAAAICJY/se28dtf7PP47b9B7YP237U9muK7JcRIGC1rfVwg6L6LfzOvD5eyIQgPNdKymafSRetn9+SlrU2pG99renMQvZGLhghKZIk5c7quEICqhBOMIx+IQ/5wIO0rNVOT2JrIW0TU6fTtjP7TPra5tpirs32/d2f9BCEUeB9EFg7RvsW9hFJd0n6aJ/H3yBpV+/2Wkkf7H29JEaAAAAAAJQTQ9yW23XElyQ9c4lN9kj6aHR9RdKLbe9Ybr+MAAEAAAAoZchJB1ttH1hyf19E7Bvg+VdKemrJ/SO9sqOXehIdIAAAAADlDNcBOhkRu4d4fm5O9LI1YgocAAAAgCo6IunqJfevkvT0ck9iBAhAebmFzMNe4b2TfnCTW1A+dWYhKVt3ajopW8iEIDy/Kd1ucToTeNDMBDI4DUuQpGbmM6jcYv21HlCw0voFHuS3TctanbQ9Li6k519n07K5Z9JjrzuVnv9cW8yHIIzg3BMmAGCcCq7lGaH7JL3D9ifUDT94LiIuOf1NogMEAAAAoKRRfq5n++OSblB3rdARSb8haVqSIuJuSfdLulnSYUnzkt5aZL90gAAAAACUM8IOUES8ZZnHQ9LbB90vHSAAAAAApVRxZjchCAAAAABqgxEgAJMlMh8ltdOF5435xaRs7mRatm7zXFJ2YWP62c/CTCYYYSo9rmfS6vUqmZQ0G+nPQjBCf0UDD9qd/Hatdhp4sbiY/plrPZ+e69lTaZtYdyw9L7k2lmuLuTabbdsAUHUVfGujAwQAAABgcONPgSuFDhAAAACAgVn5K5FOOtYAAQAAAKgNRoAAAAAAlMMUOAC1l7syfWOAweZsCEK6T59PF55PnT6flG04llkEvz4ta8+kZQvN2bQuG9MiSVI2HCET3pCZK5B5xbKqGpZQNNwgMtt1Mj9yLuxAkhYygQcLZ9NzOHUq3W595rrhG4610udm2liuLeba7NAhCLnfLQAYsyr+aaIDBAAAAKAcOkAAAAAAaqOCHSBCEAAAAADUBiNAAAAAAAYXrAECAAAAUCd0gAAgYxTJcK00Yc3zC0nZ7Mk0MWzj7PqkLJcCN9+YTsrSI7yw07QoMslwU430tWhmXgpnPlIrmqaWM2yC3DDHzsklvrU7aVmrk744i5m0N6lP4tvJ9ByuP5oeZ+PTaeLb7MlzSVmujeXaIolvAOqCESAAAAAA9VHBDhAhCAAAAABqgxEgAAAAAKVUcQrcsiNAtq+2/UXbh2w/ZvudvfItth+w/UTv62Wjry4AAACAiRBD3sakyAhQS9IvR8TDtjdJOmj7AUn/RdKDEXGn7b2S9kp6z+iqCgBLZBaZe/FCUtY4nS5kXzeVfvYTjXWZg6QBCvOddFG9JC20M4v6N6T1mZlNy6aauWCEXFmxvxYrHaDQTy7IICcXbtDOhBu02pnAg4X09W49nz8HU6fSP2m5wIPNT6ahBeuOnU/Kcm0n18aGDjwAgCqr4FvgsiNAEXE0Ih7ufX9G0iFJV0raI2l/b7P9km4dVSUBAAAAYCUMtAbI9jWSrpX0kKTtEXFU6naSbF/R5zl3SLpDknbu3DlMXQEAAABMCGuNrgF6ge2Nkj4t6V0Rcbro8yJiX0Tsjojd27ZtK1NHAAAAAJNoja4Bku1pdTs/H4uIz/SKj9ne0Rv92SHp+KgqCQAAAGDyuILrIJftANm2pA9LOhQR71/y0H2Sbpd0Z+/rvSOpIYC1qd+V7hsFB6Zzb7jtdJ+5RevNZ+eTsvWZQ7gzlz53MQ1GkKT5xfTtdPHy9Gc5tzEtm5ptJWXT0+lC/WYuLCEz9yAXgpDTb7ui4Qa57dq5sky4wYUL6evYWsj8STqbBh7Mnsq3kfVH07KNT6evbS7wINcmsoEHmTY2dAhCv98FAJh0Yx7JKavICND1kn5O0jdsP9Ire5+6HZ9P2n6bpCclvXk0VQQAAACAlbFsBygivqzuGqecG1e2OgAAAACqooohCAOlwAEAAADA99ABAgAAAFAXjAABwLByC8KHCUZopWEC1mJS1nwmfe76ViZ0YGFd9tBT8+nb6bmz6UL/81sygQCb0+deWJcJQZgtFozQaBQLRhgkBCFX1ukUCzxoL2SCI86lZVOn07K5Z9JjrDuWr/eGY2ngwezJc0lZ43Ralg08yLSdoQIPCDsAsBZVsANU+DpAAAAAAFB1jAABAAAAGFwwBQ4AAABAndABAgAAAFAHFiNAADB5igYjZLZrPJeWzV3ILIyXNPX8XFI2e3omKTt3ebr0cmFL+la8uCkNBGivT+tzYS4ti6nMYvtm5nXod4W33B+zdrqxW+nP0jifbjcznyk7k5bNZoIo1p1KX++5k2mIhSRNnT6f1nF+IS3LBR60M6/ZMIEHAICJRQcIAAAAQDkV/FiKbJYAAAzPSURBVLCIDhAAAACAUpgCBwAAAKAeQoQgAAAAAKgPV/Aaz3SAAEy+TubdtTHEdZxz85Uzi+B9Pl1s70yAgiRNL7aSsuZ8GoIw++xsUrb4olwIQvrzLW5MgwNa69PtOjNpgEInLep/Kezcy535sRuZLIKp+fS1nTmbKTuT7nDmufQ1nDqThhg05vMhCLnzlQu8yJ7/lZ7DnmuzAICJQAcIAAAAQDlMgQMAAABQF4QgAAAAAKiHEDHYAAAAAOqDESAAWC2rEYyQK7uQf6d3pj7NhQtJWW4B//R3p5OyuXVpWWcuTTJoZco6M2lYQmcqLVOmSFJ2PnejlRY2FtOyqfNp6EAjU9Y8l742zrxeuTK180EUuSCLVflkksADAKgUOkAAAAAAymEECAAAAEAdWEyBAwAAAFAXEZUMQRhiwjwAAAAAVAsjQADWjqKL0Vc6LEHKL8DvpNs6s4Df59NgBJ/N1LGZBh7MTKVlMZV5rtPEg8iUSZILBkK4lfmZW5mAgtzPnH29ir2Gfc/BSn8KSbgBACyLKXAAAAAA6oMOEAAAAIC6YAQIAAAAQD2E8lOVJxwhCAAAAABqgxEgAPWTW9w+TDCClF+AnyvLrat3WuhWJqCg0SpUFefCDfoEHhRW9OfLKRpkMM4oVQIPAKCc6g0A0QECAAAAUA5rgAAAAADUBxdCBQAAAFAXjvK3Qvu3b7L9bduHbe/NPH6D7edsP9K7/fpy+2QECAAAAMDEsd2U9IeSXi/piKSv2r4vIh6/aNO/johbiu6XESAAAAAAg4shb8u7TtLhiPhORCxK+oSkPcNWmxEgAJCKp4ANmxaXM1SC3JDpbitt0uaCk+4GACNjSR7t+/6Vkp5acv+IpNdmtnud7a9LelrSuyPisUvtlA4QAAAAgHKG+5xpq+0DS+7vi4h9S+7nPuW7uMf1sKSXRsRZ2zdL+pykXZc6KB0gAAAAAONwMiJ2X+LxI5KuXnL/KnVHeb4nIk4v+f5+2x+wvTUiTvbbKWuAAAAAAJTiiNK3Ar4qaZftl9mekXSbpPu+7/j2S9y7Arjt69Tt35y61E4ZAQIAAAAwuOJhBuV2H9Gy/Q5JX5DUlHRPRDxm++d7j98t6U2SfsF2S9I5SbdFXLp3RQcIAAbRb1H9KMIRipi00IFxIvAAAFZZjPzvUETcL+n+i8ruXvL9XZLuGmSfdIAAAAAAlFL0gqaThDVAAAAAAGqDESAAAAAA5VRwKvayI0C277F93PY3l5Rtsf2A7Sd6Xy8bbTUBAAAATJSQ3Cl/G5ciU+A+Iummi8r2SnowInZJerB3HwDqq9MpdsNgir6uvLYAMB4R5W9jsmwHKCK+JOmZi4r3SNrf+36/pFtXuF4AAAAAsOLKhiBsj4ijktT7ekW/DW3fYfuA7QMnTpwoeTgAAAAAEyeGuI3JyFPgImJfROyOiN3btm0b9eEAAAAArBJHlL6NS9kO0DHbOySp9/X4ylUJAAAAQCVUcA1Q2Rjs+yTdLunO3td7V6xGALCWDbNYv1HRS7cRUAAAa1NIquBbfJEY7I9L+htJP2T7iO23qdvxeb3tJyS9vncfAAAAACbasiNAEfGWPg/duMJ1AQAAAFAR1njX8pRVdgocAAAAgLqjAwQAAACgNugAAQBGJhcmMGnBCAQeAEB9rNUQBAAAAABYKxgBAgAAAFAKIQgAAAAA6oMOEAAAAIB6iEp2gFgDBAAAAKA2GAECAAAAMLhQJUeA6AABAAAAKKeCMdh0gAAAAACUQgocAAAAgPqgAwQAGJlGBXJrcnXsVHB+BABgzaIDBAAAAGBwIanDCBAAAACAWqjmdYDoAAEAAAAohw4QAAAAgNqgAwQAWBFVCDwoimAEAMAEoQMEAAAAYHCEIAAAAACoj5CieiP6dIAAAAAAlFPBNUBraJI5AAAAAFwaI0AAAAAABscaIAAAAAC1UsEpcHSAAAAAAJRDBwgAAABAPUQlO0CEIAAAAACoDUaAAAAAAAwuJHW4DhAAAACAuqjgFDg6QAAAAADKoQMEAAAAoB6iktcBIgQBAAAAQG0wAgQAAABgcCFFEIIAAAAAoC4qOAWODhAAAACAcioYgsAaIAAAAAC1wQgQAAAAgMFFcCFUAAAAADVSwSlwdIAAYNwaNZyNnPuZK/gpIgDUXVTwvZsOEAAAAIASopIjQDX82BEAAABAXTECBAAAAGBwoUpeB2ioESDbN9n+tu3DtveuVKUAAAAAVEB0yt/GpPQIkO2mpD+U9HpJRyR91fZ9EfH4SlUOAAAAwGQKSVGzEaDrJB2OiO9ExKKkT0jaszLVAgAAADDRIkY+ArTcjDN3/UHv8Udtv2a5fQ7TAbpS0lNL7h/plV1cqTtsH7B94MSJE0McDgAAAEBdLJlx9gZJL5f0Ftsvv2izN0ja1bvdIemDy+13mA6QM2XJGFhE7IuI3RGxe9u2bUMcDgAAAMAkiU6UvhVQZMbZHkkfja6vSHqx7R2X2ukwHaAjkq5ecv8qSU8PsT8AAAAAVTLaKXBFZpwVmpW21DAx2F+VtMv2yyT9k6TbJP2nSz3h4MGDJ23/4xDHRHlbJZ0cdyWQ4LxMHs7JZOK8TKYVOy/2XSuxG3Tx+zKZljsvL12tiqyUM3r2C38Rn9o6xC7mbB9Ycn9fROxbcr/IjLNCs9KWKt0BioiW7XdI+oKkpqR7IuKxZZ7DHLgxsX0gInaPux74fpyXycM5mUycl8nEeZlMnJfJtBbPS0TcNOJDFJlxNvCstKEuhBoR90u6f5h9AAAAAEBGkRln90l6h+1PSHqtpOci4uildjpUBwgAAAAARqHfjDPbP997/G51B2NulnRY0rykty63XzpA9bFv+U0wBpyXycM5mUycl8nEeZlMnJfJxHkpITfjrNfxeeH7kPT2Qfbp7nMAAAAAYO0bJgYbAAAAACqFDtAaZvt3bH/L9qO2P2v7xUsee6/tw7a/bfunx1nPOrJ9U++1P2x777jrU1e2r7b9RduHbD9m+5298i22H7D9RO/rZeOua93Ybtr+mu0/7d3nnEwA2y+2/ane35ZDtl/HuRkv27/Ue//6pu2P257jnIyH7XtsH7f9zSVlfc8F/4uNDx2gte0BSa+MiFdJ+n+S3itJtl+uborGKyTdJOkDtptjq2XN9F7rP5T0Bkkvl/SW3jnB6mtJ+uWI+BFJPybp7b1zsVfSgxGxS9KDvftYXe+UdGjJfc7JZPh9SX8WET8s6dXqniPOzZjYvlLSf5e0OyJeqe4i8dvEORmXj6j7f9VS2XPB/2LjRQdoDYuIP4+IVu/uV9TNRZekPZI+ERELEfH36qZmXDeOOtbUdZIOR8R3ImJR0ifUPSdYZRFxNCIe7n1/Rt1/5q5U93zs7222X9Kt46lhPdm+StLPSPrQkmLOyZjZ3izpxyV9WJIiYjEivivOzbhNSVpne0rSenWvf8I5GYOI+JKkZy4q7ncu+F9sjOgA1cd/lfT53vdXSnpqyWNHemVYHbz+E8j2NZKulfSQpO0vXEOg9/WK8dWsln5P0q9I6iwp45yM3w9KOiHpj3rTEz9ke4M4N2MTEf8k6XclPSnpqLrXP/lzcU4mSb9zwf8CY0QHqOJs/0Vv3u/Ftz1LtvlVdaf6fOyFosyuiANcPbz+E8b2RkmflvSuiDg97vrUme1bJB2PiIPjrgsSU5JeI+mDEXGtpOfF1Kqx6q0n2SPpZZJ+QNIG2z873lqhIP4XGCOuA1RxEfFTl3rc9u2SbpF0Y/xL5vkRSVcv2ewqdYfMsTp4/SeI7Wl1Oz8fi4jP9IqP2d4REUdt75B0fHw1rJ3rJb3R9s2S5iRttv3H4pxMgiOSjkTEQ737n1K3A8S5GZ+fkvT3EXFCkmx/RtK/EedkkvQ7F/wvMEaMAK1htm+S9B5Jb4yI+SUP3SfpNtuztl8maZekvx1HHWvqq5J22X6Z7Rl1F0HeN+Y61ZJtq7ue4VBEvH/JQ/dJur33/e2S7l3tutVVRLw3Iq6KiGvU/d34y4j4WXFOxi4i/lnSU7Z/qFd0o6THxbkZpycl/Zjt9b33sxvVXcvIOZkc/c4F/4uNERdCXcNsH5Y0K+lUr+grEfHzvcd+Vd11QS11p/18Pr8XjELv0+3fUzex556I+J9jrlIt2f63kv5a0jf0L+tN3qfuOqBPStqp7j8Yb46Iixe2YsRs3yDp3RFxi+3LxTkZO9s/qm44xYyk70h6q7ofpnJuxsT2b0r6j+r+Pf+apP8maaM4J6vO9scl3SBpq6Rjkn5D0ufU51zwv9j40AECAAAAUBtMgQMAAABQG3SAAAAAANQGHSAAAAAAtUEHCAAAAEBt0AECAAAAUBt0gAAAAADUBh0gAAAAALVBBwgAAABAbfx/VsQ2+H59nOoAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 1152x432 with 2 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "init()\n",
+    "time_loop(1000)\n",
+    "plot_ρs()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.7.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/doc/notebooks/08_tutorial_shanchen_twocomponent.ipynb b/doc/notebooks/08_tutorial_shanchen_twocomponent.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..81329d401fceb34b3b35ad69af670d2b1907827c
--- /dev/null
+++ b/doc/notebooks/08_tutorial_shanchen_twocomponent.ipynb
@@ -0,0 +1,392 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Shan-Chen Two-Component Lattice Boltzmann"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from lbmpy.session import *\n",
+    "from lbmpy.updatekernels import create_stream_pull_with_output_kernel\n",
+    "from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This is based on section 9.3.3 of Krüger et al.'s \"The Lattice Boltzmann Method\", Springer 2017 (http://www.lbmbook.com).\n",
+    "Sample code is available at [https://github.com/lbm-principles-practice/code/](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "N = 64\n",
+    "omega_a = 1.\n",
+    "omega_b = 1.\n",
+    "g_aa = 0.\n",
+    "g_ab = g_ba = 6.\n",
+    "g_bb = 0.\n",
+    "rho0 = 1."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Data structures"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dh = ps.create_data_handling((N, N), periodicity=True, default_target='cpu')\n",
+    "\n",
+    "method_a = create_lb_method(relaxation_rate=omega_a, compressible=True)\n",
+    "method_b = create_lb_method(relaxation_rate=omega_b, compressible=True)\n",
+    "\n",
+    "src_a = dh.add_array('src_a', values_per_cell=len(method_a.stencil))\n",
+    "dst_a = dh.add_array_like('dst_a', 'src_a')\n",
+    "\n",
+    "src_b = dh.add_array('src_b', values_per_cell=len(method_b.stencil))\n",
+    "dst_b = dh.add_array_like('dst_b', 'src_b')\n",
+    "\n",
+    "ρ_a = dh.add_array('rho_a')\n",
+    "ρ_b = dh.add_array('rho_b')\n",
+    "u_a = dh.add_array('u_a', values_per_cell=dh.dim)\n",
+    "u_b = dh.add_array('u_b', values_per_cell=dh.dim)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Force & combined velocity"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The force between the two components is\n",
+    "$\\vec{F}_k(\\vec{x})=-\\psi(\\rho_k(\\vec{x}))\\sum\\limits_{k^\\prime\\in\\{A,B\\}}g_{kk^\\prime}\\sum\\limits_{i=1}^{19}w_i\\psi(\\rho_{k^\\prime}(\\vec{x}+\\vec{c}_i))\\vec{c}_i$\n",
+    "for $k\\in\\{A,B\\}$\n",
+    "and with \n",
+    "$\\psi(\\rho)=\\rho_0\\left[1-\\exp(-\\rho/\\rho_0)\\right]$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def psi(dens):\n",
+    "    return rho0 * (1. - sp.exp(-dens / rho0));"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "zero_vec = sp.Matrix([0] * dh.dim) \n",
+    "\n",
+    "force_a = zero_vec\n",
+    "for factor, ρ in zip([g_aa, g_ab], [ρ_a, ρ_b]):\n",
+    "    stencil, weights = method_a.stencil, method_a.weights\n",
+    "    force_a += sum((psi(ρ[d]) * w_d * sp.Matrix(d)\n",
+    "                    for d, w_d in zip(stencil, weights)), \n",
+    "                   zero_vec) * psi(ρ_a.center) * -1 * factor\n",
+    "\n",
+    "force_b = zero_vec\n",
+    "for factor, ρ in zip([g_ba, g_bb], [ρ_a, ρ_b]):\n",
+    "    stencil, weights = method_b.stencil, method_b.weights\n",
+    "    force_b += sum((psi(ρ[d]) * w_d * sp.Matrix(d)\n",
+    "                    for d, w_d in zip(stencil, weights)), \n",
+    "                   zero_vec) * psi(ρ_b.center) * -1 * factor"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The barycentric velocity, which is used in place of the individual components' velocities in the equilibrium distribution and Guo force term, is\n",
+    "$\\vec{u}=\\frac{1}{\\rho_a+\\rho_b}\\left(\\rho_a\\vec{u}_a+\\frac{1}{2}\\vec{F}_a+\\rho_b\\vec{u}_b+\\frac{1}{2}\\vec{F}_b\\right)$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "u_full = [(ρ_a.center * u_a(i) + force_a[i]/2 + \n",
+    "           ρ_b.center * u_b(i) + force_b[i]/2) / (ρ_a.center + ρ_b.center)\n",
+    "          for i in range(dh.dim)]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Kernels"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "stream_a = create_stream_pull_with_output_kernel(method_a, src_a, dst_a, {'density': ρ_a, 'velocity': u_a})\n",
+    "stream_b = create_stream_pull_with_output_kernel(method_b, src_b, dst_b, {'density': ρ_b, 'velocity': u_b})\n",
+    "\n",
+    "# TODO use method above\n",
+    "collision_a = create_lb_update_rule(relaxation_rate=omega_a, \n",
+    "                                    compressible=True,\n",
+    "                                    velocity_input=u_full, density_input=ρ_a,\n",
+    "                                    force_model='guo', \n",
+    "                                    force=force_a,\n",
+    "                                    kernel_type='collide_only',\n",
+    "                                    optimization={'symbolic_field': src_a})\n",
+    "\n",
+    "collision_b = create_lb_update_rule(relaxation_rate=omega_b,\n",
+    "                                    compressible=True,\n",
+    "                                    velocity_input=u_full, density_input=ρ_b,\n",
+    "                                    force_model='guo', \n",
+    "                                    force=force_b,\n",
+    "                                    kernel_type='collide_only',\n",
+    "                                    optimization={'symbolic_field': src_b})\n",
+    "\n",
+    "opts = {'cpu_openmp': False}\n",
+    "stream_a_kernel = ps.create_kernel(stream_a, **opts).compile()\n",
+    "stream_b_kernel = ps.create_kernel(stream_b, **opts).compile()\n",
+    "collision_a_kernel = ps.create_kernel(collision_a, **opts).compile()\n",
+    "collision_b_kernel = ps.create_kernel(collision_b, **opts).compile()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Initialization"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "init_a = macroscopic_values_setter(method_a, velocity=(0, 0), \n",
+    "                                   pdfs=src_a.center_vector, density=ρ_a.center)\n",
+    "init_b = macroscopic_values_setter(method_b, velocity=(0, 0), \n",
+    "                                   pdfs=src_b.center_vector, density=ρ_b.center)\n",
+    "init_a_kernel = ps.create_kernel(init_a, ghost_layers=0).compile()\n",
+    "init_b_kernel = ps.create_kernel(init_b, ghost_layers=0).compile()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def init():\n",
+    "    dh.fill(ρ_a.name, 0.1, slice_obj=ps.make_slice[:, :0.5])\n",
+    "    dh.fill(ρ_a.name, 0.9, slice_obj=ps.make_slice[:, 0.5:])\n",
+    "\n",
+    "    dh.fill(ρ_b.name, 0.9, slice_obj=ps.make_slice[:, :0.5])\n",
+    "    dh.fill(ρ_b.name, 0.1, slice_obj=ps.make_slice[:, 0.5:])\n",
+    "\n",
+    "    dh.run_kernel(init_a_kernel)\n",
+    "    dh.run_kernel(init_b_kernel)\n",
+    "    \n",
+    "    dh.fill(u_a.name, 0.0)\n",
+    "    dh.fill(u_b.name, 0.0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Timeloop"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sync_pdfs = dh.synchronization_function([src_a.name, src_b.name])\n",
+    "sync_ρs = dh.synchronization_function([ρ_a.name, ρ_b.name])\n",
+    "\n",
+    "def time_loop(steps):\n",
+    "    dh.all_to_gpu()\n",
+    "    for i in range(steps):\n",
+    "        sync_ρs()\n",
+    "        dh.run_kernel(collision_a_kernel)\n",
+    "        dh.run_kernel(collision_b_kernel)\n",
+    "        \n",
+    "        sync_pdfs()\n",
+    "        dh.run_kernel(stream_a_kernel)\n",
+    "        dh.run_kernel(stream_b_kernel)\n",
+    "        \n",
+    "        dh.swap(src_a.name, dst_a.name)\n",
+    "        dh.swap(src_b.name, dst_b.name)\n",
+    "    dh.all_to_cpu()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_ρs():\n",
+    "    plt.subplot(1,2,1)\n",
+    "    plt.title(\"$\\\\rho_A$\")\n",
+    "    plt.scalar_field(dh.gather_array(ρ_a.name), vmin=0, vmax=2)\n",
+    "    plt.colorbar()\n",
+    "    plt.subplot(1,2,2)\n",
+    "    plt.title(\"$\\\\rho_B$\")\n",
+    "    plt.scalar_field(dh.gather_array(ρ_b.name), vmin=0, vmax=2)\n",
+    "    plt.colorbar()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Run the simulation\n",
+    "### Initial state"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA48AAAF0CAYAAACDowz8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X2wZHd93/n3RyOJcQasBwRYJQkjEtkWJEjgKRmXUubJCJkQpNTCrhTiTCi5VEnAwV7vJsK7gV1sV5nNlolTAWPFKJJTBkGEMbMuGaEVsCSxwRqBAnoAa5BZNB5FYzESYBlJq5nv/tFnTOtq7nTPvbf7dn/n/ao61X0euvt34Eoffc/v/H4nVYUkSZIkSUdy3GY3QJIkSZK0+CweJUmSJEkTWTxKkiRJkiayeJQkSZIkTWTxKEmSJEmayOJRkiRJkjSRxaMkSZIkaSKLR0mSJEnSRBaPai/JM5JcneShJPuS/Pxmt0mSpGOZ2SwtJ4tHHQt+D/ga8APAZcD/meQHNrdJkiQd08xmaQlZPKq1JK8DqKp3V9VjVfUp4M+AHxr2vzDJgSRnbmY7JUk6Vhwpm5P8xyT/ZVg+lWTr5rZW0jiLR3X3euDjh1aSHAecBDwwbPoXwH8Azp1/0yRJOiYdKZt/CHhZVV0IfBt40aa0UNJhWTyqux8Dvjm2/krgwar6apIXAfcDN2HxKEnSvBw2m4E/BY6rqieSPB04HfiTTWifpFVYPKqtJCcA5wBvSLI1yQuB9zHqbQT4eeDdwF1YPEqSNHMTsvlc4AeSfAa4B3hvVT28aY2V9BTHb3YDpBk6F/g6cAejW2H2Ab9cVTckOR+4EPggsGVYJEnSbB0pm98E/EZVvSPJScDngN/etJZKegqLR3X2IuDuqvqXwL9cse8q4Meq6iGAJH8878ZJknQMOlI2/03g1uH9KcC35tkwSZNZPKqz84C7V25M8qPAdw8VjoNHkzyzqr658nhJkrRhDpvNgxcCFyV5C3AQ+Mdza5WkqTjmUZ29CPjKyo1VdVtVvXnFtp+wcJSOLMlZST6d5O4kdyZ522GOSZJ/k2R3ki8lecnYvh1J7hmWHfNtvaQFcdhsBqiq11fVj1bVq6rq1VV1+5zbJi2deWdzqmqjz0GS1FCS04HTq+oLSZ4B3AZcWlV3jR3zWuBngdcymlHx16vqx5KcCuwCtgM1fPZHV9wBIEmSjsK8s9meR0nSVKrq/qr6wvD+O4xuPTtjxWGXAL9dI58DTh6C7TXAzVW1fwilm4GL59h8SZLamXc2WzxKko5akucBLwY+v2LXGcB9Y+t7hm2rbZckSRtgHtk81wlztmzbViecfOo8f1KSFtL/9/B+DjzySDbyO1/zim31zf0H1vz527702J3Ao2Obrq6qq1ceNzy8+6PAz1XVt1fuPsxX1xG2a5OdeNzW+r7jnrHZzZCkTffdg9/h8YOPms1HMNfi8YSTT+Wsf/rz8/xJSVpI973vPRv+nQ/uP8DnbzpzzZ8/4fSvPVpV2490zPCA748Cv1NVv3uYQ/YAZ42tnwnsHba/fMX2z6y5sdow33fcM/jxk/7eZjdDkjbdH33rYxv+nd2y2dtWJUlTSRLgA4ye0fZrqxy2E/iHw8xuLwW+VVX3AzcxmoL/lCSnABcN2yRJ0hrNO5t9zqMktVEcqIOz/IELgZ8Gvpzk0BT6vwg8F6Cq3g/cyGg2t93AXwJvHvbtT/JLfO8B4O+qqv2zbKwkSZuvVzZbPEpSEwUcnOEwwqr6zxx+fMT4MQW8ZZV91wDXzKBpkiQtpG7ZbPEoSY0cZKZXNyVJ0lHqlM2OeZQkSZIkTWTPoyQ1URQHyqdfSJK0KLpls8WjJDUyy3EVkiTp6HXKZotHSWqigAONAkqSpGXXLZstHiWpkU5XNyVJ6qBTNjthjiRJkiRpInseJamJglaD8iVJWnbdstniUZIa6fMkKUmSeuiUzVPdtprk5CQ3JPlKkruT/HiSU5PcnOSe4fWUWTdWkrS6ojiwjkXLxWyWpMXXLZunHfP468AnqupHgPOAu4GrgFuq6hzglmFdkrRZCg6sY9HSMZsladE1y+aJxWOS7wd+AvgAQFU9XlUPA5cA1w2HXQdcOqtGSpKk7zGbJUmbYZqex+cDfw78+yRfTPJbSbYBz6mq+wGG12cf7sNJrkyyK8muA488smENlyQ9WTEaV7HWRUtlw7L58Xp0fq2WpGNMt2yepng8HngJ8BtV9WLgEY7iNpiqurqqtlfV9i3btq2xmZKkycKBdSxaKhuWzSdm66zaKElqls3TFI97gD1V9flh/QZGgfVAktMBhtd9s2miJGkaBRystS9aKmazJC2Bbtk8sXisqv8G3Jfkh4dNrwLuAnYCO4ZtO4CPz6SFkiTpScxmSdJmmPY5jz8L/E6SE4F7gTczKjw/kuQK4BvAG2fTREnStBbxFhfNjNksSUugUzZPVTxW1e3A9sPsetXGNkeStFZFr4DSkZnNkrT4umXztD2PkqQlcLD6BJQkSR10ymaLR0lqotvVTUmSll23bJ5mtlVJkiRJ0jHOnkdJaqIIB7wmKEnSwuiWzRaPktRIp3EVkiR10CmbLR4lqYlu4yokSVp23bLZ4lGS2ggHqs+tMZIkLb9e2dznTCRJkiRJM2PPoyQ1UcBBrwlKkrQwumWzxaMkNdJpXIUkSR10ymaLR0lqoqrXuApJkpZdt2zucyaSJEmSpJmx51GSGjnY6NYYSZI66JTNFo+S1MToWVLeUCJJ0qLols0Wj5LUxmzHVSS5BngdsK+q/uZh9v/PwJuG1eOBc4FnVdX+JF8HvgMcAJ6oqu0za6gkSQujVzb3KYMl6Rh3aDrwtS5TuBa4eNXfr/pXVXV+VZ0PvB34f6pq/9ghrxj2WzhKko4J3bLZ4lGSNJWq+iywf+KBI5cDH5phcyRJOubNO5stHiWpkQOVNS/AaUl2jS1XrqUNSf4ao6ugHx3bXMAnk9y21u+VJGkZdcpmxzxKUhNF1jso/8ENuqX07wL/ZcVtMRdW1d4kzwZuTvKV4WqpJEltdctmi0dJauTgYjyI+DJW3BZTVXuH131JPgZcAFg8SpLa65TNC3EmkqT1OzQd+FqXjZDkJOBlwMfHtm1L8oxD74GLgDs25AclSVpg3bLZnkdJ0lSSfAh4OaPxF3uAdwInAFTV+4fD/h7wyap6ZOyjzwE+lgRGufPBqvrEvNotSVJX885mi0dJaqL4q8H1s/n+qsunOOZaRtOGj2+7FzhvNq2SJGlxdctmi0dJamTKZ0JJkqQ56ZTNFo+S1EQVHFiMQfmSJIl+2dznTCRJkiRJM2PPoyS1EQ4yu3EVkiTpaPXKZotHSWqi6HVrjCRJy65bNls8SlIjG/VMKEmStDE6ZbPFoyQ1UYSDM5wOXJIkHZ1u2dynDJYkSZIkzYw9j5LUSKdbYyRJ6qBTNls8SlITBRxsNChfkqRl1y2bLR4lqY1woNF04JIkLb9e2WzxKElNdLu6KUnSsuuWzX3ORJIkSZI0M/Y8SlIjnW6NkSSpg07ZbPEoSU1UpdWtMZIkLbtu2WzxKEmNHGgUUJIkddApm/uciSRJkiRpZux5lKQmCjjYaFyFJEnLrls2WzxKUhtpdWuMJEnLr1c2T1U8Jvk68B3gAPBEVW1PcirwYeB5wNeB/76qHppNMyVJk4yeJdXn6qaOzGyWpMXXLZuPpgx+RVWdX1Xbh/WrgFuq6hzglmFdkrSJDnDcmhctJbNZkhZcp2xeT4suAa4b3l8HXLr+5kiSpHUwmyVJMzNt8VjAJ5PcluTKYdtzqup+gOH12Yf7YJIrk+xKsuvAI4+sv8WSpMMqwsFa+6KlsyHZ/Hg9OqfmStKxp1s2TzthzoVVtTfJs4Gbk3xl2h+oqquBqwG2nnFWraGNkqQpHVzAW1w0MxuSzScd/yyzWZJmqFM2T1U8VtXe4XVfko8BFwAPJDm9qu5Pcjqwb4btlCRNUAUHFvAqpWbDbJakxdctmyeWwUm2JXnGoffARcAdwE5gx3DYDuDjs2qkJGk6nW6N0erMZklaHp2yeZqex+cAH0ty6PgPVtUnktwKfCTJFcA3gDfOrpmSJGmM2SxJmruJxWNV3Qucd5jt3wReNYtGSZKO3mhQfp9xFVqd2SxJy6FbNk87YY4kaQkcYPFucZEk6VjWKZstHiWpiYKFHB8hSdKxqls29+lDlSRJkiTNjD2PktRGr3EVkiQtv17Z3OdMJEkcJGteJklyTZJ9Se5YZf/Lk3wrye3D8o6xfRcn+WqS3Umu2sBTliRpoXXKZnseJamJOTyI+Frg3wK/fYRj/lNVvW58Q5ItwHuBVwN7gFuT7Kyqu2bVUEmSFkG3bLZ4lKRGZnlrTFV9Nsnz1vDRC4Ddw+MlSHI9cAlg8ShJaq9TNnvbqiTpkNOS7BpbrlzDd/x4kv+a5A+SvHDYdgZw39gxe4ZtkiTpyBYqm+15lKQmRg8iXtetMQ9W1fZ1fP4LwA9W1V8keS3we8A5cNhBG7WO35EkaSl0y2Z7HiWpkVkOyp+kqr5dVX8xvL8ROCHJaYyuZp41duiZwN51/6AkSUugUzbb8yhJTWz2g4iT/ADwQFVVkgsYXaD8JvAwcE6Ss4E/Ay4D/v6mNVSSpDnpls0Wj5KkqST5EPByRuMv9gDvBE4AqKr3A28A/kmSJ4DvApdVVQFPJHkrcBOwBbimqu7chFOQJKmVeWezxaMkNTLjGd0un7D/3zKaLvxw+24EbpxFuyRJWmSdstniUZK6qHUPypckSRupWTZbPEpSEwUbMrhekiRtjG7ZbPEoSY10uropSVIHnbLZR3VIkiRJkiay51GSmtjs6cAlSdKTdctmi0dJaqRTQEmS1EGnbLZ4lKQmil4zukmStOy6ZbPFoyQ10mlGN0mSOuiUzU6YI0mSJEmayJ5HSeqieo2rkCRp6TXLZotHSWqi24xukiQtu27ZbPEoSY10CihJkjrolM2OeZQkSZIkTWTPoyQ10W06cEmSll23bLZ4lKRGqlFASZLUQadstniUpEY6PUtKkqQOOmWzxaMkNVHNpgOXJGnZdctmJ8yRJEmSJE1kz6MkNdJpXIUkSR10ymaLR0lqo9eMbpIkLb9e2WzxKEmNdLq6KUlSB52y2eJRkpooeg3KlyRp2XXLZifMkSRJkiRNZM+jJHVRoynBJUnSgmiWzRaPktRIpwcRS5LUQadstniUpCaKXoPyJUladt2y2TGPkiRJkqSJ7HmUpDZ6PUtKkqTl1yubLR4lqZFOg/IlSeqgUzZPfdtqki1Jvpjk94f1s5N8Psk9ST6c5MTZNVOSNI2qrHnR8jGbJWnxdcrmoxnz+Dbg7rH1dwPvqapzgIeAKzayYZKko1PVK6A0FbNZkhZYt2yeqnhMcibwd4DfGtYDvBK4YTjkOuDSWTRQkiQ9ldksSZq3aXse/zXwz4GDw/ozgYer6olhfQ9wxuE+mOTKJLuS7DrwyCPraqwk6cgOVta8TJLkmiT7ktyxyv43JfnSsPxhkvPG9n09yZeT3J5k1wae8rFsQ7L58Xp09i2VpGNYp2yeWDwmeR2wr6puG998mEMPOxS0qq6uqu1VtX3Ltm3TtEmStEaj22PWtkzhWuDiI+z/U+BlVfUi4JeAq1fsf0VVnV9V29dybvqejczmE7N1Jm2UJI10yuZpZlu9EHh9ktcCW4HvZ3S18+Qkxw9XOM8E9k7zg5Kk2Znl+Iiq+myS5x1h/x+OrX6OUTZoNsxmSVoSnbJ5Ys9jVb29qs6squcBlwGfqqo3AZ8G3jActgP4+HoaIklan2LtA/KHYDvt0K2Mw3LlOppzBfAHT2oefDLJbev8XmE2S9Ky6JbN63nO478Ark/yy8AXgQ+s47skSZvvwY24pTTJKxgF1N8e23xhVe1N8mzg5iRfqarPrve39BRmsyT1slDZfFTFY1V9BvjM8P5e4IKj+bwkabY2+znESV7EaPbPn6qqbx7aXlV7h9d9ST7GKD8sHjeA2SxJi61TNh/Ncx4lSYtsk58lleS5wO8CP11VfzK2fVuSZxx6D1wEHHZWOEmSWmmWzeu5bVWStGhmeHkzyYeAlzMaf7EHeCdwAkBVvR94B6PHRbxv9MhBnhhutXkO8LFh2/HAB6vqE7NrqSRJC6RRNls8SpKmUlWXT9j/M8DPHGb7vcB5T/2EJElaj3lns8WjJDUyy+nAJUnS0euUzRaPktTIlA8UliRJc9Ipmy0eJamJotfVTUmSll23bLZ4lKQuCmgUUJIkLb1m2eyjOiRJkiRJE9nzKEmNdBpXIUlSB52y2eJRkjppFFCSJLXQKJstHiWpjbQalC9J0vLrlc0Wj5LUSaOrm5IktdAom50wR5IkSZI0kT2PktRF9XqWlCRJS69ZNls8SlInjW6NkSSphUbZbPEoSa30ubopSVIPfbLZMY+SJEmSpInseZSkThrdGiNJUguNstniUZI6aRRQkiS10CibLR4lqYsCGs3oJknS0muWzRaPktRINbq6KUlSB52y2QlzJEmSJEkT2fMoSZ00uropSVILjbLZ4lGSOmk0rkKSpBYaZbPFoyQ1kkZXNyVJ6qBTNls8SlIXRatbYyRJWnrNstkJcyRJkiRJE9nzKEltpNW4CkmSll+vbLZ4lKROGt0aI0lSC42y2eJRkjppFFCSJLXQKJsd8yhJkiRJmsieR0nqpNHVTUmSWmiUzRaPktRF0WpQviRJS69ZNls8SlIjnR5ELElSB52y2TGPktRJrWOZIMk1SfYluWOV/Unyb5LsTvKlJC8Z27cjyT3DsmM9pyhJ0lJplM0Wj5KkaV0LXHyE/T8FnDMsVwK/AZDkVOCdwI8BFwDvTHLKTFsqSdKx4VrmmM0Wj5KkqVTVZ4H9RzjkEuC3a+RzwMlJTgdeA9xcVfur6iHgZo4cdJIkaQrzzmbHPEpSI5s8ruIM4L6x9T3DttW2S5LUXqdsnmvx+LQHHuWv/9pX5vmTkrSQ9n3r0dl88fpmdDstya6x9aur6uqj+PzhfryOsF0L4LHnbOVr//RHNrsZkrTpHnvf1tl8caNstudRkrqYcnD9ETxYVdvX8fk9wFlj62cCe4ftL1+x/TPr+B1JkpZDs2x2zKMkaaPsBP7hMLPbS4FvVdX9wE3ARUlOGQbjXzRskyRJs7Wh2WzPoyR1MsObQZN8iNFVytOS7GE0S9sJAFX1fuBG4LXAbuAvgTcP+/Yn+SXg1uGr3lVVRxrcL0lSH42yeWLxmGQr8FngacPxN1TVO5OcDVwPnAp8Afjpqnp8+lOVJG20WQ7Kr6rLJ+wv4C2r7LsGuGYW7ToWmc2StDw6ZfM0t60+Bryyqs4DzgcuHro83w28p6rOAR4CrjiaH5YkzcAMH0SshWI2S9KyaJTNE4vH4ZkgfzGsnjAsBbwSuGHYfh1w6UxaKEmSnsRsliRthqkmzEmyJcntwD5GD5D8GvBwVT0xHLLqc0GSXJlkV5Jdj9eMpqaXJI00urqpI9uobD7wyCPzabAkHasaZfNUxWNVHaiq8xlN4XoBcO7hDlvls1dX1faq2n5iZvTsFEkSqfUtWi4blc1btm2bZTMl6ZjWLZuParbVqno4yWeAlwInJzl+uMJ56HkhkqTNtL4HEWsJmc2StOAaZfPEnsckz0py8vD++4CfBO4GPg28YThsB/DxWTVSkjSlRrfGaHVmsyQtkUbZPE3P4+nAdUm2MCo2P1JVv5/kLuD6JL8MfBH4wAzbKUmSvsdsliTN3cTisaq+BLz4MNvvZTTGQpK0IBZxfIQ2ntksScujUzYf1ZhHSdKCaxRQkiS10CibLR4lqYsFnZlNkqRjVrNsnupRHZIkSZKkY5s9j5LUSaOrm5IktdAomy0eJamTRgElSVILjbLZ4lGSGuk0rkKSpA46ZbNjHiVJkiRJE1k8SpIkSZIm8rZVSeqk0a0xkiS10CibLR4lqYtmz5KSJGnpNctmi0dJ6qRRQEmS1EKjbLZ4lKROGgWUJEktNMpmJ8yRJEmSJE1kz6MkNRF6jauQJGnZdctmi0dJ6qRRQEmS1EKjbLZ4lKQums3oJknS0muWzY55lCRJkiRNZM+jJHXS6OqmJEktNMpmi0dJ6qRRQEmS1EKjbLZ4lKRGOo2rkCSpg07ZbPEoSZ00CihJklpolM1OmCNJkiRJmsieR0nqomh1dVOSpKXXLJstHiWpkU7jKiRJ6qBTNnvbqiR1UutYppDk4iRfTbI7yVWH2f+eJLcPy58keXhs34GxfTvXc5qSJC2NRtlsz6MkNTLLq5tJtgDvBV4N7AFuTbKzqu46dExV/fzY8T8LvHjsK75bVefProWSJC2eTtlsz6MkaVoXALur6t6qehy4HrjkCMdfDnxoLi2TJOnYNNdstniUpE7Wd2vMaUl2jS1Xrvj2M4D7xtb3DNueIskPAmcDnxrbvHX43s8luXQ9pylJ0tJolM3etipJXax/RrcHq2r7EfZnlV89nMuAG6rqwNi251bV3iTPBz6V5MtV9bW1NlaSpIXXLJvteZSkJrLOZQp7gLPG1s8E9q5y7GWsuC2mqvYOr/cCn+HJYy4kSWqnWzZbPEqSpnUrcE6Ss5OcyCiEnjIzW5IfBk4B/mhs2ylJnja8Pw24ELhr5WclSdJRmWs2e9uqJHUywxndquqJJG8FbgK2ANdU1Z1J3gXsqqpDYXU5cH1VjbfmXOA3kxxkdOHyV8dngpMkqa1G2WzxKEmNzPpBxFV1I3Djim3vWLH+vx3mc38I/K2ZNk6SpAXUKZstHiWpkxkHlCRJOkqNstniUZI6aRRQkiS10CibnTBHkiRJkjSRPY+S1EXNflyFJEk6Cs2y2eJRkjppFFCSJLXQKJstHiWpkU5XNyVJ6qBTNls8SlInjQJKkqQWGmWzE+ZIkiRJkiay51GSGul0a4wkSR10yuaJPY9Jzkry6SR3J7kzyduG7acmuTnJPcPrKbNvriRpVbXORUvDbJakJdEsm6e5bfUJ4Beq6lzgpcBbkrwAuAq4parOAW4Z1iVJm6lRQOmIzGZJWhaNsnli8VhV91fVF4b33wHuBs4ALgGuGw67Drh0Vo2UJEnfYzZLkjbDUY15TPI84MXA54HnVNX9MAqxJM9e5TNXAlcCbD3u6etpqyTpCEKvcRWaznqz+fiTvLNVkmalWzZPPdtqkqcDHwV+rqq+Pe3nqurqqtpeVdtPzNa1tFGSNK1Gt8Zoso3I5i3bts2ugZKkVtk8Vc9jkhMYhdPvVNXvDpsfSHL6cGXzdGDfrBopSZpOagGTRjNhNkvScuiUzdPMthrgA8DdVfVrY7t2AjuG9zuAj2988yRJU2s2o5tWZzZL0pJols3T9DxeCPw08OUktw/bfhH4VeAjSa4AvgG8cTZNlCRJK5jNkqS5m1g8VtV/ZjTW83BetbHNkSStR6dB+Vqd2SxJy6NTNh/VbKuSpAXXKKAkSWqhUTZbPEpSI52ubkqS1EGnbLZ4lKROGgWUJEktNMrmqZ/zKEmSJEk6dtnzKEldVK9bYyRJWnrNstniUZI6aRRQkiS10CibLR4lqYnQ6+qmJEnLrls2O+ZRkiRJkjSRPY+S1Ek1urwpSVIHjbLZ4lGSGul0a4wkSR10ymaLR0nqomg1KF+SpKXXLJstHiWpkRzc7BZIkqRxnbLZCXMkSZIkSRPZ8yhJnTS6NUaSpBYaZbM9j5LUSGrty1Tfn1yc5KtJdie56jD7/1GSP09y+7D8zNi+HUnuGZYdG3fWkiQtrk7ZbM+jJHVRzHQ68CRbgPcCrwb2ALcm2VlVd6049MNV9dYVnz0VeCewfWjpbcNnH5pZgyVJ2mzNstmeR0lqZMZXNy8AdlfVvVX1OHA9cMmUTXsNcHNV7R9C6Wbg4rWcoyRJy6RTNls8SpIOOS3JrrHlyhX7zwDuG1vfM2xb6b9L8qUkNyQ56yg/K0mSvmehstnbViWpk/XdGfNgVW0/wv5M8Yv/F/ChqnosyT8GrgNeOeVnJUnqp1E22/MoSU2Emd8aswc4a2z9TGDv+AFV9c2qemxY/XfAj077WUmSuumWzRaPktRF1fqWyW4FzklydpITgcuAneMHJDl9bPX1wN3D+5uAi5KckuQU4KJhmyRJfTXLZm9blSRNpaqeSPJWRsGyBbimqu5M8i5gV1XtBP5ZktcDTwD7gX80fHZ/kl9iFHIA76qq/XM/CUmSGpl3Nls8SlIj0z4Taq2q6kbgxhXb3jH2/u3A21f57DXANTNtoCRJC6ZTNls8SlInTkEjSdJiaZTNFo+S1Misr25KkqSj0ymbLR4lqYsCDjZKKEmSll2zbHa2VUmSJEnSRPY8SlInfS5uSpLUQ6NstniUpEY6jauQJKmDTtls8ShJnUz3QGFJkjQvjbLZ4lGSGul0dVOSpA46ZbMT5kiSJEmSJrLnUZK6KFoNypckaek1y2aLR0lqIkAajauQJGnZdctmi0dJ6uTgZjdAkiQ9SaNsdsyjJEmSJGkiex4lqZFOt8ZIktRBp2y2eJSkLpoNypckaek1y2aLR0lqo1o9iFiSpOXXK5stHiWpkU4PIpYkqYNO2eyEOZIkSZKkiex5lKROGt0aI0lSC42yeWLPY5JrkuxLcsfYtlOT3JzknuH1lNk2U5I0UUEOrn3R8jCbJWlJNMvmaW5bvRa4eMW2q4Bbquoc4JZhXZK02arWvmiZXIvZLEnLoVE2Tyweq+qzwP4Vmy8BrhveXwdcusHtkiRJqzCbJUmbYa1jHp9TVfcDVNX9SZ692oFJrgSuBNh63NPX+HOSpKks3kVKzc+asvn4k7y7VZJmqlE2z3zCnKq6Grga4KTjn9XofzpJWjxZwFtctHjGs3nrGWf5RyNJM9Qpm9f6qI4HkpwOMLzu27gmSZLWrNG4Ch01s1mSFlGjbF5r8bgT2DG83wF8fGOaI0laswIOrmPRsjObJWnRNMvmaR7V8SHgj4AfTrInyRXArwKvTnIP8OphXZIkzYHZLEnaDBPHPFbV5avsetUGt0WStA6hWo2r0OrMZklaDt2yeeYT5kiS5qhRQEmS1EKjbLZ4lKROGgWUJEktNMpmi0dJ6uLQoHxJkrQYmmXzWmdblSQN3gb7AAAKoUlEQVRJkiQdQyweJamRVK15mer7k4uTfDXJ7iRXHWb//5jkriRfSnJLkh8c23cgye3DsnMDT1uSpIXVKZu9bVWSOpnhuIokW4D3MnoMxB7g1iQ7q+quscO+CGyvqr9M8k+A/wP4H4Z9362q82fWQEmSFlGjbLbnUZLaqFFArXWZ7AJgd1XdW1WPA9cDlzypBVWfrqq/HFY/B5y5oacoSdJS6ZXNFo+SpENOS7JrbLlyxf4zgPvG1vcM21ZzBfAHY+tbh+/9XJJLN6jNkiR1tlDZ7G2rktRFsd5bYx6squ1H2J9VfvWpByb/ANgOvGxs83Oram+S5wOfSvLlqvra2psrSdKCa5bNFo+S1MlspwPfA5w1tn4msHflQUl+EvhfgJdV1WOHtlfV3uH13iSfAV4MWDxKknprlM3etipJjcx4RrdbgXOSnJ3kROAy4EkzsyV5MfCbwOurat/Y9lOSPG14fxpwITA+mF+SpJY6ZbM9j5LUyQxndKuqJ5K8FbgJ2AJcU1V3JnkXsKuqdgL/Cng68B+TAHyjql4PnAv8ZpKDjC5c/uqKmeAkSeqpUTZbPEqSplZVNwI3rtj2jrH3P7nK5/4Q+FuzbZ0kSceeeWazxaMkdVHAwdld3ZQkSUepWTZbPEpSG1M/E0qSJM1Fr2y2eJSkThoFlCRJLTTKZotHSeqkUUBJktRCo2z2UR2SJEmSpInseZSkLpoNypckaek1y2aLR0lqo6AObnYjJEnSX+mVzRaPktRJo3EVkiS10CibHfMoSZIkSZrInkdJ6qLZuApJkpZes2y2eJSkThrdGiNJUguNstniUZI6aRRQkiS10CibLR4lqY1qFVCSJC2/XtnshDmSJEmSpInseZSkLgo42OdZUpIkLb1m2WzxKEmdNLo1RpKkFhpls8WjJHXSKKAkSWqhUTZbPEpSG9XqWVKSJC2/XtnshDmSJEmSpInseZSkLgqq+gzKlyRp6TXLZotHSeqk0a0xkiS10CibLR4lqZNGg/IlSWqhUTY75lGSJEmSNJE9j5LURVWrBxFLkrT0mmWzxaMkddLo1hhJklpolM0Wj5LUSDW6uilJUgedstniUZLaqFZXNyVJWn69stkJcyRJkiRJE9nzKEldFK2eJSVJ0tJrls3r6nlMcnGSrybZneSqjWqUJGmN6uDaF7VgNkvSgmmUzWvueUyyBXgv8GpgD3Brkp1VdddGNU6SNL0CqtHVTR09s1mSFku3bF5Pz+MFwO6qureqHgeuBy7ZmGZJko5a1cyvbk7q1UrytCQfHvZ/Psnzxva9fdj+1SSv2bDz1jizWZIWSbNsXk/xeAZw39j6nmHbysZemWRXkl2P16Pr+DlJ0mYa69X6KeAFwOVJXrDisCuAh6rqbwDvAd49fPYFwGXAC4GLgfcN36eNddTZfOCRR+bWOEnSxpp3Nq+neMxhtj2lT7aqrq6q7VW1/cRsXcfPSZImqYO15mUK0/RqXQJcN7y/AXhVkgzbr6+qx6rqT4Hdw/dpYx11Nm/Ztm0OzZKkY1enbF5P8bgHOGts/Uxg7zq+T5K0XrO9NWaaXq2/OqaqngC+BTxzys9q/cxmSVo0jbJ5PY/quBU4J8nZwJ8x6vL8+0f6wLcPPPjgTfv/3f8LnAY8uI7fXiSey2Lqci5dzgM8l5V+cCMaMu47PHTT/103nLaOr9iaZNfY+tVVdfXY+jS9WqsdM1WPmNbtqLP5sb17Htz9v/6C2by4PJfF0+U8wHNZyWyekM1rLh6r6okkbwVuArYA11TVnRM+8yyAJLuqavtaf3uReC6Lqcu5dDkP8FzmoaounvFPTNOrdeiYPUmOB04C9k/5Wa2T2TziuSymLufS5TzAc5mHbtm8ruc8VtWNVfVDVfXXq+pX1vNdkqSF91e9WklOZNSrtXPFMTuBHcP7NwCfqqoatl82zPh2NnAO8MdzavcxxWyWpGPKXLN5PbetSpKOIav1aiV5F7CrqnYCHwD+Q5LdjK5qXjZ89s4kHwHuAp4A3lJVBzblRCRJamLe2bxZxePVkw9ZGp7LYupyLl3OAzyXFqrqRuDGFdveMfb+UeCNq3z2VwB7whZXp79rz2UxdTmXLucBnksL88zmjHosJUmSJEla3brGPEqSJEmSjg1zLx6TXJzkq0l2J7lq3r+/HkmuSbIvyR1j205NcnOSe4bXUzazjdNIclaSTye5O8mdSd42bF/Gc9ma5I+T/NfhXP73YfvZST4/nMuHhwHECy/JliRfTPL7w/qynsfXk3w5ye2Hppdexr8vgCQnJ7khyVeGf2Z+fFnPRVqN2bz5zObFZTYvHrN588y1eEyyBXgv8FPAC4DLk7xgnm1Yp2uBldPtXgXcUlXnALcM64vuCeAXqupc4KXAW4b/H5bxXB4DXllV5wHnAxcneSnwbuA9w7k8BFyxiW08Gm8D7h5bX9bzAHhFVZ0/Nm32Mv59Afw68Imq+hHgPEb//yzruUhPYTYvDLN5cZnNi8ds3iTz7nm8ANhdVfdW1ePA9cAlc27DmlXVZxnNUDTuEuC64f11wKVzbdQaVNX9VfWF4f13GP0DdwbLeS5VVX8xrJ4wLAW8Erhh2L4U55LkTODvAL81rIclPI8jWLq/ryTfD/wEo1nKqKrHq+phlvBcpCMwmxeA2byYzObFYzZvrnkXj2cA942t7xm2LbPnVNX9MPoXP/DsTW7PUUnyPODFwOdZ0nMZbie5HdgH3Ax8DXi4qp4YDlmWv7N/Dfxz4OCw/kyW8zxg9B8Jn0xyW5Irh23L+Pf1fODPgX8/3LL0W0m2sZznIq3GbF4wZvNCMZsXj9m8ieZdPOYw25zudZMkeTrwUeDnqurbm92etaqqA1V1PnAmoyvo5x7usPm26ugkeR2wr6puG998mEMX+jzGXFhVL2F0G9xbkvzEZjdojY4HXgL8RlW9GHgEb4NRP8v875p2zObFYTYvLLN5E827eNwDnDW2fiawd85t2GgPJDkdYHjdt8ntmUqSExiF0+9U1e8Om5fyXA4Zbln4DKOxIicnOfQc02X4O7sQeH2SrzO6ZeyVjK52Ltt5AFBVe4fXfcDHGP2HwzL+fe0B9lTV54f1GxgF1jKei7Qas3lBmM0Lx2xeTGbzJpp38XgrcM4wS9WJwGXAzjm3YaPtBHYM73cAH9/EtkxluF//A8DdVfVrY7uW8VyeleTk4f33AT/JaJzIp4E3DIct/LlU1dur6syqeh6jfy4+VVVvYsnOAyDJtiTPOPQeuAi4gyX8+6qq/wbcl+SHh02vAu5iCc9FOgKzeQGYzYvHbF5MZvPmStV8e9qTvJbRVZstwDVV9StzbcA6JPkQ8HLgNOAB4J3A7wEfAZ4LfAN4Y1WtHLi/UJL8beA/AV/me/fw/yKjsRXLdi4vYjQoegujiyEfqap3JXk+o6uEpwJfBP5BVT22eS2dXpKXA/9TVb1uGc9jaPPHhtXjgQ9W1a8keSZL9vcFkOR8RhMlnAjcC7yZ4W+NJTsXaTVm8+Yzmxeb2bxYzObNM/fiUZIkSZK0fOZ926okSZIkaQlZPEqSJEmSJrJ4lCRJkiRNZPEoSZIkSZrI4lGSJEmSNJHFoyRJkiRpIotHSZIkSdJEFo+SJEmSpIn+f7Vr/lffs2kTAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x7f1214072438>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "init()\n",
+    "plot_ρs()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Check the first time step against reference data\n",
+    "\n",
+    "The reference data was obtained with the [sample code](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp) after making the following changes:\n",
+    "```c++\n",
+    "const int nsteps = 1000;\n",
+    "const int noutput = 1;\n",
+    "const int nfluids = 2;\n",
+    "const double gA = 0;\n",
+    "```\n",
+    "\n",
+    "Remove the next cell if you changed the parameters at the beginning of this notebook."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "init()\n",
+    "time_loop(1)\n",
+    "ref_a = np.array([0.133183, 0.0921801, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.0921801, 0.133183, 0.719568, 1.05507, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 1.05507, 0.719568])\n",
+    "ref_b = np.array([0.719568, 1.05507, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 1.05507, 0.719568, 0.133183, 0.0921801, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.0921801, 0.133183])\n",
+    "assert np.allclose(dh.gather_array(ρ_a.name)[0], ref_a)\n",
+    "assert np.allclose(dh.gather_array(ρ_b.name)[0], ref_b)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Run the simulation until converged"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA48AAAF0CAYAAACDowz8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XuwZWV55/Hv0zfQ9gLYghSNgjPEoAkX7UJTpBRFEY0DTI1mIMYQhxSVlGbMPZipkRlMqjCZijEVo/ZoB0wpaFCkJ4UigzLEqIRGCXIRaQmRThNabBAF6bbPeeaPvTpuD+ecvc/Ze+3Lc76fqlVnr8ve+116OL9+1rved0VmIkmSJEnSYlaNuwGSJEmSpMln8ShJkiRJ6sniUZIkSZLUk8WjJEmSJKkni0dJkiRJUk8Wj5IkSZKkniweJUmSJEk9WTxKkiRJknqyeFR5EfHUiNgcEQ9FxK6I+M1xt0mSpJXMbJamk8WjVoJPAd8EngWcDfyviHjWeJskSdKKZjZLU8jiUaVFxOsAMvNdmbknMz8H/AvwE83+F0TETERsHGc7JUlaKRbL5oj4m4j4+2b5XEQcON7WSupm8ajqzgCu2r8SEauApwMPNJt+H/hr4NjRN02SpBVpsWz+CeBlmXky8Ahw3FhaKGleFo+q7sXAd7rWXwE8mJl3RcRxwP3ANVg8SpI0KvNmM/BPwKrM3BcRTwEOB74xhvZJWoDFo8qKiLXAMcDrI+LAiHgB8Jd0ehsBfhN4F3AHFo+SJLWuRzYfCzwrIq4H7gbem5kPj62xkp5gzbgbILXoWOBe4DY6t8LsAv4wM6+IiBOAk4GPAqubRZIktWuxbH4j8L7MfEdEPB34MvDhsbVU0hNYPKqy44A7M/O/A/99zr4LgBdn5kMAEfEPo26cJEkr0GLZ/FPATc3rg4HvjrJhknqzeFRlxwN3zt0YES8CfrC/cGw8HhHPyMzvzD1ekiQNzbzZ3HgBcFpEvAWYBX51ZK2S1BfHPKqy44Cvz92YmTdn5pvnbHuphaO0uIg4MiI+HxF3RsTtEfG2eY6JiPjziNgeEbdGxAu79p0bEXc3y7mjbb2kCTFvNgNk5hmZ+aLMPDUzX5WZt4y4bdLUGXU2R2YO+xwkSQVFxOHA4Zn5lYh4KnAzcFZm3tF1zGuBXwdeS2dGxfdk5osj4hBgG7AJyOa9L5pzB4AkSVqCUWezPY+SpL5k5v2Z+ZXm9ffo3Hp2xJzDzgQ+nB1fBg5qgu3VwLWZubsJpWuB00fYfEmSyhl1Nls8SpKWLCKOAk4Ebpyz6wjgvq71Hc22hbZLkqQhGEU2j3TCnLUHrM8D1h8yyq+UpIm059Hd/HDPozHMz3z1y9fnd3bPLPv9N9+653bg8a5NmzNz89zjmod3fwL4jcx8ZO7ueT46F9muMVu3+sn5pLVPG3czJGnsfvDDR9g785jZvIiRFo8HrD+E4059whhOSVpxbr3uPUP/zAd3z3DjNRuX/f61h3/z8czctNgxzQO+PwF8JDM/Oc8hO4Aju9Y3Ajub7afM2X79shuroXnS2qfxM89x/iJJ+tI/Xzr0z6yWzd62KknqS0QE8CE6z2j70wUO2wr8UjOz20uA72bm/cA1dKbgPzgiDgZOa7ZJkqRlGnU2+5xHSSojmcnZNr/gZOBNwNciYv8U+n8APBsgM98PXE1nNrftwGPAm5t9uyPinfzoAeAXZebuNhsrSdL41crmkRaPMwfAQ89bPcqvlKSJNPOF4X9mArMtDiPMzC8w//iI7mMSeMsC+7YAW1pomgaxbx888OC4WyFJ47dv39A/slo22/MoSYXM0urVTUmStESVstkxj5IkSZKknux5lKQikmQmffqFJEmTolo2WzxKUiFtjquQJElLVymbR1o8zq6BPYfUuedXkpZrtoW/vgnMFAoojchsknv3jrsVkjR+s8PP0GrZbM+jJBVS6eqmJEkVVMpmJ8yRJEmSJPVkz6MkFZFQalC+JEnTrlo2WzxKUiGOKpckabJUyua+iseIOAj4IPBTdAro/wLcBXwMOAq4F/j5zHxo8Q+CtFyVJIjhf2SSpQbla3FDy2ZJUmuqZXO/Yx7fA3wmM38SOB64E7gAuC4zjwGua9YlSeOSMDPAoqljNkvSpCuWzT2Lx4h4GvBS4EMAmbk3Mx8GzgQubQ67FDirrUZKkqQfMZslSePQT8/jc4FvA38VEV+NiA9GxHrgsMy8H6D5eeh8b46I8yNiW0Rsm/n+o0NruCTpxyWdcRXLXTRVhpbNe9kzulZL0gpTLZv7KR7XAC8E3peZJwKPsoTbYDJzc2ZuysxNq5+yfpnNlCT1FswMsGiqDC2b13FAW22UJBXL5n6Kxx3Ajsy8sVm/gk5gPRARhwM0P3e100RJUj8SmM3lL5oqZrMkTYFq2dyzeMzMfwXui4jnNZtOBe4AtgLnNtvOBa5qpYWSJOnHmM2SpHHo98EZvw58JCLWAfcAb6ZTeH48Is4DvgW8oZ0mSpL6NYm3uKg1ZrMkTYFK2dxX8ZiZtwCb5tl16nCbI0larqRWQGlxZrMkTb5q2dxvz6MkaQrMZp2AkiSpgkrZbPEoSUVUu7opSdK0q5bNoy8eJ3DWIEmSJEnS4ux5lKQikmCmrycwSZKkUaiWzRaPklRIpXEVkiRVUCmbLR4lqYhq4yokSZp21bLZ4lGSyghmss6tMZIkTb9a2Tza4jFxwhxJAv8WSpKkqWPPoyQVkcBsoUH5kiRNu2rZbPEoSYVUGlchSVIFlbLZ4lGSisisNa5CkqRpVy2b65yJJEmSJKk1I+95jNlRf6MkrRyzhW6N0YhEwCqvJUsS0U6GVspmb1uVpCI6z5KyCJAkaVJUy2aLR0kqo91xFRGxBXgdsCszf2qe/b8LvLFZXQMcCzwzM3dHxL3A94AZYF9mbmqtoZIkTYxa2VynDJakFW7/dODLXfpwCXD6gt+f+SeZeUJmngC8Hfh/mbm765CXN/stHCVJK0K1bLZ4lCT1JTNvAHb3PLDjHOCyFpsjSdKKN+psHv1tq1lnwKgkTZqZwf7GboiIbV3rmzNz81I/JCKeTOcq6Fu7Nifw2YhI4APL+VxJkqZRpWx2zKMkFZHEoIPyHxzSLaX/Afj7ObfFnJyZOyPiUODaiPh6c7VUkqSyqmWzxaMkFTI7GQ8iPps5t8Vk5s7m566IuBI4CbB4lCSVVymbJ+JMJEmD2z8d+HKXYYiIpwMvA67q2rY+Ip66/zVwGnDbUL5QkqQJVi2b7XmUJPUlIi4DTqEz/mIHcCGwFiAz398c9h+Bz2bmo11vPQy4MjoPX14DfDQzPzOqdkuSVNWos3m0xWNCzI70GyVpMmUbHxmDDspf/PMzz+njmEvoTBveve0e4Ph2WiVJ0uSqls32PEpSIX0+E0qSJI1IpWy2eJSkIjJhZjIG5UuSJOplc50zkSRJkiS1xp5HSSojmKW9cRWSJGmpamXz6IvHFiaJkCQ104EXujVGIxR1/mEjSZOkWjbb8yhJhQzrmVCSJGk4KmWzxaMkFZEEsy1OBy5JkpamWjbXKYMlSZIkSa2x51GSCql0a4wkSRVUyuaRFo/RLJK00rXxtzCB2UKD8iVJmnbVstmeR0kqI5jxEp0kSROkVjZbPEpSEdWubkqSNO2qZXOdM5EkSZIktcaeR0kqpNKtMZIkVVApm0dfPObIv1GSVoTMKHVrjCRJ065aNtvzKEmFzBQKKEmSKqiUzXXORJIkSZLUGnseJamIBGYLjauQJGnaVctmi0dJKiNK3RojSdL0q5XNfRWPEXEv8D1gBtiXmZsi4hDgY8BRwL3Az2fmQz0/zAlzJKkVnWdJ1bm6qcUNNZslSa2ols1LKYNfnpknZOamZv0C4LrMPAa4rlmXJI3RDKuWvWgqmc2SNOEqZfMgLToTuLR5fSlw1uDNkSRJAzCbJUmt6bd4TOCzEXFzRJzfbDssM+8HaH4eOt8bI+L8iNgWEdtmHn108BZLkuaVBLO5/EVTZyjZvDcfH1FzJWnlqZbN/U6Yc3Jm7oyIQ4FrI+Lr/X5BZm4GNgMcuPFIRzxKUotmJ/AWF7VmKNn89NUbzGZJalGlbO6reMzMnc3PXRFxJXAS8EBEHJ6Z90fE4cCuFtspSeohE2Ym8Cql2mE2S9Lkq5bNPYvHiFgPrMrM7zWvTwMuArYC5wIXNz+varOhkqTeJvEWFw3f0LN5VZ2r4pI0aSplcz89j4cBV0bE/uM/mpmfiYibgI9HxHnAt4A3tNdMSZLUxWyWJI1cz+IxM+8Bjp9n+3eAU9tolCRp6TqD8u1BWgnMZkmaDtWyud8JcyRJU2CGOrfGSJJUQaVstniUpCKSWuMqJEmadtWyeeTFYzghuCRJkiRNHXseJamMWuMqJEmafrWyuc6ZSJKYJZa99BIRWyJiV0TctsD+UyLiuxFxS7O8o2vf6RFxV0Rsj4gLhnjKkiRNtErZbM+jJBUxggcRXwL8BfDhRY75u8x8XfeGiFgNvBd4FbADuCkitmbmHW01VJKkSVAtmy0eJamQNm+NycwbIuKoZbz1JGB783gJIuJy4EzA4lGSVF6lbPa2VUnSfhsiYlvXcv4yPuNnIuIfI+LTEfGCZtsRwH1dx+xotkmSpMVNVDbb8yhJRXQeRDzQrTEPZuamAd7/FeA5mfn9iHgt8CngGJh30IZzb0uSyquWzfY8SlIhbQ7K7yUzH8nM7zevrwbWRsQGOlczj+w6dCOwc+AvlCRpClTKZnseJamIcT+IOCKeBTyQmRkRJ9G5QPkd4GHgmIg4GvgX4GzgF8bWUEmSRqRaNls8SpL6EhGXAafQGX+xA7gQWAuQme8HXg/8WkTsA34AnJ2ZCeyLiLcC1wCrgS2ZefsYTkGSpFJGnc0jLx7HWHhLUnktz+h2To/9f0FnuvD59l0NXN1GuyRJmmSVstmeR0mqIgcelC9JkoapWDZbPEpSEQlDGVwvSZKGo1o2WzxKUiGVrm5KklRBpWz2UR2SJEmSpJ7seZSkIsY9Hbim2OzsuFsgSSVVy2aLR0kqpFJASZJUQaVstniUpCKSWjO6SZI07apls8WjJBVSaUY3SZIqqJTNTpgjSZIkSepp9D2PdQpvSZosWWtchSRJU69YNnvbqiQVUW1GN0mSpl21bLZ4lKRCKgWUJEkVVMpmxzxKkiRJknqy51GSiqg2HbgkSdOuWjY7YY4kFZKFAkqSpAoqZbM9j5JUSKVnSUmSVEGlbLZ4lKQisth04JIkTbtq2eyEOZIkSZKknux5lKRCKo2rkCSpgkrZPNLiMZtFkla6dv4W1prRTZKk6Vcrm+15lKRCKl3dlCSpgkrZbPEoSUUktQblS5I07aplsxPmSJIkSZJ6sudRkqrIzpTgkiRpQhTL5tEXj3V6bSVp4lR6ELFGqNK/bCRpwlTKZnseJamIpNagfEmSpl21bHbMoyRJkiSpJ3seJamMWs+SkiRp+tXKZotHSSrEoWuSJE2WStnc922rEbE6Ir4aEX/brB8dETdGxN0R8bGIWNdeMyVJ/ciMZS+aPmazJE2+Stm8lJ7HtwF3Ak9r1t8FvDszL4+I9wPnAe9b9BMC0lGWktTKzNOZtQblqy+DZ7MkqTXVsrmvUi4iNgI/B3ywWQ/gFcAVzSGXAme10UBJkvREZrMkadT67Qf8M+D3gNlm/RnAw5m5r1nfARwx3xsj4vyI2BYR22YefXSgxkqSFjebseyll4jYEhG7IuK2Bfa/MSJubZYvRsTxXfvujYivRcQtEbFtiKe8kg0lm/fm4+23VJJWsErZ3LN4jIjXAbsy8+buzfMcOu9Q0MzcnJmbMnPT6vXr+2mTJGmZOrfHLG/pwyXA6Yvs/yfgZZl5HPBOYPOc/S/PzBMyc9Nyzk0/MsxsXhcHttJGSVJHpWzuZ8zjycAZEfFa4EA64yr+DDgoItY0Vzg3Ajv7+UJJUnvaHFeRmTdExFGL7P9i1+qX6WSD2mE2S9KUqJTNPYvHzHw78HaAiDgF+J3MfGNE/A3weuBy4Fzgqr6+MQrNVStJEyQZeGa2DXNuW9mcmXOvUPbrPODTXesJfDYiEvjAAJ8rWshmSVIrqmXzIM95/H3g8oj4Q+CrwIcG+CxJ0vg9OIxbSiPi5XQC6me7Np+cmTsj4lDg2oj4embeMOh36QnMZkmqZaKyeUnFY2ZeD1zfvL4HOGkp75cktWvc93ZExHF0Zv98TWZ+Z//2zNzZ/NwVEVfSyQ+LxyEwmyVpslXKZp+6KElV5HgfRBwRzwY+CbwpM7/RtX19RDx1/2vgNGDeWeEkSSqlWDYPctuqJGnStHh5MyIuA06hM/5iB3AhsBYgM98PvIPO4yL+svPIQfY1t9ocBlzZbFsDfDQzP9NeSyVJmiCFsnnkxWPa1ylJUykzz+mx/1eAX5ln+z3A8U98hyZCJszO9j5Okqrr89kYk2TU2WzPoyQV0uZ04JIkaekqZbPFoyQVMoUXTSVJKq1SNls8SlIRSa2rm5IkTbtq2WzxKElVJFAooCRJmnrFsnm0xWM0iyStdP4tlCRJU8aeR0kqpNK4CkmSKqiUzRaPklRJoYCSJKmEQtls8ShJZUSpQfmSJE2/Wtls8ShJlRS6uilJUgmFsnn0xWOdwluSJEmSVgx7HiWpiqz1LClJkqZesWy2eJSkSgrdGiNJUgmFstniUZJKqXN1U5KkGupk86pxN0CSJEmSNPnseZSkSgrdGiNJUgmFstniUZIqKRRQkiSVUCibLR4lqYoECs3oJknS1CuWzRaPklRIFrq6KUlSBZWy2QlzJEmSJEk9jbbnMSH2jfQbJWkytXUVstDVTUmSSiiUzd62KkmVFBpXIUlSCYWy2eJRkgqJQlc3JUmqoFI2WzxKUhVJqVtjJEmaesWy2QlzJEmSJEk9jbTncdU+OGC39aokrWpl8rAoNa5CI7IqiHXrxt0KSRq/vW1kaK1s9rZVSaqk0K0xkiSVUCibLR4lqZJCASVJUgmFstl7SCVJkiRJPdnzKEmVFLq6KUlSCYWyeaTF4+o9cPBdM6P8SkmaSPftaeFDk1KD8jUia9bAYRvG3QpJGr/HWyiNimWzPY+SVEilBxFLklRBpWx2zKMkVZIDLD1ExJaI2BURty2wPyLizyNie0TcGhEv7Np3bkTc3SznDnKKkiRNlULZbPEoSerXJcDpi+x/DXBMs5wPvA8gIg4BLgReDJwEXBgRB7faUkmSVoZLGGE2WzxKkvqSmTcAuxc55Ezgw9nxZeCgiDgceDVwbWbuzsyHgGtZPOgkSVIfRp3NjnmUpELGPK7iCOC+rvUdzbaFtkuSVF6lbB7tbKuP/ZCn3frtUX6lJE2k1Y/9sJ0PHmxGtw0Rsa1rfXNmbl7C++f78lxkuybAzJPX8shxzxx3MyRp7Ga+vbadDy6UzfY8SlIVfQ6uX8SDmblpgPfvAI7sWt8I7Gy2nzJn+/UDfI8kSdOhWDY75lGSNCxbgV9qZnZ7CfDdzLwfuAY4LSIObgbjn9ZskyRJ7RpqNtvzKEmVtHgzaERcRucq5YaI2EFnlra1AJn5fuBq4LXAduAx4M3Nvt0R8U7gpuajLsrMxQb3S5JUR6Fs7lk8RsSBwA3AAc3xV2TmhRFxNHA5cAjwFeBNmbm3/1OVJA1bm4PyM/OcHvsTeMsC+7YAW9po10pkNkvS9KiUzf30PO4BXpGZ34+ItcAXIuLTwG8B787MyyPi/cB5NM8NWdC+ffDAg0tpnyTVtG9fO5/rNDQrxdCyeeYAeOh5q9tvsSRNuJkvtPTBhbK555jH5pkg329W1zZLAq8Armi2Xwqc1UoLJUnSjzGbJUnj0NeEORGxOiJuAXbReYDkN4GHM3P/pfMFnwsSEedHxLaI2LZ39vFhtFmStJAcYNFUGVY273vs0dE0WJJWqkLZ3FfxmJkzmXkCnSlcTwKOne+wBd67OTM3ZeamdasOXH5LJUmLihxs0XQZVjavefL6NpspSStatWxe0myrmflwRFwPvAQ4KCLWNFc49z8vRJI0ToM9iFhTyGyWpAlXKJv7mW31mcAPm3B6EvBK4F3A54HX05nV7Vzgqp7fNpvkXid9kyRmW7qcOIFXKTV8w8zm2TWw55DZNpsrSVNhtq2HGBbK5n7+JzocuDQiVtO5zfXjmfm3EXEHcHlE/CHwVeBDLbZTkiT9iNksSRq5nsVjZt4KnDjP9nvojLGQJE2ISRwfoeEzmyVpelTK5rY6ZyVJ41AooCRJKqFQNls8SlIVEzozmyRJK1axbLZ4lCRpJQtI/zUgSVBnUtTWGBeSVEmhq5uSJJVQKJstHiWpkkIBJUlSCYWy2eJRkgqpNK5CkqQKKmXzqnE3QJIkSZI0+SweJUmSJEk9eduqJFVS6NYYSZJKKJTNFo+SVEWxZ0lJkjT1imWzxaMkVVIooCRJKqFQNls8SlIlhQJKkqQSCmWzxaMkSStdoX/YSJLaY/EoSUUEtcZVSJI07apls8WjJFVSKKAkSSqhUDZbPEpSFcVmdJMkaeoVy+ZV426AJEmSJGny2fMoSZUUurqpEUn8vZEkaO9vYaG/sRaPklRJoYCSJKmEQtls8ShJhVQaVyFJUgWVstniUZIqKRRQkiSVUCibnTBHkiRJktTTaHseI2CV9aokETH8z3TiEy1TzI67BZJUVLFs9rZVSSqk0rgKSZIqqJTNdgNKUiU5wNKHiDg9Iu6KiO0RccE8+98dEbc0yzci4uGufTNd+7YOcpqSJE2NQtlsz6MkFdLm1c2IWA28F3gVsAO4KSK2ZuYd+4/JzN/sOv7XgRO7PuIHmXlCey2UJGnyVMpmex4lSf06Cdiemfdk5l7gcuDMRY4/B7hsJC2TJGllGmk22/MoSZUMdnVzQ0Rs61rfnJmbu9aPAO7rWt8BvHi+D4qI5wBHA5/r2nxg8/n7gIsz81MDtVbDky1M4CRJ6iiUzRaPklTF4DO6PZiZmxbZP1+FsdA3ng1ckZkzXduenZk7I+K5wOci4muZ+c3lNlaSpIlXLJu9bVWSiogBlz7sAI7sWt8I7Fzg2LOZc1tMZu5sft4DXM+Pj7mQJKmcatls8ShJ6tdNwDERcXRErKMTQk+YmS0ingccDHypa9vBEXFA83oDcDJwx9z3SpKkJRlpNnvbqiRV0uKMbpm5LyLeClwDrAa2ZObtEXERsC0z94fVOcDlmdndmmOBD0TELJ0Llxd3zwQnSVJZhbLZ4lGSCmn7QcSZeTVw9Zxt75iz/j/med8XgZ9utXFanoSYHXcjJGkCtJShlbLZ4lGSKmk5oCRJ0hIVymaLR0mqpFBASZJUQqFsdsIcSZIkSVJP9jxKUhXZ/rgKSZK0BMWy2eJRkiopFFCSJJVQKJtHXzxGn4+7lCQtWaWrmxohf28kqTWVstmeR0mqpFBASZJUQqFsdsIcSZIkSVJP9jxKUiGVbo2RJKmCStncs+cxIo6MiM9HxJ0RcXtEvK3ZfkhEXBsRdzc/D26/uZKkBeWAi6aG2SxJU6JYNvfT87gP+O3M/EpEPBW4OSKuBX4ZuC4zL46IC4ALgN9vr6mSpJ4mMGjUiqFlczSLJK10rf0tLJTNPXseM/P+zPxK8/p7wJ3AEcCZwKXNYZcCZ7XVSEmS9CNmsyRpHJY05jEijgJOBG4EDsvM+6ETYhFx6ALvOR84H+DAWD9IWyVJiwhqjatQfwbN5jUHeWerJLWlWjb3PdtqRDwF+ATwG5n5SL/vy8zNmbkpMzetiwOX00ZJUr8KjatQb8PI5tXrvbArSa0qlM199TxGxFo64fSRzPxks/mBiDi8ubJ5OLCrrUZKkvoTOYFJo1aYzZI0HSplc8/iMSIC+BBwZ2b+adeurcC5wMXNz6taaaEkqT8TepVSwzf0bPb3RpLaUSyb++l5PBl4E/C1iLil2fYHdILp4xFxHvAt4A3tNFGSJM1hNkuSRq5n8ZiZX2DhmWtPHW5zJEmDqDQoXwszmyVpelTK5iXNtipJmnCFAkqSpBIKZbPFoyQVUunqpiRJFVTKZotHSaqkUEBphPy9kaT2FPob2/dzHiVJkiRJK5c9j5JURda6NUaSpKlXLJstHiWpkkIBJUlSCYWy2eJRkooIal3dlCRp2lXL5tEXj6scZilJkiRJ08aeR0mqJAtd3pQkqYJC2WzxKEmFVLo1RpKkCipls8WjJFWRlBqUL0nS1CuWzRaPklRIzI67BZIkqVulbLZ4lCRphat0S5UkqT0Wj5JUiUWAJEmTpVA2+9wMSSokcvlLX58fcXpE3BUR2yPignn2/3JEfDsibmmWX+nad25E3N0s5w7vrCVJmlyVstmeR0mqIml1OvCIWA28F3gVsAO4KSK2ZuYdcw79WGa+dc57DwEuBDY1Lb25ee9DrTVYkqRxK5bN9jxKUiEtX908Cdiemfdk5l7gcuDMPpv2auDazNzdhNK1wOnLOUdJkqZJpWy2eJQk7bchIrZ1LefP2X8EcF/X+o5m21z/KSJujYgrIuLIJb5XkiT9yERls7etSlIlg90Z82Bmblpkf/Txjf8HuCwz90TErwKXAq/o872SJNVTKJvteZSkIoLWb43ZARzZtb4R2Nl9QGZ+JzP3NKv/G3hRv++VJKmaatls8ShJVWQOtvR2E3BMRBwdEeuAs4Gt3QdExOFdq2cAdzavrwFOi4iDI+Jg4LRmmyRJdRXLZm9blST1JTP3RcRb6QTLamBLZt4eERcB2zJzK/BfI+IMYB+wG/jl5r27I+KddEIO4KLM3D3yk5AkqZBRZ7PFoyQV0u8zoZYrM68Grp6z7R1dr98OvH2B924BtrTaQC1LzjfqRZI0FJWy2eJRkipxChpJkiZLoWy2eJSkQtq+uilJkpamUjZbPEpSFQnMFkooSZKmXbFsdrZVSZIkSVJPo+95nJ0d+VdK0opR5+KmJEk1FMpmb1uVpEIqjauQJKmCStls8ShJlfT3QGFJkjQqhbLZ4lGSCql0dVOSpAoqZbMT5kiSJEmSerLnUZKqSEoNypckaeoVy2aLR0kqIoAoNK5CIxTjboAk1VQtmy0eJakSn4YkSdJkKZTNjnmUJEmSJPVkz6MkFVLp1hhJkiqolM0Wj5JURbE6xIOdAAAM2ElEQVRB+ZIkTb1i2WzxKEllZKkHEWuEnDBHklpSK5stHiWpkEoPIpYkqYJK2eyEOZIkSZKknux5lKRKCt0aI0lSCYWyuWfPY0RsiYhdEXFb17ZDIuLaiLi7+Xlwu82UJPWUELPLXzQ9zGZJmhLFsrmf21YvAU6fs+0C4LrMPAa4rlmXJI1b5vIXTZNLGFI2p4uLi4vLvy2tKJTNPYvHzLwB2D1n85nApc3rS4GzhtwuSZK0ALNZkjQOyx3zeFhm3g+QmfdHxKELHRgR5wPnAxwY65f5dZKkvkzeRUqNzrKyec1B3t0qSa0qlM2tT5iTmZuBzQBPX72h0P90kjR5YgJvcdHk6c7mAzYe6S+NJLWoUjYv91EdD0TE4QDNz13Da5IkadkKjavQkpnNkjSJCmXzcnsetwLnAhc3P6/q+50T+D+CJJWQwATOzKaRWX42R0stkqSVrlg29/OojsuALwHPi4gdEXEenWB6VUTcDbyqWZckSSNgNkuSxqFnz2NmnrPArlOH3BZJ0gCCLDWuQgszmyVpOlTL5tYnzJEkjVChgJIkqYRC2WzxKEmVFAooSZJKKJTNFo+SVEWxQfkakYBc7tzrklRJG5OHFctm40KSJEmS1JPFoyQVEpnLXvr6/IjTI+KuiNgeERfMs/+3IuKOiLg1Iq6LiOd07ZuJiFuaZesQT1uSpIlVKZu9bVWSKmlxXEVErAbeS+cxEDuAmyJia2be0XXYV4FNmflYRPwa8MfAf272/SAzT2itgZIkTaJC2WzPoySVkZ2AWu7S20nA9sy8JzP3ApcDZ/5YCzI/n5mPNatfBjYO9RQlSZoqtbLZnkdJ0n4bImJb1/rmzNzctX4EcF/X+g7gxYt83nnAp7vWD2w+fx9wcWZ+atAGa0iizkyAklTMRGWzxaMkVZEMemvMg5m5aZH9881DN+8XRsQvApuAl3VtfnZm7oyI5wKfi4ivZeY3l99cSZImXLFstniUpEranQ58B3Bk1/pGYOfcgyLilcB/A16WmXv2b8/Mnc3PeyLieuBEwOJRklRboWx2zKMkFdLyjG43AcdExNERsQ44G/ixmdki4kTgA8AZmbmra/vBEXFA83oDcDLQPZhfkqSSKmWzPY+SVEmLM7pl5r6IeCtwDbAa2JKZt0fERcC2zNwK/AnwFOBvIgLgW5l5BnAs8IGImKVz4fLiOTPBSZJUU6FsHm3xmAmz7fbbStJUaDFI2pSZVwNXz9n2jq7Xr1zgfV8Efrrd1mm50vuQJGlqjTKb7XmUpCoSmJ3OolSSpJKKZbPFoySV0fczoSRJ0kjUymaLR0mqpFBASZJUQqFstniUpEoKBZQkSSUUymaLR0mSVrJg/kdMS9JK49/CniweJamKYoPyJUmaesWy2eJRkspISB+HJEnS5KiVzRaPklRJoXEVkiSVUCibfSywJEmSJKknex4lqYpi4yo0Qk4SIUntKJbNFo+SVEmhW2MkSSqhUDZbPEpSJYUCSpKkEgpls8WjJJWRpQJKkqTpVyubnTBHkiRJktSTPY+SVEUCs3WeJSVJ0tQrls0Wj5JUSaFbYyRJKqFQNls8SlIlhQJKkqQSCmWzxaMklZGlniUlSdL0q5XNTpgjSZIkSerJnkdJqiIhs86gfEmSpl6xbLZ4lKRKCt0aoxFJiH3jboQkTYC2IrRQNls8SlIlhQblS5JUQqFsdsyjJEmSJKknex4lqYrMUg8iliRp6hXLZotHSaqk0K0xkiSVUCibR1s8rgpi3bqRfqUkTaS90crHZqGrmxqNVfvggN2OYpGkVS1NHlYpm+15lKQystTVTUmSpl+tbPZSoyRJkiSpJ3seJamKpNSzpCRJmnrFsnmgnseIOD0i7oqI7RFxwbAaJUlappxd/qISzGZJmjCFsnnZPY8RsRp4L/AqYAdwU0Rszcw7Fv62NXDYhuV+pSTV8fjwb/xIIAtd3dTSLSebV++Bg++aGVUTJWli3bdn+J9ZLZsH6Xk8Cdiemfdk5l7gcuDM4TRLkrRkma1f3ezVqxURB0TEx5r9N0bEUV373t5svysiXj2081Y3s1mSJkmxbB6keDwCuK9rfUezbW5jz4+IbRGxbe/MDwb4OknSOHX1ar0GeD5wTkQ8f85h5wEPZea/B94NvKt57/OBs4EXAKcDf9l8noZrydn8wz3fH1njJEnDNepsHqR4nO8hZU/ok83MzZm5KTM3rVv9pAG+TpLUS87mspc+9NOrdSZwafP6CuDUiIhm++WZuScz/wnY3nyehmvJ2bz2gKeMoFmStHJVyuZBiscdwJFd6xuBnQN8niRpUO3eGtNPr9a/HZOZ+4DvAs/o870anNksSZOmUDYPMmPDTcAxEXE08C90ujx/YbE3PLLngQev+cYf/zOwAXhwgO+eJJ7LZKpyLlXOAzyXuZ4zjIZ0+x4PXfN/84pBZiU7MCK2da1vzszNXev99GotdExfPWIa2JKz+dGHdjz4pSt+12yeXJ7L5KlyHuC5zGU298jmZRePmbkvIt4KXAOsBrZk5u093vNMgIjYlpmblvvdk8RzmUxVzqXKeYDnMgqZeXrLX9FPr9b+Y3ZExBrg6cDuPt+rAZnNHZ7LZKpyLlXOAzyXUaiWzQM95zEzr87Mn8jMf5eZfzTIZ0mSJt6/9WpFxDo6vVpb5xyzFTi3ef164HOZmc32s5sZ344GjgH+YUTtXlHMZklaUUaazcN/0JgkqaSFerUi4iJgW2ZuBT4E/HVEbKdzVfPs5r23R8THgTuAfcBbMtOHC0qSNIBRZ/O4isfNvQ+ZGp7LZKpyLlXOAzyXEjLzauDqOdve0fX6ceANC7z3jwB7wiZXpd9rz2UyVTmXKucBnksJo8zm6PRYSpIkSZK0sIHGPEqSJEmSVoaRF48RcXpE3BUR2yPiglF//yAiYktE7IqI27q2HRIR10bE3c3Pg8fZxn5ExJER8fmIuDMibo+ItzXbp/FcDoyIf4iIf2zO5X8224+OiBubc/lYM4B44kXE6oj4akT8bbM+redxb0R8LSJu2T+99DT+fgFExEERcUVEfL35b+ZnpvVcpIWYzeNnNk8us3nymM3jM9LiMSJWA+8FXgM8HzgnIp4/yjYM6BJg7nS7FwDXZeYxwHXN+qTbB/x2Zh4LvAR4S/P/wzSeyx7gFZl5PHACcHpEvAR4F/Du5lweAs4bYxuX4m3AnV3r03oeAC/PzBO6ps2ext8vgPcAn8nMnwSOp/P/z7Sei/QEZvPEMJsnl9k8eczmMRl1z+NJwPbMvCcz9wKXA2eOuA3Llpk30JmhqNuZwKXN60uBs0baqGXIzPsz8yvN6+/R+Q/uCKbzXDIzv9+srm2WBF4BXNFsn4pziYiNwM8BH2zWgyk8j0VM3e9XRDwNeCmdWcrIzL2Z+TBTeC7SIszmCWA2TyazefKYzeM16uLxCOC+rvUdzbZpdlhm3g+dP/zAoWNuz5JExFHAicCNTOm5NLeT3ALsAq4Fvgk8nJn7mkOm5ffsz4DfA2ab9WcwnecBnX8kfDYibo6I85tt0/j79Vzg28BfNbcsfTAi1jOd5yItxGyeMGbzRDGbJ4/ZPEajLh5jnm1O9zomEfEU4BPAb2TmI+Nuz3Jl5kxmngBspHMF/dj5Dhttq5YmIl4H7MrMm7s3z3PoRJ9Hl5Mz84V0boN7S0S8dNwNWqY1wAuB92XmicCjeBuM6pnmvzXlmM2Tw2yeWGbzGI26eNwBHNm1vhHYOeI2DNsDEXE4QPNz15jb05eIWEsnnD6SmZ9sNk/luezX3LJwPZ2xIgdFxP7nmE7D79nJwBkRcS+dW8ZeQedq57SdBwCZubP5uQu4ks4/HKbx92sHsCMzb2zWr6ATWNN4LtJCzOYJYTZPHLN5MpnNYzTq4vEm4Jhmlqp1wNnA1hG3Ydi2Auc2r88FrhpjW/rS3K//IeDOzPzTrl3TeC7PjIiDmtdPAl5JZ5zI54HXN4dN/Llk5tszc2NmHkXnv4vPZeYbmbLzAIiI9RHx1P2vgdOA25jC36/M/Ffgvoh4XrPpVOAOpvBcpEWYzRPAbJ48ZvNkMpvHKzJH29MeEa+lc9VmNbAlM/9opA0YQERcBpwCbAAeAC4EPgV8HHg28C3gDZk5d+D+RImInwX+DvgaP7qH/w/ojK2YtnM5js6g6NV0LoZ8PDMviojn0rlKeAjwVeAXM3PP+Frav4g4BfidzHzdNJ5H0+Yrm9U1wEcz848i4hlM2e8XQEScQGeihHXAPcCbaX7XmLJzkRZiNo+f2TzZzObJYjaPz8iLR0mSJEnS9Bn1bauSJEmSpClk8ShJkiRJ6sniUZIkSZLUk8WjJEmSJKkni0dJkiRJUk8Wj5IkSZKkniweJUmSJEk9WTxKkiRJknr6/+XiTYfc+O0IAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x7f120f9267b8>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "init()\n",
+    "time_loop(1000)\n",
+    "plot_ρs()"
+   ]
+  }
+ ],
+ "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.6.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/doc/sphinx/tutorials.rst b/doc/sphinx/tutorials.rst
index e66e1b5b182db3f3cdc88c50c8798c9b715c4720..359f101250f95e1c639294520e2034ae14e887d3 100644
--- a/doc/sphinx/tutorials.rst
+++ b/doc/sphinx/tutorials.rst
@@ -14,6 +14,8 @@ You can open the notebooks directly to play around with the code examples.
     /notebooks/04_tutorial_nondimensionalization_and_scaling.ipynb
     /notebooks/05_tutorial_modifying_method_smagorinsky.ipynb
     /notebooks/06_tutorial_thermal_lbm.ipynb
+    /notebooks/07_tutorial_shanchen_twophase.ipynb
+    /notebooks/08_tutorial_shanchen_twocomponent.ipynb
     /notebooks/demo_stencils.ipynb
     /notebooks/demo_create_method_from_scratch.ipynb
     /notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb