diff --git a/doc/notebooks/07_tutorial_shanchen_twocomponent.ipynb b/doc/notebooks/07_tutorial_shanchen_twocomponent.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..81329d401fceb34b3b35ad69af670d2b1907827c
--- /dev/null
+++ b/doc/notebooks/07_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..b1ccbd68c7174c835bd9357da85365ee23593abb 100644
--- a/doc/sphinx/tutorials.rst
+++ b/doc/sphinx/tutorials.rst
@@ -14,6 +14,7 @@ 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/demo_stencils.ipynb
/notebooks/demo_create_method_from_scratch.ipynb
/notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb