07_tutorial_shanchen_twophase.ipynb 34 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
{
 "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",
18
19
    "from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter\n",
    "from lbmpy.maxwellian_equilibrium import get_weights"
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
   ]
  },
  {
   "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",
46
47
48
49
    "rho0 = 1.\n",
    "\n",
    "stencil = get_stencil(\"D2Q9\")\n",
    "weights = get_weights(stencil, c_s_sq=sp.Rational(1,3))"
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data structures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
65
    "dim = len(stencil[0])\n",
66
    "\n",
67
    "dh = ps.create_data_handling((N,)*dim, periodicity=True, default_target='cpu')\n",
68
    "\n",
69
70
    "src = dh.add_array('src', values_per_cell=len(stencil))\n",
    "dst = dh.add_array_like('dst', 'src')\n",
71
    "\n",
72
    "ρ = dh.add_array('rho')"
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
   ]
  },
  {
   "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",
110
111
    "force = sum((psi(ρ[d]) * w_d * sp.Matrix(d)\n",
    "            for d, w_d in zip(stencil, weights)), zero_vec) * psi(ρ.center) * -1 * g_aa"
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kernels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
127
128
129
130
131
132
133
134
135
136
    "collision = create_lb_update_rule(stencil=stencil,\n",
    "                                  relaxation_rate=omega_a, \n",
    "                                  compressible=True,\n",
    "                                  force_model='guo', \n",
    "                                  force=force,\n",
    "                                  kernel_type='collide_only',\n",
    "                                  optimization={'symbolic_field': src})\n",
    "\n",
    "stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'density': ρ})\n",
    "\n",
137
    "\n",
138
139
    "opts = {'cpu_openmp': False, \n",
    "        'target': dh.default_target}\n",
140
    "\n",
141
142
    "stream_kernel = ps.create_kernel(stream, **opts).compile()\n",
    "collision_kernel = ps.create_kernel(collision, **opts).compile()"
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
158
159
160
161
162
163
    "method_without_force = create_lb_method(stencil=stencil, relaxation_rate=omega_a, compressible=True)\n",
    "init_assignments = macroscopic_values_setter(method_without_force, velocity=(0, 0), \n",
    "                                             pdfs=src.center_vector, density=ρ.center)\n",
    "\n",
    "\n",
    "init_kernel = ps.create_kernel(init_assignments, ghost_layers=0).compile()"
164
165
166
167
168
169
170
171
172
173
174
175
   ]
  },
  {
   "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",
176
    "                dh.fill(ρ.name, 2.1, slice_obj=[x,y])\n",
177
    "            else:\n",
178
    "                dh.fill(ρ.name, 0.15, slice_obj=[x,y])\n",
179
    "\n",
180
    "    dh.run_kernel(init_kernel)"
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Timeloop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
196
197
    "sync_pdfs = dh.synchronization_function([src.name])\n",
    "sync_ρs = dh.synchronization_function([ρ.name])\n",
198
199
200
201
202
    "\n",
    "def time_loop(steps):\n",
    "    dh.all_to_gpu()\n",
    "    for i in range(steps):\n",
    "        sync_ρs()\n",
203
    "        dh.run_kernel(collision_kernel)\n",
204
205
    "        \n",
    "        sync_pdfs()\n",
206
    "        dh.run_kernel(stream_kernel)\n",
207
    "        \n",
208
    "        dh.swap(src.name, dst.name)\n",
209
210
211
212
213
    "    dh.all_to_cpu()"
   ]
  },
  {
   "cell_type": "code",
214
   "execution_count": 10,
215
216
217
218
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_ρs():\n",
219
220
    "    plt.title(\"$\\\\rho$\")\n",
    "    plt.scalar_field(dh.gather_array(ρ.name), vmin=0, vmax=2.5)\n",
221
222
223
224
225
226
227
228
229
230
231
232
233
    "    plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the simulation\n",
    "### Initial state"
   ]
  },
  {
   "cell_type": "code",
234
   "execution_count": 11,
235
236
237
238
   "metadata": {},
   "outputs": [
    {
     "data": {
239
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAF0CAYAAAAKF1nQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAeoUlEQVR4nO3de7Bld1Un8O8yCaI8DJgGevKwsSajRkuE6oowTEkkMgakkkwJM2F8RCZWahwZwdGSgFVYWk4V1EyBWghMD0HCFMNjApIeK4AxImiVRDohPJKGSYwKbWLS4a0oVLrX/HFO8Nrc7nvPvvfcc07vz6dq1z17n332WXV376RXr99v/aq7AwAAMAbfsOgAAAAAdooECAAAGA0JEAAAMBoSIAAAYDQkQAAAwGhIgAAAgNGQAAEAAKMhAQIAAEZDAgQwMlX1iKraV1Wfq6r7qurnFx0TAOwUCRDA+LwryZ8neVySy5L896p63GJDAoCdIQECGJGqenaSdPcruvsr3f2HSf46yb9YbGQAsDMkQADjcnGS6x7cqapvSPItSe5dWEQAsIMkQADj8v1JPrNm/+lJ7u/uTy4oHgDYURIggJGoqtOSnJvkOVX10Kr67iSvSfLixUYGADvn1EUHAMCO+a4kf5nk45kMebsvya9397WLDAoAdlJ196JjAGAHVNWPJ/k33f2ji44FABbFEDiA8XhCkoOLDgIAFkkCBDAe35vkE4sOAgA2o6rOrqr3VdXBqrqtql64zjkXVNUXqurW6fayDa9rCBwAALBsqmp3kt3dfUtVPSLJzUku7e7b15xzQZJf7O5nb/a6KkAAAMDS6e57uvuW6esvZTKM+8ytXlcCBAAALLWq2pPkiUluWuftp1TVR6rq3dMlHk5oR9tgn3HGGb1nz56d/EpYOnd85FOLDgFg6Zz7hHMWHQIs1M0333x/d+9adByz+OEffFh/5rNHBn/+5o9+5bYk/7Dm0L7u3nfseVX18CTvSPKi7v7iMW/fkuTbuvtvq+pZSd6VyZp3x7WjCdCePXty4MCBnfxKWDoXPfY/LToEgKXzngOvWXQIsFBV9VeLjmFW93/2SG5671mDP3/a7j//h+7ee6Jzpot4vyPJm7v7nce+vzYh6u7rq+o1VXVGd99/vGsaAgcAACydqqokVyc52N2vPM45j5uel6o6P5P85jMnuu6OVoAAAICTRedIH53nFzw1yU8k+VhV3To99tIk5yRJd78uyXOS/ExVPZDk75Nc1hu0uZYAAQAAM+skRzO/JXW6+0+S1AbnvDrJq2e5rgQIAAAY5GjmWgGaC3OAAACA0VABAgAAZtbpHDnxdJulJAECAAAGmeccoHmRAAEAADPrJEckQAAAwFisYgVIEwQAAGA0VIAAAICZdaIJAgAAMB6rtwrQJofAVdXpVXVtVX2iqg5W1VOq6tFVdUNV3TH9+ah5BwsAACyHTufIFrZF2ewcoN9M8p7u/s4kT0hyMMlVSW7s7nOT3DjdBwAAxqCTI1vYFmXDBKiqHpnkB5JcnSTd/dXu/nySS5JcMz3tmiSXzitIAACA7bCZCtC3Jzmc5Heq6sNV9fqqeliSx3b3PUky/fmY9T5cVVdW1YGqOnD48OFtCxwAAFiczmQO0NBtUTaTAJ2a5ElJXtvdT0zyd5lhuFt37+vuvd29d9euXQPDBAAAlkvlyBa2RdlMAnQoyaHuvmm6f20mCdG9VbU7SaY/75tPiAAAwLLpJEd7+LYoGyZA3f03ST5dVd8xPXRhktuT7E9y+fTY5Umum0uEAAAA22Sz6wD95yRvrqqHJLkryfMzSZ7eXlVXJPlUkufOJ0QAAGAZLXIo21CbSoC6+9Yke9d568LtDQcAAFgFnZM4AQIAADjW0ZYAAQAAI7CqFaDNdIEDAAA4KagAAQAAM+tUjqxgPUUCBAAADGIOEAAAMAqrOgdIAgQAAAxQOdKrNwRu9SIGAAAYSAUIAACYWSc5uoL1FAkQAAAwiDlAAADAKHSbAwQAALDUVIAAAIBBjhoCBwAAjMFkHaDVG1AmAQIAAAZYzTlAEiAAAGBmq9oGe/UiBgAAGEgFCAAAGORIa4IAAACMQKc0QQAAAMbjqCYIAADAGKxqG+zVixgAAGAgFSAAAGBmndIEAQAAGI9VXAdIAgQAAMysOzmygk0QVi9iAACAgVSAAACAASpHYw4QAAAwAp3VHAInAQIAAAZZxXWAJEAAAMDMOpWjK9gGe/VSNgAAgIFUgAAAgEEMgQMAAEahkxzVBAEAABiHyhFtsAEAgDFY1QrQ6kUMAAAwkAoQAAAwiCFwAADAKHSXIXAAAMB4HOlvGLxtpKrOrqr3VdXBqrqtql64zjlVVb9VVXdW1Uer6kkbXVcFCAAAWEYPJPmF7r6lqh6R5OaquqG7b19zzjOTnDvdvj/Ja6c/j0sFCAAAmFknOZoavG14/e57uvuW6esvJTmY5MxjTrskyZt64oNJTq+q3Se6rgoQAAAwQG1qKNu2fFPVniRPTHLTMW+dmeTTa/YPTY/dc7xrbSoBqqq/TPKlJEeSPNDde6vq0UnelmRPkr9M8m+7+3ObuR4AALDaJusAbakL3BlVdWDN/r7u3nfsSVX18CTvSPKi7v7isW8fJ7TjmqUC9IPdff+a/auS3NjdL6+qq6b7L57hegAAwAo7srUZNfd3994TnVBVp2WS/Ly5u9+5zimHkpy9Zv+sJHef6JpbifiSJNdMX1+T5NItXAsAAOBrqqqSXJ3kYHe/8jin7U/yk9NucE9O8oXuPu7wt2TzFaBO8vtV1Un+x7Q09dgHL97d91TVY44T+JVJrkySc845Z5NfBwAALLNObXUI3EaemuQnknysqm6dHntpknOSpLtfl+T6JM9KcmeSLyd5/kYX3WwC9NTuvnua5NxQVZ/YbNTTZGlfkuzdu/eE4/EAAIDVcXSOTaW7+0+y/hyfted0kp+d5bqbSoC6++7pz/uq6neTnJ/k3qraPa3+7E5y3yxfDAAArK7u5Mh8K0BzsWHKVlUPmy48lKp6WJJ/neTjmYy3u3x62uVJrptXkAAAwPI52jV4W5TNVIAem+R3J3OQcmqS/93d76mqDyV5e1VdkeRTSZ47vzABAAC2bsMEqLvvSvKEdY5/JsmF8wgKAABYbpMmCDuzEOp2mmUdIAAAgK85cuIeBUtJAgQAAMysk4XO5Rlq9WpWAAAAA6kAAQAAA5gDBAAAjMhRc4AAAIAxWNWFUCVAAADAIIbAAbAyLn7/7ese3/+083Y4EgDYORIgAABgZpOFUA2BAwAARkITBAAAYBQshAoAALDkVIAAltDxGhScLN+t0QLAyUEXOAAAYBxaEwQAAGAkOpogAAAAI7KKFaDVG7QHAAAwkAoQwA5aZHODZbLZ34NmCQDLa1XbYEuAAACAQSRAAADAKHR0gQMAAEZkFbvAaYIAAACMhgoQwDbQ3GA+Zvm9apgAsMPaHCAAAGAkdIEDAABGZRUTIHOAAACA0VABAgAAZqYNNsAIaHawvNa7NxojAMxXS4AAAICxWMV1gCRAAADAzHpF22BrggAAAIyGChAAADCIOUAAJxEND1afxggA86QLHAAAMCIqQAAAwCh0NEEAAABYaipAAADA7HrSCnvVSIAAouHBmGiMALB9LIQKAACMQmc1myCYAwQAAIyGChAAADCAdYAAAIARWcUmCJseAldVp1TVh6vq96b7j6+qm6rqjqp6W1U9ZH5hAgAAy6a7Bm+LMsscoBcmObhm/xVJXtXd5yb5XJIrtjMwAABgeXWfxAlQVZ2V5EeSvH66X0menuTa6SnXJLl0HgECAABsl83OAfqNJL+U5BHT/W9N8vnufmC6fyjJmet9sKquTHJlkpxzzjnDIwUAAJbKKjZB2LACVFXPTnJfd9+89vA6p647Baq793X33u7eu2vXroFhAgAAy2YyDG7YtiibqQA9NcnFVfWsJA9N8shMKkKnV9Wp0yrQWUnunl+YAADAslnFhVA3TIC6+yVJXpIkVXVBkl/s7h+rqv+T5DlJ3prk8iTXzTFOgG1z8ftvX3QILJn1/kzsf9p5C4gEYHV0FtvMYKhZusAd68VJ/ktV3ZnJnKCrtyckAACA+ZhpIdTu/qMkfzR9fVeS87c/JAAAYBWs4DqosyVAAAAASZI+SecAAQAArGsFS0BbmQMEAAAwF1X1hqq6r6o+fpz3L6iqL1TVrdPtZZu5rgoQAAAwyJyHwL0xyauTvOkE5/xxdz97lotKgAAAgEHmuaBpd3+gqvZs93UNgQMAAGbWmVSAhm5JzqiqA2u2KweE8ZSq+khVvbuqvnszH1ABAgAAZtdJtjYE7v7u3ruFz9+S5Nu6+2+r6llJ3pXk3I0+JAECTmoXv//2RYfAilrvz87+p523gEgAWE93f3HN6+ur6jVVdUZ333+iz0mAAACAQeY5B2gjVfW4JPd2d1fV+ZlM7/nMRp+TAAEAAMPMMQGqqrckuSCTuUKHkvxKktOSpLtfl+Q5SX6mqh5I8vdJLuveOCWTAAEAAAPUXNtgd/fzNnj/1Zm0yZ6JBAgAABhmgUPghtIGGwAAGA0VIAAAYHaduQ6BmxcJEAAAMMwKDoGTAAEAAAOtXgXIHCAAAGA0VIAAAIBhDIEDAABGQwIEAACMQifRBQ4AABiLXsEKkCYIAADAaKgAAQAAw6xgBUgCBAAADGMOEAAAMBalAgQAAIxCZyWHwGmCAAAAjIYKEAAAMECZAwQAAIzICg6BkwABAADDrGACZA4QAAAwGipAAADAMCtYAZIAASeNi99/+6JD4CS33p+x/U87bwGRACyBjiYIAADAeFgIFQAAGI8VTIA0QQAAAEZDAgQAAIyGIXAAAMAg5gABLNB63bh0hmM76fgGcAxd4AAAgFHoaIIAAACwzFSAAACAYU7GClBVPbSq/qyqPlJVt1XVr06PP76qbqqqO6rqbVX1kPmHCwAALIvq4duibGYI3FeSPL27n5Dk+5JcVFVPTvKKJK/q7nOTfC7JFfMLEwAAWDq9hW1BNkyAeuJvp7unTbdO8vQk106PX5Pk0rlECAAAsE021QShqk6pqluT3JfkhiR/nuTz3f3A9JRDSc48zmevrKoDVXXg8OHD2xEzAACwDE7GClCSdPeR7v6+JGclOT/Jd6132nE+u6+793b33l27dg2PFAAAWBpbmf+zyDlAM3WB6+7PV9UfJXlyktOr6tRpFeisJHfPIT4AAGBZreBCqJvpArerqk6fvv6mJD+U5GCS9yV5zvS0y5NcN68gAQCAJbSCQ+A2UwHaneSaqjolk4Tp7d39e1V1e5K3VtWvJ/lwkqvnGCcAAMCWbZgAdfdHkzxxneN3ZTIfCAAAGKFFzuUZaqY5QAAAAF8jAQIAAEZhwd3chtpUG2wAAICTgQoQAAAwzApWgCRAAADAMBIgAABgLMwBAgAAWGISIAAAYDQMgQMAAIZZwSFwEiAAAGB2K7oOkAQIAAAYRgIEAACMhgQIYLnsf9p5X3fs4vffvoBIWDXr/dkBYPVJgAAAgJlVzAECAADGRAIEAACMwop2gbMQKgAAsHSq6g1VdV9Vffw471dV/VZV3VlVH62qJ23muhIgAABgmN7CtrE3JrnoBO8/M8m50+3KJK/dzEUlQAAAwDBzTIC6+wNJPnuCUy5J8qae+GCS06tq90bXNQcIAAAYZItzgM6oqgNr9vd1974ZPn9mkk+v2T80PXbPiT4kAQIAAIbZWgJ0f3fv3cLna51jG0ZkCBwAALCKDiU5e83+WUnu3uhDKkDA6Ox/2nlfd+zi99++gEhYFuv9mQBgA5tvZjAv+5O8oKremuT7k3yhu084/C2RAAEAAAPNcx2gqnpLkgsymSt0KMmvJDktSbr7dUmuT/KsJHcm+XKS52/muhIgAABgmDkmQN39vA3e7yQ/O+t1JUAAAMAg86wAzYsmCAAAwGioAAEAAMOsYAVIAgQAAMxu8V3gBpEAAQAAM6usvxLpsjMHCAAAGA0VIAAAYBhD4ABW0/6nnfd1xy5+/+0LiIR5W+9eAzDMKrbBlgABAADDSIAAAIDRWMEESBMEAABgNFSAAACA2bU5QAAAwJhIgABOHjrDrT4d3wDmSwUIAAAYjxVMgDRBAAAARkMFCAAAGGQVh8BtWAGqqrOr6n1VdbCqbquqF06PP7qqbqiqO6Y/HzX/cAEAgKXQW9wWZDMVoAeS/EJ331JVj0hyc1XdkOSnktzY3S+vqquSXJXkxfMLFWDxjjepXnOExdPwAGABTsYKUHff0923TF9/KcnBJGcmuSTJNdPTrkly6byCBAAA2A4zzQGqqj1JnpjkpiSP7e57kkmSVFWPOc5nrkxyZZKcc845W4kVAABYEpWTdA7Qg6rq4UnekeRF3f3FzX6uu/d1997u3rtr164hMQIAAMvoJJ0DlKo6LZPk583d/c7p4Xurave0+rM7yX3zChIAAFg+1atXAtowAaqqSnJ1koPd/co1b+1PcnmSl09/XjeXCAFWwGYn4GuWMBuNDQCW2IIrOUNtpgL01CQ/keRjVXXr9NhLM0l83l5VVyT5VJLnzidEAACA7bFhAtTdf5LJHKf1XLi94QAAAKtiFZsgzNQFDgAA4GskQAAAwFioAAFwQpolTGhuAHCSWMEEaNPrAAEAAKw6FSAAAGB2bQgcAAAwJhIgAABgDCoqQABsk51oEnC8RgsaFABwMpMAAQAAw/TqlYAkQAAAwCCGwAEAAOPQ0QQBAAAYjzq66AhmJwECGCnNDgAYIwkQAAAwjCFwAADAWGiCAAAAjENHG2wAAGA8VrEC9A2LDgAAAGCnqAABAADDrGAFSAIEAADMrLKaQ+AkQAAAwOy6V7IJgjlAAADAaKgAAQAAgxgCBwAAjIcECAAAGAsVIAAAYBw6ydHVy4A0QQAAAEZDBQgAABhm9QpAEiAAAGAYc4AAAIDxsBAqAAAwFtXDt01dv+qiqvpkVd1ZVVet8/5PVdXhqrp1uv30RtdUAQIAAJZOVZ2S5LeTPCPJoSQfqqr93X37Mae+rbtfsNnrqgABAACz6y1uGzs/yZ3dfVd3fzXJW5NcstWwJUAAAMDMKkl1D9424cwkn16zf2h67Fg/WlUfraprq+rsjS4qAQIAAIY5uoUtOaOqDqzZrjzm6rXONx6bOf3fJHu6+3uT/EGSazYK2RwgAABgEe7v7r0neP9QkrUVnbOS3L32hO7+zJrd/5nkFRt9qQoQAAAwyJyHwH0oyblV9fiqekiSy5Ls/yffX7V7ze7FSQ5udFEVIAAAYHabb2Yw7PLdD1TVC5K8N8kpSd7Q3bdV1a8lOdDd+5P8XFVdnOSBJJ9N8lMbXVcCBAAADNBzXwi1u69Pcv0xx1625vVLkrxklmtKgAAAgEE2u6DpMjEHCAAAGA0VIAAAYJg5D4Gbhw0rQFX1hqq6r6o+vubYo6vqhqq6Y/rzUfMNEwAAWCqd1NHh26JsZgjcG5NcdMyxq5Lc2N3nJrlxug8AAIxJ9/BtQTZMgLr7A5m0lFvrkvzjKqvXJLl0m+MCAADYdkObIDy2u+9JkunPxxzvxKq6sqoOVNWBw4cPD/w6AABg6fQWtgWZexe47t7X3Xu7e++uXbvm/XUAAMAOqe7B26IMTYDurardSTL9ed/2hQQAAKyEk3EO0HHsT3L59PXlSa7bnnAAAICV0EmObmFbkM20wX5Lkj9N8h1Vdaiqrkjy8iTPqKo7kjxjug8AALDUNlwItbufd5y3LtzmWAAAgBVRWexcnqE2TIAAAADWJQECAABGQwIEAACMwoNNEFbM3NcBAgAAWBYqQAAAwCCaIAAAAOMhAQIAAMahVzIBMgcIAAAYDRUgAABgdp2VrABJgAAAgGFWsA22BAgAABhEFzgAAGA8VjAB0gQBAAAYDRUgAABgdp3k6OpVgCRAAADAAKu5DpAECAAAGEYCBAAAjMYKJkCaIAAAAKOhAgQAAMxOEwQAAGA8Oumjiw5iZhIgAABgGHOAAAAAlpcKEAAAMDtzgAAAgFFZwSFwEiAAAGAYCRAAADAOvZIJkCYIAADAaKgAAQAAs+skR60DBAAAjMUKDoGTAAEAAMNIgAAAgHHolVwHSBMEAABgNFSAAACA2XXSrQkCAAAwFis4BE4CBAAADLOCTRDMAQIAAEZDBQgAAJhdt4VQAQCAEVnBIXASIAAAYJBWAQIAAMahV7ICpAkCAAAwGipAAADA7DoruQ7QlipAVXVRVX2yqu6sqqu2KygAAGAF9NHh24IMrgBV1SlJfjvJM5IcSvKhqtrf3bdvV3AAAMBy6iQ9sgrQ+Unu7O67uvurSd6a5JLtCQsAAFhq3XOvAG004qyqvrGq3jZ9/6aq2rPRNbeSAJ2Z5NNr9g9Njx0b1JVVdaCqDhw+fHgLXwcAAIzFmhFnz0xyXpLnVdV5x5x2RZLPdfc/T/KqJK/Y6LpbSYBqnWNfVwPr7n3dvbe79+7atWsLXwcAACyTPtqDt03YzIizS5JcM319bZILq2q9POVrtpIAHUpy9pr9s5LcvYXrAQAAq2S+Q+A2M+Lsa+d09wNJvpDkW0900a20wf5QknOr6vFJ/jrJZUn+/Yk+cPPNN99fVX+1he9kmDOS3L/oIPg67stycl+Wk/uynLbtvlS9djsuw4TnZTltdF++bacC2S5fyufe+wd97RlbuMRDq+rAmv193b1vzf5mRpxtalTaWoMToO5+oKpekOS9SU5J8obuvm2DzxgDtwBVdaC79y46Dv4p92U5uS/LyX1ZTu7LcnJfltPJeF+6+6I5f8VmRpw9eM6hqjo1ybck+eyJLrqlhVC7+/ok12/lGgAAAOvYzIiz/UkuT/KnSZ6T5A+7ez4VIAAAgHk53oizqvq1JAe6e3+Sq5P8r6q6M5PKz2UbXVcCNA77Nj6FBXBflpP7spzcl+Xkviwn92U5uS8DrDfirLtftub1PyR57izXrA0qRAAAACeNrbTBBgAAWCkSoJNYVf23qvpEVX20qn63qk5f895LqurOqvpkVf3wIuMco6q6aPq7v7Oqrlp0PGNVVWdX1fuq6mBV3VZVL5wef3RV3VBVd0x/PmrRsY5NVZ1SVR+uqt+b7j++qm6a3pO3VdVDFh3jGFXV6VV17fT/LQer6imel8Wqqp+f/vfr41X1lqp6qOdlMarqDVV1X1V9fM2xdZ+Pmvit6d8DPlpVT1pc5OMjATq53ZDke7r7e5P8vyQvSZKqOi+TCWLfneSiJK+pqlMWFuXITH/Xv53kmUnOS/K86T1h5z2Q5Be6+7uSPDnJz07vxVVJbuzuc5PcON1nZ70wycE1+69I8qrpPflckisWEhW/meQ93f2dSZ6QyT3yvCxIVZ2Z5OeS7O3u78lkkvhl8bwsyhsz+XvVWsd7Pp6Z5NzpdmUSC2HtIAnQSay7f3+6Im6SfDCT3ulJckmSt3b3V7r7L5LcmeT8RcQ4UucnubO77+ruryZ5ayb3hB3W3fd09y3T11/K5C9zZ2ZyP66ZnnZNkksXE+E4VdVZSX4kyeun+5Xk6UmunZ7inixAVT0yyQ9k0nEp3f3V7v58PC+LdmqSb5quf/LNSe6J52UhuvsD+fr1Z473fFyS5E098cEkp1fV7p2JFAnQePyHJO+evj4zyafXvHdoeoyd4fe/hKpqT5InJrkpyWO7+55kkiQlecziIhul30jyS0mOTve/Ncnn1/yDjmdmMb49yeEkvzMdnvj6qnpYPC8L091/neS/J/lUJonPF5LcHM/LMjne8+HvAgskAVpxVfUH03G/x26XrDnnlzMZ6vPmBw+tcyntAHeO3/+SqaqHJ3lHkhd19xcXHc+YVdWzk9zX3TevPbzOqZ6ZnXdqkicleW13PzHJ38Vwt4Wazie5JMnjk/yzJA/LZGjVsTwvy8d/1xbIOkArrrt/6ETvV9XlSZ6d5MI1q+IeSnL2mtPOSnL3fCJkHX7/S6SqTssk+Xlzd79zevjeqtrd3fdMhyTct7gIR+epSS6uqmcleWiSR2ZSETq9qk6d/qu2Z2YxDiU51N03TfevzSQB8rwszg8l+YvuPpwkVfXOJP8ynpdlcrznw98FFkgF6CRWVRcleXGSi7v7y2ve2p/ksqr6xqp6fCYT8P5sETGO1IeSnDvt0vOQTCas7l9wTKM0nVtydZKD3f3KNW/tT3L59PXlSa7b6djGqrtf0t1ndfeeTJ6NP+zuH0vyviTPmZ7mnixAd/9Nkk9X1XdMD12Y5PZ4XhbpU0meXFXfPP3v2YP3xPOyPI73fOxP8pPTbnBPTvKFB4fKMX8WQj2JVdWdSb4xyWemhz7Y3f9x+t4vZzIv6IFMhv28e/2rMA/Tf93+jUw69ryhu//rgkMapar6V0n+OMnH8o/zTV6ayTygtyc5J5O/YDy3u4+d2MqcVdUFSX6xu59dVd+eScOQRyf5cJIf7+6vLDK+Maqq78ukOcVDktyV5PmZ/GOq52VBqupXk/y7TP5//uEkP53JXBLPyw6rqrckuSDJGUnuTfIrSd6VdZ6PacL66ky6xn05yfO7+8Ai4h4jCRAAADAahsABAACjIQECAABGQwIEAACMhgQIAAAYDQkQAAAwGhIgAABgNCRAAADAaEiAAACA0fj/0wN5PNSa6lIAAAAASUVORK5CYII=\n",
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
      "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",
272
   "execution_count": 12,
273
274
275
276
277
   "metadata": {},
   "outputs": [],
   "source": [
    "init()\n",
    "time_loop(1)\n",
278
279
280
    "ref = 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",
    "\n",
    "assert np.allclose(dh.gather_array(ρ.name)[N//2], ref)"
281
282
283
284
285
286
287
288
289
290
291
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the simulation until converged"
   ]
  },
  {
   "cell_type": "code",
292
   "execution_count": 13,
293
294
295
296
   "metadata": {},
   "outputs": [
    {
     "data": {
297
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAF0CAYAAAAKF1nQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df6ycV33n8c9nZu4PO3YIiU1w86NOtdm2tIIGWQGWVZWSsg00Ilk17Ibtj5RNFW1VWugPlYRKrYq6EmgrYLspYb0kxaxYAhso8VahNE2pKFJJcUIaSAyNm7aJG9exnR+2c+17PTPf/WMm6a3PGd9nZu7cZ577vF/S6HrOfeZ5zsxzZq7PnHM+jyNCAAAAAFAHjbIrAAAAAABrhQ4QAAAAgNqgAwQAAACgNugAAQAAAKgNOkAAAAAAaoMOEAAAAIDaoAMEAAAAoDboAAEAAACoDTpAAFAztjfb3mn7WdtP2/7lsusEAMBaoQMEAPXzBUl/K+mVkq6X9Lu2X1lulQAAWBt0gACgRmxfLUkR8cGIWIyIP5P0j5L+dbk1AwBgbdABAoB6eZuku1+8Y7sh6WWSDpZWIwAA1hAdIACol9dJOrLs/pskHY6I75RUHwAA1hQdIACoCdszki6VdJ3teds/IOmjkt5bbs0AAFg7rbIrAABYM98v6e8lfUu9KW9PS/qdiLirzEoBALCWHBFl1wEAsAZs/5Skfx8RP1F2XQAAKAtT4ACgPl4jaW/ZlQAAoEx0gACgPl4t6dtlVwIAgCJsX2T7y7b32n7E9rsz21xh+3nbD/Vvv7nifpkCBwAAAGDa2N4maVtEPGh7s6QHJF0bEY8u2+YKSb8WEVcX3S8jQAAAAACmTkQciIgH+/8+pt407gvG3S8dIAAAAABTzfZ2SZdJuj/z6zfY/mvbX+xf4uGM1jQGe8uWLbF9+/a1PCQwdR57+Imyq4Cp47IrMCWYkl1nl7764rKrAJTqgQceOBwRW8uuxzB+7EfOiiPPdEZ+/AMPLz4i6eSyop0RsfP07WxvkvQ5Se+JiKOn/fpBSd8dEcdtv1XSF9S75t1Aa9oB2r59u/bs2bOWhwSmzlsu+MWyq4Bp02AwXpLU7ZZdA5Toi3v+R9lVAEpl+x/KrsOwDj/T0f1funDkx89s+9uTEbHjTNv0L+L9OUmfiojPn/775R2iiLjH9kdtb4mIw4P2yV9dAAAAAFPHtiXdLmlvRHxowDav7G8n25er1785cqb9rukIEAAAAID1ItSJiY7ev1HST0v6pu2H+mXvk3SxJEXExyRdJ+nnbbclnZB0fawQc00HCADKlpv6td6nxTHdDQAqLyR1J7h+MyK+qhUWykbErZJuHWa/dIAAAAAAjKSr6n2htc6/YgQAAACAf8YIEAAAAIChhUKdMy+3mUp0gAAAAACMZJJrgCaFDhAAAACAoYWkDh0gAAAAAHVRxREgQhAAAAAA1AYjQAAAAACGFhIhCAAAAADqo3pXASo4Bc72Obbvsv1t23ttv8H2ubbvtf1Y/+fLJ11ZAAAAANMhFOqMcStL0TVA/13SH0fE90l6jaS9km6WdF9EXCrpvv59AAAAAHUQUmeMW1lW7ADZPlvSD0u6XZIiYikinpN0jaRd/c12Sbp2UpUEAAAAgNVQZAToeyQdkvQHtr9h++O2z5J0fkQckKT+z1fkHmz7Jtt7bO85dOjQqlUcAAAAQHlCvTVAo97KUqQD1JL0Wkm3RcRlkl7QENPdImJnROyIiB1bt24dsZoAAAAApovVGeNWliIdoP2S9kfE/f37d6nXITpoe5sk9X8+PZkqAgAAAJg2Iakbo9/KsmIHKCL+SdKTtr+3X3SlpEcl7ZZ0Q7/sBkl3T6SGAAAAALBKil4H6Bclfcr2rKTHJb1Tvc7TZ23fKOkJSW+fTBUBAAAATKMyp7KNqlAHKCIekrQj86srV7c6AAAAAKogtI47QAAAAABwum7QAQIAAABQA4wAAQBWTzdzhYRGkeDOKZR7LgAAlIQOEAAAAIChhaxOoavqTBc6QAAAAABGwhogAAAAALXAGiAAAAAANWJ1gilwAIBJqUIwAoEHAIApRwcIAAAAwNBCUpcQBAAAAAB1wRogAAAAALUQUc01QNWrMQAAAACMiBEgAAAAACPpMgUOAAAAQB30rgNUvQlldIAAAAAAjKCaa4DoAAEAAAAYWlVjsKtXYwAAAAAYESNAAFAVjQp8Z5WrY7e79vUAAKyJThCCAAAAAKAGQiYEAQAAAEB9dAlBAAAAAFAHVY3Brl6NAQAAAGBEjAABwFqqQpDBahvnOROgAABTK2RCEAAAAADURxWvA0QHCAAAAMDQIqROBUMQqldjAAAAABgRI0AAAAAARmB1xRogAKinOoYbrIVhXlcCEwBgTYWqOQWODhAAAACAkVTxOkB0gAAAAAAMLWR1KxiDXb0uGwAAAACMiBEgAAAAACNhChwArHfTFnbgKZt6EFHesXPnhmAEAJiYkNQlBAEAAABAPVgdYrABAAAA1EFVR4CqV2MAAAAAGBEjQAAAAABGwhQ4AKiqMsMNckEGubJGwT8yRfc3jFy4QdHAg27Bx04iQKHoeSUsAQCGFmGmwAEAAACoj040Rr6txPZFtr9se6/tR2y/O7ONbf+e7X22H7b92pX2ywgQAAAAgGnUlvSrEfGg7c2SHrB9b0Q8umybt0i6tH97naTb+j8HYgQIAAAAwNBCUlce+bbi/iMORMSD/X8fk7RX0gWnbXaNpE9Gz9cknWN725n2ywgQAAAAgBG40FS2VTmSvV3SZZLuP+1XF0h6ctn9/f2yA4P2VagDZPvvJR2T1JHUjogdts+V9BlJ2yX9vaT/EBHPFtkfAJRqEoEH4wQZZOoTzUwdm820rJWWRSvz2ExdYkAwggsGFLidCQ5od9KyTlrmTuaxuSCCogEKZyovItcmCEYAgDPqXQdorJCdLbb3LLu/MyJ2nr6R7U2SPifpPRFx9PRfD6jaQMOMAP1IRBxedv9mSfdFxAds39y//94h9gcAAACgwjrjrag5HBE7zrSB7Rn1Oj+fiojPZzbZL+miZfcvlPTUmfY5To2vkbSr/+9dkq4dY18AAAAA8BLblnS7pL0R8aEBm+2W9DP9NLjXS3o+IgZOf5OKjwCFpD+xHZL+Z39o6vwXdx4RB2y/YkDFb5J0kyRdfPHFBQ8HAAAAYJqFPO4UuJW8UdJPS/qm7Yf6Ze+TdLEkRcTHJN0j6a2S9klakPTOlXZatAP0xoh4qt/Judf2t4vWut9Z2ilJO3bsmMBV7gAAAACUoTvBUOmI+Krya3yWbxOSfmGY/RbqAEXEU/2fT9v+Q0mXSzpoe1t/9GebpKeHOTAAAACA6oqQOpMdAZqIFTtAts+S1IiIY/1//ztJ71dvvt0Nkj7Q/3n3JCsKACuaRLrb6QYkp6lgalvMzRQq62xIy7rz6f7ambLubFrHbiuXUpcW9SqUFjXaaWFjKS1rnUwT3xqZsuaJU2l1FouV5VLleuWZ1LbVTobLIS0OQI1NeArcRBQZATpf0h/21iCpJen/RMQf2/66pM/avlHSE5LePrlqAgAAAMD4VuwARcTjkl6TKT8i6cpJVAoAAADAdOuFIKzNhVBX0zDXAQIAAACAl3TOnFEwlegAAQAAABhaaP2uAQKA6bPagQe5cINcWSsNHZCkmJ9Nyrob07L25rmkbOll6Ufx0ub0+S1tSuvT3pgJPEgPq26u2oNewsya/kYmd6CxlJa1FtLnMns8DSKYPZa+DrPPt9P9HVtMj7uQObAkn8yUtzMVzwUjrHZYAsEIADC16AABAAAAGAFrgAAAAADUSJc1QAAAAADqYN1eCBUAAAAAcpgCBwCTsBaBB830GDE7k5ZtTBfvS1L77Pmk7OSWNI3gxHnpcRbPTeuztDldlN/ZmJZ159PF9tHKLMBvZhb5D/rSLpcH0Ek3djt9Lo2T6XbNhbRs9lj652fumTSpYcOR9BzMH86HILSOnkzruJCGKHjpVPrgTuY1IxgBANYlOkAAAAAAhta7ECpT4AAAAADUBCEIAAAAAGqhqhdCrd6qJQAAAAAYESNAANa3XOBBK11snws86J69ISlb3JKWSdIL56cfpyfOT4998tx0YX377Ha6ww2dpKg5l5bNNNOF9Y1Gegy7WJkkRebbvFxZt5uWdTrp92pLi+nrvXQiLTt5blq2mCnbkAmckKSzDmaCFQ6nj28cPZGUZYMR2unrPVYwAgCsQ6TAAQAAAKiHIAQBAAAAQE2ECEEAAAAAUCNVHAGq3qQ9AAAAABgRI0AApktjjO9ligYezM0mZZ2Xb0zKTpyfLrY//l3p/iRpYVtatnReZhH9pnSx/cxcGoIwM5MJQcgEHjSHCDcoul0u8KDodp1c2Vx6Tk9tSF/H9lnpn6QXNqfhFKc25dtIe2P6+E1z6Xnd0Eof33x2ISmzljIHGSMYYVDb7qbnFQCqoKox2HSAAAAAAIyEDhAAAACAWgiRAgcAAACgRqqYAkcIAgAAAIDaYAQIQDnGCTuQ8oEHzXSfMZsuos8FHixs25CUHb04Xai/sC2/4L19Xhpk0DorDTyYnUvLWrlwg0aubLxwg+KKPT4XgpD7o9LJPJfcc25ngh+WZtLtFjPnVJI6s+nRc2XRSM912iKk5jOZgIlc4EEnE2JQNBhByr8XCEYAUAXBGiAAAAAANUEKHAAAAIBaqWIHiDVAAAAAAGqDESAAAAAAQyMGGwDWUiYEIRd40D07XfB+4vz5pCwbeHBBupC9vSUNMZCkuU2LSdnsbCYYoWC4wfhBBv9SY8z95f7AFa1jq5lul3vO2eCHTFjCUisNS5CkxeZcUrbQyAUmpOfa3bRNbGynx248nzlXJ5fSQwwTggAAFZYLxJl2dIAAAAAAjKSK1wGiAwQAAABgaFHRGGxCEAAAAADUBiNAAAAAAEbCGiAAyMld6X4YmcADtdKF7LExXQS/uCUNQTj+XZnAg23FAg9yYQeSNJcLPGimi/UbmadSNExg3CCDcYxz7KIBCs1MM7HT19CzAw60KS3Kna2FbhqM0FxK20RzMW0786cy9WlnQhlOZV6vYYIRcu+ZbhrKAADlIgUOAAAAQI0wAgQAAACgFkKEIAAAAADAVGMECAAAAMDwoprXfaYDBGB1TSLwILM6PubTlfDts+eTshfOTz/mFralh2ifl4YY5AIPcmEHUj7woNmY/nCDtVD0+eWW+DezMysyoQOSlAtHyAUjdNKdLiyl7aS1kCl7IW1jM0tpm3AusKAzIMSg6P8eCEYAMIW4ECoAAACAWghVMwSBNUAAAAAAaoMRIAAAAAAj4DpAAAAAAGqkiiEIhafA2W7a/obtP+rfv8T2/bYfs/0Ze+C1uQEAAACsQxEe+VaWYUaA3i1pr6Sz+/c/KOnDEXGn7Y9JulHSbatcPwB1k02BayZF3Y3pdy4nt6RlJ85P97d0Xpoi1jrrVFI2m0l8y6W9SSS+rYbca5ObWjH4tU7PTWS+mutkzvXSeen3gSeOp+1u7mi6w+ZCpmwxPYa6A+pdxa9PAUC9j691G4Jg+0JJPy7p4/37lvQmSXf1N9kl6dpJVBAAAAAAVkvREaCPSPp1SZv798+T9FxEvPj16H5JF+QeaPsmSTdJ0sUXXzx6TQEAAABMlSqGIKw4AmT7aklPR8QDy4szm2bH8CNiZ0TsiIgdW7duHbGaAAAAAKZNbxrcaLeyFBkBeqOkt9l+q6R59dYAfUTSObZb/VGgCyU9NblqAgAAAJg2VVwDtGIHKCJukXSLJNm+QtKvRcRP2v6/kq6TdKekGyTdPcF6AphGjQlcS7mRfpDG3ExS1t48l5SdyCxkP3lu5iumTZnAg7m0rNXoFqneQAQejK9oMEJv27Qsdw5z5/rEplzbyQQjZNrY3LNpW2wsLCVl7uQDNJRWsbjce7A7zg4BoLhQuWluoxrnfy/vlfQrtveptybo9tWpEgAAAABMxlAXQo2IP5f05/1/Py7p8tWvEgAAAIAqqOJch6E6QAAAAAAgSarodYDoAAEAAAAYTQWHgOgAAWutaHDAel/I7AHfGGVen1wIwtLL0o+vxXPTfbbPbidlM3NpWauZvt7NRvqp7gHBBtMUeLBW12Qo6zkPOm7uHdPMvN1y57qVaROnzs61sbQs1xZnnkvbrE+mwQi9X2RqXmY+7FqYRIAKgHXH9h2SXrwkzw9mfn+FekFsf9cv+nxEvH+l/dIBAgAAADCSCU+B+4SkWyV98gzb/EVEXD3MTukAAQAAABjJJAesI+Irtrev9n4ZgwYAAAAwtFBvBGjUm6Qttvcsu900QjXeYPuvbX/R9g8UeQAjQAAAAACGF5LGmwJ3OCJ2jPH4ByV9d0Qct/1WSV+QdOlKD6IDBEyrSSwSnqZghQEhCJFZtd7ZkAlB2Jxut7Q5Mw6/oZMUzcykZc1G+toMCjwoy1qFGxRVtD5lBkTkzmHuXOfaxKlM21na3MyUpW1xPtNmfTz/nnY78zpOUwgCgQUAplREHF3273tsf9T2log4fKbH0QECAAAAMJIyv6+x/UpJByMibF+u3vKeIys9jg4QAAAAgNFMsANk+9OSrlBvrdB+Sb8laUaSIuJjkq6T9PO225JOSLo+YuUuGR0gAAAAACPwRGOwI+IdK/z+VvVisodCBwgAAADAaKZoyWJRdICAOsktZi4ajLDaC6EbA74xaqaLzLvzmYXnm9LHdzZmFrzPZQIPmulzbjaKfYKv1YL+aQs8GEfuuUzidcztM3fs3LnOtolM28m1sVxbzLXZXNuWJDXaadk4eSXT9D4HgClEBwgAAADA8EITnQI3KXSAAAAAAIyGKXAAAAAA6qN6I0BM9gUAAABQG4wAAXW3FouePcS3Q610oXg7s6C8vTG38Dxd6D2TW9y+RkEGqI5cm8gFI5yaT7drb0zfQ7k2O5tp2wPl3jPjXG2QcAMAk1LBP6l0gAAAAACMhg4QAAAAgFoISaTAAQAAAKiLcWbnloVJwQAAAABqgxEgAOUYEIwQrfR7me5sJvBgNvfYdNF6o5F+NeXMgvdc2VrpVnD6wLhyz7mxRueg6PnPtZ1cG+vOpuEGuTaba9u9Y9fv/ANYRyo4AkQHCAAAAMBoKvglHh0gAAAAACOp4pUl6AABAAAAGF6oklPgCEEAAAAAUBuMAAEox6CF35nybitTlq47l5qrG26wVovyMTm5c1g0dCLbdjJtLNcWc212mDYPANVg1gABAAAAqJEKfldIBwgAAADAaCrYAWINEAAAAIDaYAQIAAAAwGgqOAJEBwjAVIncgvDc+src+HV23fnqBiNgfSrcTsZoi9m2DQBVFiIEAQAAAEB9VPE7RTpAAAAAAEZTwQ4QIQgAAAAAaoMOEAAAAIDaYAocAAAAgJGwBggAxuTIfJLmPly7mbLcQzPpNLmySk5ixqop3E7GaIvZtg0AVUcKHAAAAIBaCFXy+0PWAAEAAACoDUaAAAAAAIxmPY4A2Z63/Ve2/9r2I7Z/u19+ie37bT9m+zO2ZydfXQAAAADTwjH6rSxFRoAWJb0pIo7bnpH0VdtflPQrkj4cEXfa/pikGyXdNsG6AlhPBi0Iz5Q32pmyTuaxnaKBB8V0M49tVDHupsZy57CobNvJtLFcW8y12WHaPABURgU/wlYcAYqe4/27M/1bSHqTpLv65bskXTuRGgIAAADAKikUgmC7afshSU9LulfS30p6LiLa/U32S7pgwGNvsr3H9p5Dhw6tRp0BAAAATIMY41aSQh2giOhExA9JulDS5ZK+P7fZgMfujIgdEbFj69ato9cUAAAAwNQYZ/3PtK8BeklEPGf7zyW9XtI5tlv9UaALJT01gfoBAAAAmFbr8UKotrdKOtXv/GyQ9KOSPijpy5Kuk3SnpBsk3T3JigJYZwYs/Ha7m5Q1ljIhCEu5x6aD2t1usWCEXJnX6OupXLDCOIv3q6DMMImi5z/XdnJtLNcWc20217b7B8+XA0AVVPAjrMgI0DZJu2w31Zsy99mI+CPbj0q60/bvSPqGpNsnWE8AAAAAGNuKHaCIeFjSZZnyx9VbDwQAAACghqp4dYih1gABAAAAwEvoAAEAAACohZLT3EZFBwiou+6AhdmnaxRKzc8bZpF3u5MUtU5myhbSj6/GyXTReqeT1ruTWfDOh2G95dpEru3k2lhrIW3fuTaba9sDrXYwwlq8zwGgIvibDwAAAGA0jAABAAAAqA06QAAAAADqooprgJjsCwAAAKA2GAEC6qToQuiijx1nwXR3wFdGnXSheCOzoHz2ePr45kK6QH1psZkeYi4TjNBIn1+rmR6jm1ksL0mNVf4KLLe/Qceedqv92gxS9PXpdAuGZWTazmymjeXaYq7N5tq2pMHvhVFN0/scAKYQHSAAAAAAo6ngFDg6QAAAAACGx3WAAAAAANQKHSAAAAAAtUEHCMCqGWchcxUMuNK9O+nzbp44lZTNHpvLlKUfaUsn0oXspzakZa1m5riNtI4ucay/aJjAWoUlrFW4wTgi81p0uumi/lOn0jahTNuZPZYJQTiWhhvk2myubfcrmS+fFpP4LCJYAUCJ6AABAAAAGJrFGiAAAAAAdUIHCAAAAEAtVDQFjkm4AAAAAKaO7TtsP237WwN+b9u/Z3uf7Ydtv7bIfhkBAtbaeg83KGrQwu/M6+PFTAjC8+2kbO6ZdNH6yXPTsvZZ6UdfeyazkL2RC0ZIiiRJubNaVkhAFcIJxjEo5CEfeJCWtTvpSWwvpm2idTRtO3PPpK9tri3m2uzA9/60hyBMAp+DwPox2Y+wT0i6VdInB/z+LZIu7d9eJ+m2/s8zYgQIAAAAwGhijNtKu474iqRnzrDJNZI+GT1fk3SO7W0r7ZcRIAAAAAAjGXPSwRbbe5bd3xkRO4d4/AWSnlx2f3+/7MCZHkQHCAAAAMBoxusAHY6IHWM8PjcnesUaMQUOAAAAQBXtl3TRsvsXSnpqpQcxAgRgdLmFzONe4b2bfnGTW1DeOraYlG04MpOULWZCEF7YnG63NJMJPGhmAhmchiVIUjPzHVRusf56DyhYbYMCD/LbpmXtbtoelxbT86/jadn8M+mxNxxJz3+uLeZDECZw7gkTAFCmgmt5Jmi3pHfZvlO98IPnI+KM098kOkAAAAAARjTJ7/Vsf1rSFeqtFdov6bckzUhSRHxM0j2S3ippn6QFSe8ssl86QAAAAABGM8EOUES8Y4Xfh6RfGHa/dIAAAAAAjKSKM7sJQQAAAABQG4wAAZgukfkqqZMuPG8sLCVl84fTsg1nzydlpzal3/0szmaCEVrpcT2bVq9fyaSk2UifC8EIgxUNPOh089u1O2ngxdJS+meu/UJ6rueOpG1iw8H0vOTaWK4t5tpstm0DQNVV8KONDhAAAACA4ZWfAjcSOkAAAAAAhmblr0Q67VgDBAAAAKA2GAECAAAAMBqmwAGovdyV6RtDDDZnQxDSffpkuvC8dfRkUnbWwcwi+I1pWWc2LVtszqV12ZQWSZKy4QiZ8IbMXIHMK5ZV1bCEouEGkdmum3nKubADSVrMBB4sHk/PYetIut3GzHXDzzrYTh+baWO5tphrs2OHIOTeWwBQsir+aaIDBAAAAGA0dIAAAAAA1EYFO0CEIAAAAACoDUaAAAAAAAwvWAMEAAAAoE7oAAFAxiSS4dppwpoXFpOyucNpYtimuY1JWS4FbqExk5SlR3hxp2lRZJLhWo30tWhmXgpnvlIrmqaWM26C3DjHzsklvnW6aVm7m744S5m0N2lA4tvh9BxuPJAeZ9NTaeLb3OETSVmujeXaIolvAOqCESAAAAAA9VHBDhAhCAAAAABqgxEgAAAAACOp4hS4FUeAbF9k+8u299p+xPa7++Xn2r7X9mP9ny+ffHUBAAAATIUY81aSIiNAbUm/GhEP2t4s6QHb90r6WUn3RcQHbN8s6WZJ751cVQFgmcwicy+dSsoaR9OF7Bta6Xc/0diQOUgaoLDQTRfVS9JiJ7Oo/6y0PrNzaVmrmQtGyJUV+2ux2gEKg+SCDHJy4QadTLhBu5MJPFhMX+/2C/lz0DqS/knLBR6c/UQaWrDh4MmkLNd2cm1s7MADAKiyCn4ErjgCFBEHIuLB/r+PSdor6QJJ10ja1d9sl6RrJ1VJAAAAAFgNQ60Bsr1d0mWS7pd0fkQckHqdJNuvGPCYmyTdJEkXX3zxOHUFAAAAMCWsdboG6EW2N0n6nKT3RMTRoo+LiJ0RsSMidmzdunWUOgIAAACYRut0DZBsz6jX+flURHy+X3zQ9rb+6M82SU9PqpIAAAAApo8ruA5yxQ6QbUu6XdLeiPjQsl/tlnSDpA/0f949kRoCWJ8GXem+UXBgOveB20n3mVu03nx2ISnbmDmEu/PpY5fSYARJWlhKP06Xzkufy4lNaVlrrp2UzcykC/WbubCEzNyDXAhCzqDtioYb5Lbr5Moy4QanTqWvY3sx8yfpeBp4MHck30Y2HkjLNj2Vvra5wINcm8gGHmTa2NghCIPeCwAw7UoeyRlVkRGgN0r6aUnftP1Qv+x96nV8Pmv7RklPSHr7ZKoIAAAAAKtjxQ5QRHxVvTVOOVeubnUAAAAAVEUVQxCGSoEDAAAAgJfQAQIAAABQF4wAAcC4cgvCxwlGaKdhAtZSUtZ8Jn3sxnYmdGBxQ/bQrYX04/TE8XSh/8lzM4EAZ6ePPbUhE4IwVywYodEoFowwTAhCrqzbLRZ40FnMBEecSMtaR9Oy+WfSY2w4mK/3WQfTwIO5wyeSssbRtCwbeJBpO2MFHhB2AGA9qmAHqPB1gAAAAACg6hgBAgAAADC8YAocAAAAgDqhAwQAAACgDixGgABg+hQNRshs13g+LZs/lVkYL6n1wnxSNnd0Nik7cV669HLx3PSjeGlzGgjQ2ZjW59R8WhatzGL7ZuZ1GHSFt9wfs066sdvpc2mcTLebXciUHUvL5jJBFBuOpK/3/OE0xEKSWkdPpnVcWEzLcoEHncxrNk7gAQBgatEBAgAAADCaCn5ZRAcIAAAAwEiYAgcAAACgHkKEIAAAAACoD1fwGs90gABMv27m07UxxnWcc/OVM4vgfTJdbO9MgIIkzSy1ky4mC8IAAA0YSURBVLLmQhqCMPfsXFK29LJcCEL6/JY2pcEB7Y3pdt3ZNEChmxYNvhR27uXOPO1GJougtZC+trPHM2XH0h3OPp++hq1jaYhBYyEfgpA7X7nAi+z5X+057Lk2CwCYCnSAAAAAAIyGKXAAAAAA6oIQBAAAAAD1ECIGGwAAAEB9MAIEAGtlLYIRcmWn8p/0ztSnuXgqKcst4J95biYpm9+QlnXn0ySDdqasO5uGJXRbaZkyRZKy87kb7bSwsZSWtU6moQONTFnzRPraOPN65crUyQdR5IIs1uSbSQIPAKBS6AABAAAAGA0jQAAAAADqwGIKHAAAAIC6iKhkCMIYE+YBAAAAoFoYAQKwfhRdjL7aYQlSfgF+N93WmQX8PpkGI/h4po7NNPBgtpWWRSvzWKeJB5EpkyQXDIRwO/Oc25mAgtxzzr5exV7Dgedgtb+FJNwAAFbEFDgAAAAA9UEHCAAAAEBdMAIEAAAAoB5C+anKU44QBAAAAAC1wQgQgPrJLW4fJxhByi/Az5Xl1tU7LXQ7E1DQaBeqinPhBgMCDwor+vxyigYZlBmlSuABAIymegNAdIAAAAAAjIY1QAAAAADqgwuhAgAAAKgLx+i3Qvu3r7L9Hdv7bN+c+f3P2j5k+6H+7edW2icjQAAAAACmju2mpN+X9GZJ+yV93fbuiHj0tE0/ExHvKrpfRoAAAAAADC/GvK3sckn7IuLxiFiSdKeka8atNiNAACAVTwEbNy0uZ6wEuTHT3VbbtM0FJ90NACbGkjzZz/0LJD257P5+Sa/LbPcTtn9Y0t9I+uWIeDKzzUsYAQIAAAAwmu4YN2mL7T3Lbjedtvfct3yn97j+n6TtEfFqSX8qaddKVWYECAAAAEAZDkfEjjP8fr+ki5bdv1DSU8s3iIgjy+7+L0kfXOmgjAABAAAAGIkjRr4V8HVJl9q+xPaspOsl7f4Xx7e3Lbv7Nkl7V9opI0AAAAAAhlc8zGC03Ue0bb9L0pckNSXdERGP2H6/pD0RsVvSL9l+m6S2pGck/exK+6UDBADDGLSofhLhCEVMW+hAmQg8AIA1FhP/OxQR90i657Sy31z271sk3TLMPukAAQAAABhJ0QuaThPWAAEAAACoDUaAAAAAAIymglOxVxwBsn2H7adtf2tZ2bm277X9WP/nyydbTQAAAABTJSR3R7+VpcgUuE9Iuuq0spsl3RcRl0q6r38fAOqr2y12w3CKvq68tgBQjojRbyVZsQMUEV9RL1JuuWv0z1dZ3SXp2lWuFwAAAACsulFDEM6PiAOS1P/5ikEb2r7J9h7bew4dOjTi4QAAAABMnRjjVpKJp8BFxM6I2BERO7Zu3TrpwwEAAABYI44Y+VaWUTtAB21vk6T+z6dXr0oAAAAAKqGCa4BGjcHeLekGSR/o/7x71WoEAOvZOIv1GxW9dBsBBQCwPoWkCn7EF4nB/rSkv5T0vbb3275RvY7Pm20/JunN/fsAAAAAMNVWHAGKiHcM+NWVq1wXAAAAABVhlbuWZ1SjToEDAAAAUHd0gAAAAADUBh0gAMDE5MIEpi0YgcADAKiP9RqCAAAAAADrBSNAAAAAAEZCCAIAAACA+qADBAAAAKAeopIdINYAAQAAAKgNRoAAAAAADC9UyREgOkAAAAAARlPBGGw6QAAAAABGQgocAAAAgPqgAwQAmJhGBXJrcnXsVnB+BABg3aIDBAAAAGB4IanLCBAAAACAWqjmdYDoAAEAAAAYDR0gAAAAALVBBwgAsCqqEHhQFMEIAIApQgcIAAAAwPAIQQAAAABQHyFF9Ub06QABAAAAGE0F1wCto0nmAAAAAHBmjAABAAAAGB5rgAAAAADUSgWnwNEBAgAAADAaOkAAAAAA6iEq2QEiBAEAAABAbTACBAAAAGB4IanLdYAAAAAA1EUFp8DRAQIAAAAwGjpAAAAAAOohKnkdIEIQAAAAANQGI0AAAAAAhhdSBCEIAAAAAOqiglPg6AABAAAAGE0FQxBYAwQAAACgNhgBAgAAADC8CC6ECgAAAKBGKjgFjg4QAJStUcPZyLnnXMFvEQGg7qKCn910gAAAAACMICo5AlTDrx0BAAAA1BUjQAAAAACGF6rkdYDGGgGyfZXt79jeZ/vm1aoUAAAAgAqI7ui3kow8AmS7Ken3Jb1Z0n5JX7e9OyIeXa3KAQAAAJhOISlqNgJ0uaR9EfF4RCxJulPSNatTLQAAAABTLWLiI0ArzTizPWf7M/3f3297+0r7HKcDdIGkJ5fd398vO71SN9neY3vPoUOHxjgcAAAAgLpYNuPsLZJeJekdtl912mY3Sno2Iv6VpA9L+uBK+x2nA+RMWTIGFhE7I2JHROzYunXrGIcDAAAAME2iGyPfCigy4+waSbv6/75L0pW2c/2Ul4zTAdov6aJl9y+U9NQY+wMAAABQJZOdAldkxtlL20REW9Lzks47007HicH+uqRLbV8i6R8lXS/pP53pAQ888MBh2/8wxjExmi2SDpddCSQ4L9OJ8zKdOC/TadXOi33rauwGPbxfptNK5+W716oiq+WYnv3Sn8ZdW8bYxbztPcvu74yIncvuF5lxVmhW2nIjd4Aiom37XZK+JKkp6Y6IeGSFxzAHrgS290TEjrLrgX+J8zKdOC/TifMynTgv04nzMp3W43mJiKsmfIgiM85e3Ga/7Zakl0l65kw7HetCqBFxj6R7xtkHAAAAAGQUmXG2W9INkv5S0nWS/iwiJjMCBAAAAACTMmjGme33S9oTEbsl3S7pf9vep97Iz/Ur7ZcOUD3sXHkTlIDzMp04L9OJ8zKdOC/TifMynTgvI8jNOIuI31z275OS3j7MPr3CCBEAAAAArBvjxGADAAAAQKXQAVrHbP8329+2/bDtP7R9zrLf3WJ7n+3v2P6xMutZR7av6r/2+2zfXHZ96sr2Rba/bHuv7Udsv7tffq7te20/1v/58rLrWje2m7a/YfuP+vcvsX1//5x8xvZs2XWsI9vn2L6r/7dlr+038H4pl+1f7n9+fcv2p23P834ph+07bD9t+1vLyrLvD/f8Xv//AQ/bfm15Na8fOkDr272SfjAiXi3pbyTdIkm2X6XeArEfkHSVpI/abpZWy5rpv9a/L+ktkl4l6R39c4K115b0qxHx/ZJeL+kX+ufiZkn3RcSlku7r38faerekvcvuf1DSh/vn5FlJN5ZSK/x3SX8cEd8n6TXqnSPeLyWxfYGkX5K0IyJ+UL1F4teL90tZPqHe/6uWG/T+eIukS/u3myTdtkZ1hOgArWsR8Sf9K+JK0tfUy06XpGsk3RkRixHxd5L2Sbq8jDrW1OWS9kXE4xGxJOlO9c4J1lhEHIiIB/v/Pqbef+YuUO987OpvtkvSteXUsJ5sXyjpxyV9vH/fkt4k6a7+JpyTEtg+W9IPq5e4pIhYiojnxPulbC1JG/rXP9ko6YB4v5QiIr6i9Pozg94f10j6ZPR8TdI5tretTU1BB6g+/rOkL/b/fYGkJ5f9bn+/DGuD138K2d4u6TJJ90s6PyIOSL1OkqRXlFezWvqIpF+X1O3fP0/Sc8u+0OE9U47vkXRI0h/0pyd+3PZZ4v1Smoj4R0m/K+kJ9To+z0t6QLxfpsmg9wf/FygRHaCKs/2n/Xm/p9+uWbbNb6g31edTLxZldkUc4Nrh9Z8ytjdJ+pyk90TE0bLrU2e2r5b0dEQ8sLw4synvmbXXkvRaSbdFxGWSXhDT3UrVX09yjaRLJH2XpLPUm1p1Ot4v04fPtRJxHaCKi4gfPdPvbd8g6WpJVy67Ku5+SRct2+xCSU9NpobI4PWfIrZn1Ov8fCoiPt8vPmh7W0Qc6E9JeLq8GtbOGyW9zfZbJc1LOlu9EaFzbLf632rzninHfkn7I+L+/v271OsA8X4pz49K+ruIOCRJtj8v6d+I98s0GfT+4P8CJWIEaB2zfZWk90p6W0QsLPvVbknX256zfYl6C/D+qow61tTXJV3aT+mZVW/B6u6S61RL/bUlt0vaGxEfWvar3ZJu6P/7Bkl3r3Xd6ioibomICyNiu3rvjT+LiJ+U9GVJ1/U345yUICL+SdKTtr+3X3SlpEfF+6VMT0h6ve2N/c+zF88J75fpMej9sVvSz/TT4F4v6fkXp8ph8rgQ6jpme5+kOUlH+kVfi4j/0v/db6i3Lqit3rSfL+b3gknof7v9EfUSe+6IiP9acpVqyfa/lfQXkr6pf15v8j711gF9VtLF6v0H4+0RcfrCVkyY7Ssk/VpEXG37e9QLDDlX0jck/VRELJZZvzqy/UPqhVPMSnpc0jvV+zKV90tJbP+2pP+o3t/zb0j6OfXWkvB+WWO2Py3pCklbJB2U9FuSvqDM+6PfYb1VvdS4BUnvjIg9ZdS7jugAAQAAAKgNpsABAAAAqA06QAAAAABqgw4QAAAAgNqgAwQAAACgNugAAQAAAKgNOkAAAAAAaoMOEAAAAIDaoAMEAAAAoDb+PyOlBq4QuCNmAAAAAElFTkSuQmCC\n",
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
      "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()"
   ]
  }
 ],
 "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",
331
   "version": "3.7.2"
332
333
334
335
336
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}