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/07_tutorial_shanchen_twocomponent.ipynb b/doc/notebooks/08_tutorial_shanchen_twocomponent.ipynb
similarity index 100%
rename from doc/notebooks/07_tutorial_shanchen_twocomponent.ipynb
rename to doc/notebooks/08_tutorial_shanchen_twocomponent.ipynb
diff --git a/doc/sphinx/tutorials.rst b/doc/sphinx/tutorials.rst
index b1ccbd68c7174c835bd9357da85365ee23593abb..359f101250f95e1c639294520e2034ae14e887d3 100644
--- a/doc/sphinx/tutorials.rst
+++ b/doc/sphinx/tutorials.rst
@@ -14,7 +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_twocomponent.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