08_tutorial_shanchen_twocomponent.ipynb 42.8 KB
Newer Older
Michael Kuron's avatar
Michael Kuron committed
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-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"
Michael Kuron's avatar
Michael Kuron committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
   ]
  },
  {
   "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": [
43
44
45
46
47
    "N = 64       # domain size\n",
    "omega_a = 1. # relaxation rate of first component\n",
    "omega_b = 1. # relaxation rate of second component\n",
    "\n",
    "# interaction strength\n",
Michael Kuron's avatar
Michael Kuron committed
48
49
50
    "g_aa = 0.\n",
    "g_ab = g_ba = 6.\n",
    "g_bb = 0.\n",
51
52
53
54
55
    "\n",
    "rho0 = 1.\n",
    "\n",
    "stencil = get_stencil(\"D2Q9\")\n",
    "weights = get_weights(stencil, c_s_sq=sp.Rational(1,3))"
Michael Kuron's avatar
Michael Kuron committed
56
57
58
59
60
61
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
62
63
64
65
66
    "## Data structures\n",
    "\n",
    "We allocate two sets of PDF's, one for each phase. Additionally, for each phase there is one field to store its density and velocity.\n",
    "\n",
    "To run the simulation on GPU, change the `default_target` to gpu"
Michael Kuron's avatar
Michael Kuron committed
67
68
69
70
71
72
73
74
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
75
76
    "dim = len(stencil[0])\n",
    "dh = ps.create_data_handling((N, ) * dim, periodicity=True, default_target='cpu')\n",
Michael Kuron's avatar
Michael Kuron committed
77
    "\n",
78
    "src_a = dh.add_array('src_a', values_per_cell=len(stencil))\n",
Michael Kuron's avatar
Michael Kuron committed
79
80
    "dst_a = dh.add_array_like('dst_a', 'src_a')\n",
    "\n",
81
    "src_b = dh.add_array('src_b', values_per_cell=len(stencil))\n",
Michael Kuron's avatar
Michael Kuron committed
82
83
84
85
86
87
88
89
90
91
92
93
    "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": [
94
95
96
97
    "## Force & combined velocity\n",
    "\n",
    "The two LB methods are coupled using a force term. Its symbolic representation is created in the next cells.\n",
    "The force value is not written to a field, but directly evaluated inside the collision kernel."
Michael Kuron's avatar
Michael Kuron committed
98
99
100
101
102
103
104
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The force between the two components is\n",
105
    "$\\vec{F}_k(\\vec{x})=-\\psi(\\rho_k(\\vec{x}))\\sum\\limits_{k^\\prime\\in\\{A,B\\}}g_{kk^\\prime}\\sum\\limits_{i=1}^{q}w_i\\psi(\\rho_{k^\\prime}(\\vec{x}+\\vec{c}_i))\\vec{c}_i$\n",
Michael Kuron's avatar
Michael Kuron committed
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
    "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",
    "    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",
    "    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": [
174
175
    "collision_a = create_lb_update_rule(stencil=stencil,\n",
    "                                    relaxation_rate=omega_a, \n",
Michael Kuron's avatar
Michael Kuron committed
176
177
178
179
180
181
182
    "                                    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",
183
184
    "collision_b = create_lb_update_rule(stencil=stencil,\n",
    "                                    relaxation_rate=omega_b,\n",
Michael Kuron's avatar
Michael Kuron committed
185
186
187
188
189
    "                                    compressible=True,\n",
    "                                    velocity_input=u_full, density_input=ρ_b,\n",
    "                                    force_model='guo', \n",
    "                                    force=force_b,\n",
    "                                    kernel_type='collide_only',\n",
190
191
192
193
194
195
196
197
198
199
200
201
202
    "                                    optimization={'symbolic_field': src_b})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "stream_a = create_stream_pull_with_output_kernel(collision_a.method, src_a, dst_a, \n",
    "                                                 {'density': ρ_a, 'velocity': u_a})\n",
    "stream_b = create_stream_pull_with_output_kernel(collision_b.method, src_b, dst_b, \n",
    "                                                 {'density': ρ_b, 'velocity': u_b})\n",
Michael Kuron's avatar
Michael Kuron committed
203
    "\n",
204
205
    "opts = {'cpu_openmp': 1,  # number of threads when running on CPU\n",
    "        'target': dh.default_target}\n",
Michael Kuron's avatar
Michael Kuron committed
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
    "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",
221
   "execution_count": 9,
Michael Kuron's avatar
Michael Kuron committed
222
223
224
   "metadata": {},
   "outputs": [],
   "source": [
225
    "init_a = macroscopic_values_setter(collision_a.method, velocity=(0, 0), \n",
Michael Kuron's avatar
Michael Kuron committed
226
    "                                   pdfs=src_a.center_vector, density=ρ_a.center)\n",
227
    "init_b = macroscopic_values_setter(collision_b.method, velocity=(0, 0), \n",
Michael Kuron's avatar
Michael Kuron committed
228
229
230
231
232
233
234
    "                                   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",
235
   "execution_count": 10,
Michael Kuron's avatar
Michael Kuron committed
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
   "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",
262
   "execution_count": 11,
Michael Kuron's avatar
Michael Kuron committed
263
264
265
266
267
268
269
270
271
   "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",
272
    "        sync_ρs()  # collision kernel evaluates force values, that depend on neighboring ρ's\n",
Michael Kuron's avatar
Michael Kuron committed
273
274
275
276
277
278
279
280
281
282
283
284
285
286
    "        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",
287
   "execution_count": 12,
Michael Kuron's avatar
Michael Kuron committed
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
   "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",
312
   "execution_count": 13,
Michael Kuron's avatar
Michael Kuron committed
313
314
315
316
   "metadata": {},
   "outputs": [
    {
     "data": {
317
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA48AAAF0CAYAAACDowz8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfbBkd33f+fdHI4lxBqwHBFglCSMS2RYkSOApGZdS5skImRCk1MKuFOJMKLlUScDBXu8mwruBXWxXmc2WiVMBY8UoklMGQYQxsy4ZoRWwJLHBGoECegBrkFk0HkVjMRJgGUmrme/+0WdM62rudM+9t/t2f+f9qjrVfR66+3fgSh99z+/8fidVhSRJkiRJR3LcZjdAkiRJkrT4LB4lSZIkSRNZPEqSJEmSJrJ4lCRJkiRNZPEoSZIkSZrI4lGSJEmSNJHFoyRJkiRpIotHSZIkSdJEFo9qL8kzklyd5KEk+5L8/Ga3SZKkY5nZLC0ni0cdC34P+BrwA8BlwP+Z5Ac2t0mSJB3TzGZpCVk8qrUkrwOoqndX1WNV9Sngz4AfGva/MMmBJGduZjslSTpWHCmbk/zHJP9lWD6VZOvmtlbSOItHdfd64OOHVpIcB5wEPDBs+hfAfwDOnX/TJEk6Jh0pm38IeFlVXQh8G3jRprRQ0mFZPKq7HwO+Obb+SuDBqvpqkhcB9wM3YfEoSdK8HDabgT8FjquqJ5I8HTgd+JNNaJ+kVVg8qq0kJwDnAG9IsjXJC4H3MeptBPh54N3AXVg8SpI0cxOy+VzgB5J8BrgHeG9VPbxpjZX0FMdvdgOkGToX+DpwB6NbYfYBv1xVNyQ5H7gQ+CCwZVgkSdJsHSmb3wT8RlW9I8lJwOeA3960lkp6CotHdfYi4O6q+pfAv1yx7yrgx6rqIYAkfzzvxkmSdAw6Ujb/TeDW4f0pwLfm2TBJk1k8qrPzgLtXbkzyo8B3DxWOg0eTPLOqvrnyeEmStGEOm82DFwIXJXkLcBD4x3NrlaSpOOZRnb0I+MrKjVV1W1W9ecW2n7BwlI4syVlJPp3k7iR3JnnbYY5Jkn+TZHeSLyV5ydi+HUnuGZYd8229pAVx2GwGqKrXV9WPVtWrqurVVXX7nNsmLZ15Z3OqaqPPQZLUUJLTgdOr6gtJngHcBlxaVXeNHfNa4GeB1zKaUfHXq+rHkpwK7AK2AzV89kdX3AEgSZKOwryz2Z5HSdJUqur+qvrC8P47jG49O2PFYZcAv10jnwNOHoLtNcDNVbV/CKWbgYvn2HxJktqZdzZbPEqSjlqS5wEvBj6/YtcZwH1j63uGbattlyRJG2Ae2TzXCXO2bNtWJ5x86jx/UpIW0v/38H4OPPJINvI7X/OKbfXN/QfW/PnbvvTYncCjY5uurqqrVx43PLz7o8DPVdW3V+4+zFfXEbZrk5143Nb6vuOesdnNkKRN992D3+Hxg4+azUcw1+LxhJNP5ax/+vPz/ElJWkj3ve89G/6dD+4/wOdvOnPNnz/h9K89WlXbj3TM8IDvjwK/U1W/e5hD9gBnja2fCewdtr98xfbPrLmx2jDfd9wz+PGT/t5mN0OSNt0ffetjG/6d3bLZ21YlSVNJEuADjJ7R9murHLYT+IfDzG4vBb5VVfcDNzGagv+UJKcAFw3bJEnSGs07m33OoyS1URyog7P8gQuBnwa+nOTQFPq/CDwXoKreD9zIaDa33cBfAm8e9u1P8kt87wHg76qq/bNsrCRJm69XNls8SlITBRyc4TDCqvrPHH58xPgxBbxllX3XANfMoGmSJC2kbtls8ShJjRxkplc3JUnSUeqUzY55lCRJkiRNZM+jJDVRFAfKp19IkrQoumWzxaMkNTLLcRWSJOnodcpmi0dJaqKAA40CSpKkZdctmy0eJamRTlc3JUnqoFM2O2GOJEmSJGkiex4lqYmCVoPyJUladt2y2eJRkhrp8yQpSZJ66JTNU922muTkJDck+UqSu5P8eJJTk9yc5J7h9ZRZN1aStLqiOLCORcvFbJakxdctm6cd8/jrwCeq6keA84C7gauAW6rqHOCWYV2StFkKDqxj0dIxmyVp0TXL5onFY5LvB34C+ABAVT1eVQ8DlwDXDYddB1w6q0ZKkqTvMZslSZthmp7H5wN/Dvz7JF9M8ltJtgHPqar7AYbXZx/uw0muTLIrya4DjzyyYQ2XJD1ZMRpXsdZFS2XDsvnxenR+rZakY0y3bJ6meDweeAnwG1X1YuARjuI2mKq6uqq2V9X2Ldu2rbGZkqTJwoF1LFoqG5bNJ2brrNooSWqWzdMUj3uAPVX1+WH9BkaB9UCS0wGG132zaaIkaRoFHKy1L1oqZrMkLYFu2TyxeKyq/wbcl+SHh02vAu4CdgI7hm07gI/PpIWSJOlJzGZJ0maY9jmPPwv8TpITgXuBNzMqPD+S5ArgG8AbZ9NESdK0FvEWF82M2SxJS6BTNk9VPFbV7cD2w+x61cY2R5K0VkWvgNKRmc2StPi6ZfO0PY+SpCVwsPoElCRJHXTKZotHSWqi29VNSZKWXbdsnma2VUmSJEnSMc6eR0lqoggHvCYoSdLC6JbNFo+S1EincRWSJHXQKZstHiWpiW7jKiRJWnbdstniUZLaCAeqz60xkiQtv17Z3OdMJEmSJEkzY8+jJDVRwEGvCUqStDC6ZbPFoyQ10mlchSRJHXTKZotHSWqiqte4CkmSll23bO5zJpIkSZKkmbHnUZIaOdjo1hhJkjrolM0Wj5LUxOhZUt5QIknSouiWzRaPktTGbMdVJLkGeB2wr6r+5mH2/8/Am4bV44FzgWdV1f4kXwe+AxwAnqiq7TNrqCRJC6NXNvcpgyXpGHdoOvC1LlO4Frh41d+v+ldVdX5VnQ+8Hfh/qmr/2CGvGPZbOEqSjgndstniUZI0lar6LLB/4oEjlwMfmmFzJEk65s07my0eJamRA5U1L8BpSXaNLVeupQ1J/hqjq6AfHdtcwCeT3LbW75UkaRl1ymbHPEpSE0XWOyj/wQ26pfTvAv9lxW0xF1bV3iTPBm5O8pXhaqkkSW11y2aLR0lq5OBiPIj4MlbcFlNVe4fXfUk+BlwAWDxKktrrlM0LcSaSpPU7NB34WpeNkOQk4GXAx8e2bUvyjEPvgYuAOzbkByVJWmDdstmeR0nSVJJ8CHg5o/EXe4B3AicAVNX7h8P+HvDJqnpk7KPPAT6WBEa588Gq+sS82i1JUlfzzmaLR0lqovirwfWz+f6qy6c45lpG04aPb7sXOG82rZIkaXF1y2aLR0lqZMpnQkmSpDnplM0Wj5LURBUcWIxB+ZIkiX7Z3OdMJEmSJEkzY8+jJLURDjK7cRWSJOlo9cpmi0dJaqLodWuMJEnLrls2WzxKUiMb9UwoSZK0MTpls8WjJDVRhIMznA5ckiQdnW7Z3KcMliRJkiTNjD2PktRIp1tjJEnqoFM2WzxKUhMFHGw0KF+SpGXXLZstHiWpjXCg0XTgkiQtv17ZbPEoSU10u7opSdKy65bNfc5EkiRJkjQz9jxKUiOdbo2RJKmDTtls8ShJTVSl1a0xkiQtu27ZbPEoSY0caBRQkiR10Cmb+5yJJEmSJGlm7HmUpCYKONhoXIUkScuuWzZbPEpSG2l1a4wkScuvVzZPVTwm+TrwHeAA8ERVbU9yKvBh4HnA14H/vqoemk0zJUmTjJ4l1efqpo7MbJakxdctm4+mDH5FVZ1fVduH9auAW6rqHOCWYV2StIkOcNyaFy0ls1mSFlynbF5Piy4BrhveXwdcuv7mSJKkdTCbJUkzM23xWMAnk9yW5Mph23Oq6n6A4fXZh/tgkiuT7Eqy68Ajj6y/xZKkwyrCwVr7oqWzIdn8eD06p+ZK0rGnWzZPO2HOhVW1N8mzgZuTfGXaH6iqq4GrAbaecVatoY2SpCkdXMBbXDQzG5LNJx3/LLNZkmaoUzZPVTxW1d7hdV+SjwEXAA8kOb2q7k9yOrBvhu2UJE1QBQcW8CqlZsNslqTF1y2bJ5bBSbYlecah98BFwB3ATmDHcNgO4OOzaqQkaTqdbo3R6sxmSVoenbJ5mp7H5wAfS3Lo+A9W1SeS3Ap8JMkVwDeAN86umZIkaYzZLEmau4nFY1XdC5x3mO3fBF41i0ZJko7eaFB+n3EVWp3ZLEnLoVs2TzthjiRpCRxg8W5xkSTpWNYpmy0eJamJgoUcHyFJ0rGqWzb36UOVJEmSJM2MPY+S1EavcRWSJC2/Xtnc50wkSRwka14mSXJNkn1J7lhl/8uTfCvJ7cPyjrF9Fyf5apLdSa7awFOWJGmhdcpmex4lqYk5PIj4WuDfAr99hGP+U1W9bnxDki3Ae4FXA3uAW5PsrKq7ZtVQSZIWQbdstniUpEZmeWtMVX02yfPW8NELgN3D4yVIcj1wCWDxKElqr1M2e9uqJOmQ05LsGluuXMN3/HiS/5rkD5K8cNh2BnDf2DF7hm2SJOnIFiqb7XmUpCZGDyJe160xD1bV9nV8/gvAD1bVXyR5LfB7wDlw2EEbtY7fkSRpKXTLZnseJamRWQ7Kn6Sqvl1VfzG8vxE4IclpjK5mnjV26JnA3nX/oCRJS6BTNtvzKElNbPaDiJP8APBAVVWSCxhdoPwm8DBwTpKzgT8DLgP+/qY1VJKkOemWzRaPkqSpJPkQ8HJG4y/2AO8ETgCoqvcDbwD+SZIngO8Cl1VVAU8keStwE7AFuKaq7tyEU5AkqZV5Z7PFoyQ1MuMZ3S6fsP/fMpou/HD7bgRunEW7JElaZJ2y2eJRkrqodQ/KlyRJG6lZNls8SlITBRsyuF6SJG2Mbtls8ShJjXS6uilJUgedstlHdUiSJEmSJrLnUZKa2OzpwCVJ0pN1y2aLR0lqpFNASZLUQadstniUpCaKXjO6SZK07Lpls8WjJDXSaUY3SZI66JTNTpgjSZIkSZrInkdJ6qJ6jauQJGnpNctmi0dJaqLbjG6SJC27btls8ShJjXQKKEmSOuiUzY55lCRJkiRNZM+jJDXRbTpwSZKWXbdstniUpEaqUUBJktRBp2y2eJSkRjo9S0qSpA46ZbPFoyQ1Uc2mA5ckadl1y2YnzJEkSZIkTWTPoyQ10mlchSRJHXTKZotHSWqj14xukiQtv17ZbPEoSY10uropSVIHnbLZ4lGSmih6DcqXJGnZdctmJ8yRJEmSJE1kz6MkdVGjKcElSdKCaJbNFo+S1EinBxFLktRBp2y2eJSkJopeg/IlSVp23bLZMY+SJEmSpInseZSkNno9S0qSpOXXK5stHiWpkU6D8iVJ6qBTNk9922qSLUm+mOT3h/Wzk3w+yT1JPpzkxNk1U5I0jaqsedHyMZslafF1yuajGfP4NuDusfV3A++pqnOAh4ArNrJhkqSjU9UroDQVs1mSFli3bJ6qeExyJvB3gN8a1gO8ErhhOOQ64NJZNFCSJD2V2SxJmrdpex7/NfDPgYPD+jOBh6vqiWF9D3DG4T6Y5Moku5LsOvDII+tqrCTpyA5W1rxMkuSaJPuS3LHK/jcl+dKw/GGS88b2fT3Jl5PcnmTXBp7ysWxDsvnxenT2LZWkY1inbJ5YPCZ5HbCvqm4b33yYQw87FLSqrq6q7VW1fcu2bdO0SZK0RqPbY9a2TOFa4OIj7P9T4GVV9SLgl4CrV+x/RVWdX1Xb13Ju+p6NzOYTs3UmbZQkjXTK5mlmW70QeH2S1wJbge9ndLXz5CTHD1c4zwT2TvODkqTZmeX4iKr6bJLnHWH/H46tfo5RNmg2zGZJWhKdsnliz2NVvb2qzqyq5wGXAZ+qqjcBnwbeMBy2A/j4ehoiSVqfYu0D8odgO+3QrYzDcuU6mnMF8AdPah58Mslt6/xeYTZL0rLols3rec7jvwCuT/LLwBeBD6zjuyRJm+/BjbilNMkrGAXU3x7bfGFV7U3ybODmJF+pqs+u97f0FGazJPWyUNl8VMVjVX0G+Mzw/l7ggqP5vCRptjb7OcRJXsRo9s+fqqpvHtpeVXuH131JPsYoPyweN4DZLEmLrVM2H81zHiVJi2yTnyWV5LnA7wI/XVV/MrZ9W5JnHHoPXAQcdlY4SZJaaZbN67ltVZK0aGZ4eTPJh4CXMxp/sQd4J3ACQFW9H3gHo8dFvG/0yEGeGG61eQ7wsWHb8cAHq+oTs2upJEkLpFE2WzxKkqZSVZdP2P8zwM8cZvu9wHlP/YQkSVqPeWezxaMkNTLL6cAlSdLR65TNFo+S1MiUDxSWJElz0imbLR4lqYmi19VNSZKWXbdstniUpC4KaBRQkiQtvWbZ7KM6JEmSJEkT2fMoSY10GlchSVIHnbLZ4lGSOmkUUJIktdAomy0eJamNtBqUL0nS8uuVzRaPktRJo6ubkiS10CibnTBHkiRJkjSRPY+S1EX1epaUJElLr1k2WzxKUieNbo2RJKmFRtls8ShJrfS5uilJUg99stkxj5IkSZKkiex5lKROGt0aI0lSC42y2eJRkjppFFCSJLXQKJstHiWpiwIazegmSdLSa5bNFo+S1Eg1uropSVIHnbLZCXMkSZIkSRPZ8yhJnTS6uilJUguNstniUZI6aTSuQpKkFhpls8WjJDWSRlc3JUnqoFM2WzxKUhdFq1tjJElaes2y2QlzJEmSJEkT2fMoSW2k1bgKSZKWX69stniUpE4a3RojSVILjbLZ4lGSOmkUUJIktdAomx3zKEmSJEmayJ5HSeqk0dVNSZJaaJTNFo+S1EXRalC+JElLr1k2WzxKUiOdHkQsSVIHnbLZMY+S1EmtY5kgyTVJ9iW5Y5X9SfJvkuxO8qUkLxnbtyPJPcOyYz2nKEnSUmmUzRaPkqRpXQtcfIT9PwWcMyxXAr8BkORU4J3AjwEXAO9McspMWypJ0rHhWuaYzRaPkqSpVNVngf1HOOQS4Ldr5HPAyUlOB14D3FxV+6vqIeBmjhx0kiRpCvPOZsc8SlIjmzyu4gzgvrH1PcO21bZLktRep2yea/H4tAce5a//2lfm+ZOStJD2fevR2Xzx+mZ0Oy3JrrH1q6vq6qP4/OF+vI6wXQvgseds5Wv/9Ec2uxmStOkee9/W2Xxxo2y251GSuphycP0RPFhV29fx+T3AWWPrZwJ7h+0vX7H9M+v4HUmSlkOzbHbMoyRpo+wE/uEws9tLgW9V1f3ATcBFSU4ZBuNfNGyTJEmztaHZbM+jJHUyw5tBk3yI0VXK05LsYTRL2wkAVfV+4EbgtcBu4C+BNw/79if5JeDW4aveVVVHGtwvSVIfjbJ5YvGYZCvwWeBpw/E3VNU7k5wNXA+cCnwB+Omqenz6U5UkbbRZDsqvqssn7C/gLavsuwa4ZhbtOhaZzZK0PDpl8zS3rT4GvLKqzgPOBy4eujzfDbynqs4BHgKuOJofliTNwAwfRKyFYjZL0rJolM0Ti8fhmSB/MayeMCwFvBK4Ydh+HXDpTFooSZKexGyWJG2GqSbMSbIlye3APkYPkPwa8HBVPTEcsupzQZJcmWRXkl2P14ymppckjTS6uqkj26hsPvDII/NpsCQdqxpl81TFY1UdqKrzGU3hegFw7uEOW+WzV1fV9qrafmJm9OwUSRKp9S1aLhuVzVu2bZtlMyXpmNYtm49qttWqejjJZ4CXAicnOX64wnnoeSGSpM20vgcRawmZzZK04Bpl88SexyTPSnLy8P77gJ8E7gY+DbxhOGwH8PFZNVKSNKVGt8ZodWazJC2RRtk8Tc/j6cB1SbYwKjY/UlW/n+Qu4Pokvwx8EfjADNspSZK+x2yWJM3dxOKxqr4EvPgw2+9lNMZCkrQgFnF8hDae2SxJy6NTNh/VmEdJ0oJrFFCSJLXQKJstHiWpiwWdmU2SpGNWs2ye6lEdkiRJkqRjmz2PktRJo6ubkiS10CibLR4lqZNGASVJUguNstniUZIa6TSuQpKkDjpls2MeJUmSJEkTWTxKkiRJkibytlVJ6qTRrTGSJLXQKJstHiWpi2bPkpIkaek1y2aLR0nqpFFASZLUQqNstniUpE4aBZQkSS00ymYnzJEkSZIkTWTPoyQ1EXqNq5Akadl1y2aLR0nqpFFASZLUQqNstniUpC6azegmSdLSa5bNjnmUJEmSJE1kz6MkddLo6qYkSS00ymaLR0nqpFFASZLUQqNstniUpEY6jauQJKmDTtls8ShJnTQKKEmSWmiUzU6YI0mSJEmayJ5HSeqiaHV1U5Kkpdcsmy0eJamRTuMqJEnqoFM2e9uqJHVS61imkOTiJF9NsjvJVYfZ/54ktw/LnyR5eGzfgbF9O9dzmpIkLY1G2WzPoyQ1Msurm0m2AO8FXg3sAW5NsrOq7jp0TFX9/NjxPwu8eOwrvltV58+uhZIkLZ5O2WzPoyRpWhcAu6vq3qp6HLgeuOQIx18OfGguLZMk6dg012y2eJSkTtZ3a8xpSXaNLVeu+PYzgPvG1vcM254iyQ8CZwOfGtu8dfjezyW5dD2nKUnS0miUzd62KkldrH9GtweravsR9meVXz2cy4AbqurA2LbnVtXeJM8HPpXky1X1tbU2VpKkhdcsm+15lKQmss5lCnuAs8bWzwT2rnLsZay4Laaq9g6v9wKf4cljLiRJaqdbNls8SpKmdStwTpKzk5zIKISeMjNbkh8GTgH+aGzbKUmeNrw/DbgQuGvlZyVJ0lGZazZ726okdTLDGd2q6okkbwVuArYA11TVnUneBeyqqkNhdTlwfVWNt+Zc4DeTHGR04fJXx2eCkySprUbZbPEoSY3M+kHEVXUjcOOKbe9Ysf6/HeZzfwj8rZk2TpKkBdQpmy0eJamTGQeUJEk6So2y2eJRkjppFFCSJLXQKJudMEeSJEmSNJE9j5LURc1+XIUkSToKzbLZ4lGSOmkUUJIktdAomy0eJamRTlc3JUnqoFM2WzxKUieNAkqSpBYaZbMT5kiSJEmSJrLnUZIa6XRrjCRJHXTK5ok9j0nOSvLpJHcnuTPJ24btpya5Ock9w+sps2+uJGlVtc5FS8NslqQl0Sybp7lt9QngF6rqXOClwFuSvAC4Crilqs4BbhnWJUmbqVFA6YjMZklaFo2yeWLxWFX3V9UXhvffAe4GzgAuAa4bDrsOuHRWjZQkSd9jNkuSNsNRjXlM8jzgxcDngedU1f0wCrEkz17lM1cCVwJsPe7p62mrJOkIQq9xFZrOerP5+JO8s1WSZqVbNk8922qSpwMfBX6uqr497eeq6uqq2l5V20/M1rW0UZI0rUa3xmiyjcjmLdu2za6BkqRW2TxVz2OSExiF0+9U1e8Omx9IcvpwZfN0YN+sGilJmk5qAZNGM2E2S9Jy6JTN08y2GuADwN1V9Wtju3YCO4b3O4CPb3zzJElTazajm1ZnNkvSkmiWzdP0PF4I/DTw5SS3D9t+EfhV4CNJrgC+AbxxNk2UJEkrmM2SpLmbWDxW1X9mNNbzcF61sc2RJK1Hp0H5Wp3ZLEnLo1M2H9Vsq5KkBdcooCRJaqFRNls8SlIjna5uSpLUQadstniUpE4aBZQkSS00yuapn/MoSZIkSTp22fMoSV1Ur1tjJElaes2y2eJRkjppFFCSJLXQKJstHiWpidDr6qYkScuuWzY75lGSJEmSNJE9j5LUSTW6vClJUgeNstniUZIa6XRrjCRJHXTKZotHSeqiaDUoX5Kkpdcsmy0eJamRHNzsFkiSpHGdstkJcyRJkiRJE9nzKEmdNLo1RpKkFhplsz2PktRIau3LVN+fXJzkq0l2J7nqMPv/UZI/T3L7sPzM2L4dSe4Zlh0bd9aSJC2uTtlsz6MkdVHMdDrwJFuA9wKvBvYAtybZWVV3rTj0w1X11hWfPRV4J7B9aOltw2cfmlmDJUnabM2y2Z5HSWpkxlc3LwB2V9W9VfU4cD1wyZRNew1wc1XtH0LpZuDitZyjJEnLpFM2WzxKkg45LcmuseXKFfvPAO4bW98zbFvpv0vypSQ3JDnrKD8rSZK+Z6Gy2dtWJamT9d0Z82BVbT/C/kzxi/8X8KGqeizJPwauA1455WclSeqnUTbb8yhJTYSZ3xqzBzhrbP1MYO/4AVX1zap6bFj9d8CPTvtZSZK66ZbNFo+S1EXV+pbJbgXOSXJ2khOBy4Cd4wckOX1s9fXA3cP7m4CLkpyS5BTgomGbJEl9Nctmb1uVJE2lqp5I8lZGwbIFuKaq7kzyLmBXVe0E/lmS1wNPAPuBfzR8dn+SX2IUcgDvqqr9cz8JSZIamXc2WzxKUiPTPhNqrarqRuDGFdveMfb+7cDbV/nsNcA1M22gJEkLplM2WzxKUidOQSNJ0mJplM0Wj5LUyKyvbkqSpKPTKZstHiWpiwIONkooSZKWXbNsdrZVSZIkSdJE9jxKUid9Lm5KktRDo2y2eJSkRjqNq5AkqYNO2WzxKEmdTPdAYUmSNC+NstniUZIa6XR1U5KkDjplsxPmSJIkSZImsudRkrooWg3KlyRp6TXLZotHSWoiQBqNq5Akadl1y2aLR0nq5OBmN0CSJD1Jo2x2zKMkSZIkaSJ7HiWpkU63xkiS1EGnbLZ4lKQumg3KlyRp6TXLZotHSWqjWj2IWJKk5dcrmy0eJamRTg8iliSpg07Z7IQ5kiRJkqSJ7HmUpE4a3RojSVILjbJ5Ys9jkmuS7Etyx9i2U5PcnOSe4fWU2TZTkjRRQQ6ufdHyMJslaUk0y+Zpblu9Frh4xbargFuq6hzglmFdkrTZqta+aJlci9ksScuhUTZPLB6r6rPA/hWbLwGuG95fB1y6we2SJEmrMJslSZthrWMen1NV9wNU1f1Jnr3agUmuBK4E2Hrc09f4c5KkqSzeRUrNz5qy+fiTvLtVkmaqUTbPfMKcqroauBrgpOOf1eh/OklaPFnAW1y0eMazeesZZ/lHI0kz1Cmb1/qojgeSnA4wvO7buCZJktas0bgKHTWzWZIWUaNsXmvxuBPYMbzfAXx8Y5ojSVqzAg6uY9GyM5sladE0y+ZpHtXxIeCPgB9OsifJFcCvAq9Ocg/w6mFdkiTNgdksSdoME8c8VtXlq+x61Qa3RZK0DqFajavQ6sxmSVoO3bJ55hPmSJLmqFFASZLUQqNstniUpE4aBZQkSS00ymaLR0nq4tCgfEmStBiaZfNaZ1uVJA3eBvsAAAqhSURBVEmSJB1DLB4lqZFUrXmZ6vuTi5N8NcnuJFcdZv//mOSuJF9KckuSHxzbdyDJ7cOycwNPW5KkhdUpm71tVZI6meG4iiRbgPcyegzEHuDWJDur6q6xw74IbK+qv0zyT4D/A/gfhn3frarzZ9ZASZIWUaNstudRktqoUUCtdZnsAmB3Vd1bVY8D1wOXPKkFVZ+uqr8cVj8HnLmhpyhJ0lLplc0Wj5KkQ05LsmtsuXLF/jOA+8bW9wzbVnMF8Adj61uH7/1ckks3qM2SJHW2UNnsbauS1EWx3ltjHqyq7UfYn1V+9akHJv8A2A68bGzzc6tqb5LnA59K8uWq+tramytJ0oJrls0Wj5LUyWynA98DnDW2fiawd+VBSX4S+F+Al1XVY4e2V9Xe4fXeJJ8BXgxYPEqSemuUzd62KkmNzHhGt1uBc5KcneRE4DLgSTOzJXkx8JvA66tq39j2U5I8bXh/GnAhMD6YX5Kkljplsz2PktTJDGd0q6onkrwVuAnYAlxTVXcmeRewq6p2Av8KeDrwH5MAfKOqXg+cC/xmkoOMLlz+6oqZ4CRJ6qlRNls8SpKmVlU3Ajeu2PaOsfc/ucrn/hD4W7NtnSRJx555ZrPFoyR1UcDB2V3dlCRJR6lZNls8SlIbUz8TSpIkzUWvbLZ4lKROGgWUJEktNMpmi0dJ6qRRQEmS1EKjbPZRHZIkSZKkiex5lKQumg3KlyRp6TXLZotHSWqjoA5udiMkSdJf6ZXNFo+S1EmjcRWSJLXQKJsd8yhJkiRJmsieR0nqotm4CkmSll6zbLZ4lKROGt0aI0lSC42y2eJRkjppFFCSJLXQKJstHiWpjWoVUJIkLb9e2eyEOZIkSZKkiex5lKQuCjjY51lSkiQtvWbZbPEoSZ00ujVGkqQWGmWzxaMkddIooCRJaqFRNls8SlIb1epZUpIkLb9e2eyEOZIkSZKkiex5lKQuCqr6DMqXJGnpNctmi0dJ6qTRrTGSJLXQKJstHiWpk0aD8iVJaqFRNjvmUZIkSZI0kT2PktRFVasHEUuStPSaZbPFoyR10ujWGEmSWmiUzRaPktRINbq6KUlSB52y2eJRktqoVlc3JUlafr2y2QlzJEmSJEkT2fMoSV0UrZ4lJUnS0muWzevqeUxycZKvJtmd5KqNapQkaY3q4NoXtWA2S9KCaZTNa+55TLIFeC/wamAPcGuSnVV110Y1TpI0vQKq0dVNHT2zWZIWS7dsXk/P4wXA7qq6t6oeB64HLtmYZkmSjlrVzK9uTurVSvK0JB8e9n8+yfPG9r192P7VJK/ZsPPWOLNZkhZJs2xeT/F4BnDf2PqeYdvKxl6ZZFeSXY/Xo+v4OUnSZhrr1fop4AXA5UlesOKwK4CHqupvAO8B3j189gXAZcALgYuB9w3fp4111Nl84JFH5tY4SdLGmnc2r6d4zGG2PaVPtqqurqrtVbX9xGxdx89Jkiapg7XmZQrT9GpdAlw3vL8BeFWSDNuvr6rHqupPgd3D92ljHXU2b9m2bQ7NkqRjV6dsXk/xuAc4a2z9TGDvOr5PkrRes701Zpperb86pqqeAL4FPHPKz2r9zGZJWjSNsnk9j+q4FTgnydnAnzHq8vz7R/rAtw88+OBN+//d/wucBjy4jt9eJJ7LYupyLl3OAzyXlX5wIxoy7js8dNP/XTecto6v2Jpk19j61VV19dj6NL1aqx0zVY+Y1u2os/mxvXse3P2//oLZvLg8l8XT5TzAc1nJbJ6QzWsuHqvqiSRvBW4CtgDXVNWdEz7zLIAku6pq+1p/e5F4Loupy7l0OQ/wXOahqi6e8U9M06t16Jg9SY4HTgL2T/lZrZPZPOK5LKYu59LlPMBzmYdu2byu5zxW1Y1V9UNV9der6lfW812SpIX3V71aSU5k1Ku1c8UxO4Edw/s3AJ+qqhq2XzbM+HY2cA7wx3Nq9zHFbJakY8pcs3k9t61Kko4hq/VqJXkXsKuqdgIfAP5Dkt2MrmpeNnz2ziQfAe4CngDeUlUHNuVEJElqYt7ZvFnF49WTD1kansti6nIuXc4DPJcWqupG4MYV294x9v5R4I2rfPZXAHvCFlenv2vPZTF1OZcu5wGeSwvzzOaMeiwlSZIkSVrdusY8SpIkSZKODXMvHpNcnOSrSXYnuWrev78eSa5Jsi/JHWPbTk1yc5J7htdTNrON00hyVpJPJ7k7yZ1J3jZsX8Zz2Zrkj5P81+Fc/vdh+9lJPj+cy4eHAcQLL8mWJF9M8vvD+rKex9eTfDnJ7Yeml17Gvy+AJCcnuSHJV4Z/Zn58Wc9FWo3ZvPnM5sVlNi8es3nzzLV4TLIFeC/wU8ALgMuTvGCebVina4GV0+1eBdxSVecAtwzri+4J4Beq6lzgpcBbhv8flvFcHgNeWVXnAecDFyd5KfBu4D3DuTwEXLGJbTwabwPuHltf1vMAeEVVnT82bfYy/n0B/Drwiar6EeA8Rv//LOu5SE9hNi8Ms3lxmc2Lx2zeJPPuebwA2F1V91bV48D1wCVzbsOaVdVnGc1QNO4S4Lrh/XXApXNt1BpU1f1V9YXh/XcY/QN3Bst5LlVVfzGsnjAsBbwSuGHYvhTnkuRM4O8AvzWshyU8jyNYur+vJN8P/ASjWcqoqser6mGW8FykIzCbF4DZvJjM5sVjNm+ueRePZwD3ja3vGbYts+dU1f0w+hc/8OxNbs9RSfI84MXA51nScxluJ7kd2AfcDHwNeLiqnhgOWZa/s38N/HPg4LD+TJbzPGD0HwmfTHJbkiuHbcv49/V84M+Bfz/csvRbSbaxnOcircZsXjBm80IxmxeP2byJ5l085jDbnO51kyR5OvBR4Oeq6tub3Z61qqoDVXU+cCajK+jnHu6w+bbq6CR5HbCvqm4b33yYQxf6PMZcWFUvYXQb3FuS/MRmN2iNjgdeAvxGVb0YeARvg1E/y/zvmnbM5sVhNi8ss3kTzbt43AOcNbZ+JrB3zm3YaA8kOR1geN23ye2ZSpITGIXT71TV7w6bl/JcDhluWfgMo7EiJyc59BzTZfg7uxB4fZKvM7pl7JWMrnYu23kAUFV7h9d9wMcY/YfDMv597QH2VNXnh/UbGAXWMp6LtBqzeUGYzQvHbF5MZvMmmnfxeCtwzjBL1YnAZcDOObdho+0EdgzvdwAf38S2TGW4X/8DwN1V9Wtju5bxXJ6V5OTh/fcBP8lonMingTcMhy38uVTV26vqzKp6HqN/Lj5VVW9iyc4DIMm2JM849B64CLiDJfz7qqr/BtyX5IeHTa8C7mIJz0U6ArN5AZjNi8dsXkxm8+ZK1Xx72pO8ltFVmy3ANVX1K3NtwDok+RDwcuA04AHgncDvAR8Bngt8A3hjVa0cuL9Qkvxt4D8BX+Z79/D/IqOxFct2Li9iNCh6C6OLIR+pqncleT6jq4SnAl8E/kFVPbZ5LZ1ekpcD/1NVvW4Zz2No88eG1eOBD1bVryR5Jkv29wWQ5HxGEyWcCNwLvJnhb40lOxdpNWbz5jObF5vZvFjM5s0z9+JRkiRJkrR85n3bqiRJkiRpCVk8SpIkSZImsniUJEmSJE1k8ShJkiRJmsjiUZIkSZI0kcWjJEmSJGkii0dJkiRJ0kQWj5IkSZKkif5/tWv+V9+zaRMAAAAASUVORK5CYII=\n",
Michael Kuron's avatar
Michael Kuron committed
318
      "text/plain": [
319
       "<Figure size 1152x432 with 4 Axes>"
Michael Kuron's avatar
Michael Kuron committed
320
321
      ]
     },
322
323
324
     "metadata": {
      "needs_background": "light"
     },
Michael Kuron's avatar
Michael Kuron committed
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
     "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",
352
   "execution_count": 14,
Michael Kuron's avatar
Michael Kuron committed
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
   "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",
373
   "execution_count": 15,
Michael Kuron's avatar
Michael Kuron committed
374
375
376
377
   "metadata": {},
   "outputs": [
    {
     "data": {
378
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA48AAAF0CAYAAACDowz8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df7DldX3n+ee7bzd0AioggbA0BpwlGUwGwXShU0z5MyI6rjA1OgvjZDouKWpmNWsyv4LZGt0lSRXubMVkKkbToz3glIIOiaE3RUQWZZ1MRkKjjPyS0BJHbpqhhQYVpMF773v/ON+Oh8u995x7z/l+zznv+3xUnerz/Z7vOefzbS731e/v58c3MhNJkiRJktayZdINkCRJkiRNP4tHSZIkSdJAFo+SJEmSpIEsHiVJkiRJA1k8SpIkSZIGsniUJEmSJA1k8ShJkiRJGsjiUZIkSZI0kMWjyouIF0TE7oh4PCIORsSvTLpNkiRtZmazNJssHrUZ/BHwDeDHgUuA/zsifnyyTZIkaVMzm6UZZPGo0iLirQCZ+cHMfCYzvwD8FfCTzes/HRGLEbFjku2UJGmzWCubI+I/RsR/bh5fiIjtk22tpH4Wj6rubcANRzYiYgvwIuCRZtevAv8BOKv7pkmStCmtlc0/CbwmM88HvgucPZEWSlqRxaOqeyXwWN/264FHM/P+iDgbeBi4CYtHSZK6smI2A38JbMnMhYg4FjgF+IsJtE/SKiweVVZEbAPOBN4eEdsj4qeB36PX2wjwK8AHgXuxeJQkqXUDsvks4Mcj4lbgAeDDmfnExBor6Xm2TroBUovOAr4J3E1vKMxB4Dcy8/qIOAc4H/gUMNc8JElSu9bK5ncCH8nM90fEi4AvA5+YWEslPY/Foyo7G7gvM/818K+XvXYF8MrMfBwgIv6868ZJkrQJrZXNPwPc3jw/HvhOlw2TNJjFoyp7OXDf8p0R8bPA00cKx8bhiHhxZj62/HhJkjQ2K2Zz46eBCyLi3cAS8E86a5WkoTjnUZWdDXx9+c7MvCMz37Vs36stHKW1RcRpEfHFiLgvIu6JiPeucExExL+NiP0R8bWIeEXfa7si4oHmsavb1kuaEitmM0Bmvi0zfzYz35CZb8zMOztumzRzus7myMxxn4MkqaCIOAU4JTO/EhEvAO4ALs7Me/uOeQvwS8Bb6K2o+DuZ+cqIOAHYB+wEsnnvzy4bASBJktah62y251GSNJTMfDgzv9I8/x69oWenLjvsIuAT2fNl4Lgm2N4E3JyZh5pQuhm4sMPmS5JUTtfZbPEoSVq3iDgdOBe4bdlLpwIP9W3PN/tW2y9Jksagi2zudMGcuWOPya0vPqHLr5SkqbTw2CEWn3wqxvmZb3rdMfnYocUNv/+Orz1zD3C4b9fuzNy9/Ljm5t1/APxyZn53+csrfHSusV8TdlRsz+1xzKSbIUkTdzif4tk8bDavodPiceuLT+CUK543h1OSNp2Hr/qdsX/mo4cWue2mHRt+/7ZTvnE4M3eudUxzg+8/AD6ZmX+4wiHzwGl92zuAA83+1y7bf+uGG6ux2R7H8Kqj3zzpZkjSxH35mT8Z+2dWy2aHrUqShhIRAXyc3j3afmuVw/YC/7hZ2e1VwHcy82HgJnpL8B8fEccDFzT7JEnSBnWdzd7nUZLKSBZzqc0vOB/4eeCuiDiyhP6vAS8ByMyPAjfSW81tP/B94F3Na4ci4tf54Q3Ar8zMQ202VpKkyauVzd0Xj85wkaRWJLDU4i/ZzPxTVp4f0X9MAu9e5bU9wJ4WmiZJ0lSqls32PEpSIUu0enVTkiStU6Vsds6jJEmSJGkgex4lqYgkWUznBkiSNC2qZbPFoyQV0ua8CkmStH6Vsrnb4rE3Y1SS1EKOJLBYKKDUjQB6K71L0ubWxm/Catlsz6MkFVLp6qYkSRVUymYXzJEkSZIkDWTPoyQVkVBqUr4kSbOuWjZbPEpSIU4rlyRpulTK5qGKx4g4DvgY8DP0Cuj/Bbgf+DRwOvBN4B9k5uMDP2vRSfmS1IYkS03K19rGmc2SpHZUy+Zh5zz+DvC5zPybwMuB+4ArgFsy80zglmZbkjQpCYsjPDRzzGZJmnbFsnlg8RgRLwReDXwcIDOfzcwngIuAa5rDrgEubquRkiTph8xmSdIkDNPz+FLg28C/j4ivRsTHIuIY4OTMfBig+fOkld4cEZdHxL6I2Lf45FNja7gk6bmO3Ep3ow/NlLFl87M8012rJWmTqZbNwxSPW4FXAB/JzHOBp1jHMJjM3J2ZOzNz59yxx2ywmZKkwYLFER6aKWPL5qM4uq02SpKKZfMwC+bMA/OZeVuzfT29gHokIk7JzIcj4hTg4MBPSohpLKElqWstzGNIYGkK50eoFePL5gjY4m2fJYkYf7FWLZsHpkVm/nfgoYj4qWbXG4B7gb3ArmbfLuCGVlooSZKew2yWJE3CsPd5/CXgkxFxFPAg8C56hednIuIy4FvAO9ppoiRpWNM4xEWtMZslaQZUyuahisfMvBPYucJLbxhvcyRJG5XUCiitzWyWpOlXLZuH7XmUJM2ApawTUJIkVVApmy0eJamIalc3JUmaddWyufvicanOX54kSZIkbRb2PEpSEUmwONTteyVJUheqZbPFoyQVUmlehSRJFVTKZotHSSqi2rwKSZJmXbVstniUpDKCxawzNEaSpNlXK5u7LR4TYrHTb5Sk6ZSTboAkSdL62PMoSUUksFRoUr4kSbOuWjZbPEpSIZXmVUiSVEGlbLZ4lKQiMmvNq5AkadZVy+Y6ZyJJkiRJak2nPY8BhItESFJrA1iWCg2NUYfCnxtJakulbHbYqiQV0buXlANKJEmaFtWy2eJRkspod15FROwB3goczMyfWeH1fwm8s9ncCpwF/FhmHoqIbwLfAxaBhczc2VpDJUmaGrWyuU4ZLEmb3JHlwDf6GMLVwIWrfn/mv8nMczLzHOB9wP+XmYf6Dnld87qFoyRpU6iWzRaPkqShZOaXgEMDD+y5FLi2xeZIkrTpdZ3NnQ9bjaWuv1GSNo/FHGlS/okRsa9ve3dm7l7vh0TEj9K7Cvqevt0JfD4iEvj9jXyuJEmzqFI2O+dRkopIYtRJ+Y+OaUjp/wT852XDYs7PzAMRcRJwc0R8vblaKklSWdWy2eJRkgpZmo4bEV/CsmExmXmg+fNgRHwWOA+weJQklVcpm6fiTCRJozuyHPhGH+MQES8CXgPc0LfvmIh4wZHnwAXA3WP5QkmSpli1bLbnUZI0lIi4FngtvfkX88AHgG0AmfnR5rC/B3w+M5/qe+vJwGejdyP6rcCnMvNzXbVbkqSqus7mbovHbB6StNm18LswiVEn5a/9+ZmXDnHM1fSWDe/f9yDw8nZapbHY4kAkSWpDtWy251GSChnynlCSJKkjlbLZ4lGSisiExemYlC9JkqiXzXXORJIkSZLUGnseJamMYIn25lVIkqT1qpXN3RePLpgjSa1Iag2NkSRp1lXLZnseJamQcd0TSpIkjUelbLZ4lKQikmCpxeXAJUnS+lTL5jplsCRJkiSpNfY8SlIhlYbGSJJUQaVs7rx4DBfMkaRWJLBUaFK+JEmzrlo22/MoSWUEi4WWA5ckafbVymaLR0kqotrVTUmSZl21bK5zJpIkSZKk1tjzKEmFVBoaI0lSBZWy2eJRkorIjFJDYyRJmnXVstniUZIKWSwUUJIkVVApm+uciSRJkiSpNfY8SlIRCSwVmlchSdKsq5bNFo+SVEaUGhojSdLsq5XNQxWPEfFN4HvAIrCQmTsj4gTg08DpwDeBf5CZjw/6rKxTeEvSVOndS8pfspvFOLNZktSOatm8njL4dZl5TmbubLavAG7JzDOBW5ptSdIELbJlww/NJLNZkqZcpWwepUUXAdc0z68BLh69OZIkaQRmsySpNcMWjwl8PiLuiIjLm30nZ+bDAM2fJ630xoi4PCL2RcS+xaeeGr3FkqQVJcFSbvyhmTOWbH42D3fUXEnafKpl87AL5pyfmQci4iTg5oj4+rBfkJm7gd0A2089LTfQRknSkJamcIiLWjOWbH7R3IlmsyS1qFI2D1U8ZuaB5s+DEfFZ4DzgkYg4JTMfjohTgIMttlOSNEAmLE7hVUq1w2yWpOlXLZsHFo8RcQywJTO/1zy/ALgS2AvsAq5q/rxhqG+s83cnSVNnGoe4aPzGns2SpNZUyuZheh5PBj4bEUeO/1Rmfi4ibgc+ExGXAd8C3tFeMyVJUh+zWZLUuYHFY2Y+CLx8hf2PAW9oo1GSpPXrTcqvM69CqzObJWk2VMvmYRfMkSTNgEXnBkiSNFUqZbPFoyQVkdSaVyFJ0qyrls3dFo+BC+ZIEvi7UNNlaWnSLZAkzQB7HiWpjFrzKiRJmn21srnOmUiSWCI2/BgkIvZExMGIuHuV118bEd+JiDubx/v7XrswIu6PiP0RccUYT1mSpKlWKZvteZSkIjq4EfHVwO8Cn1jjmP+UmW/t3xERc8CHgTcC88DtEbE3M+9tq6GSJE2Datls8ShJhbQ5NCYzvxQRp2/grecB+5vbSxAR1wEXARaPkqTyKmVz58VjoSG/klTNiRGxr297d2buXudn/O2I+K/AAeBfZOY9wKnAQ33HzAOvHK2pkiRtClOVzfY8SlIRvRsRjzQ05tHM3DnC+78C/ERmPhkRbwH+CDiTldeWzRG+R5KkmVAtm+0HlKRC2pyUP0hmfjczn2ye3whsi4gT6V3NPK3v0B30rn5KklRepWy251GSipj0jYgj4seBRzIzI+I8ehcoHwOeAM6MiDOAvwIuAf7hxBoqSVJHqmWzxaMkaSgRcS3wWnrzL+aBDwDbADLzo8DbgX8aEQvA08AlmZnAQkS8B7gJmAP2NPMtJEnSCLrO5k6LxwQmWHhL0tRoa8Jfyyu6XTrg9d+lt1z4Sq/dCNzYRrskSZpmlbLZnkdJqiJHnpQvSZLGqVg2WzxKUhEJY5lcL0mSxqNaNls8SlIhla5uSpJUQaVs9lYdkiRJkqSBuu15DCxXJQlWvjXviCa9HLhmWLa1hJMkbW7Vstlhq5JUSKWAkiSpgkrZbPEoSUUktVZ0kyRp1lXLZotHSSqk0opukiRVUCmbnYEoSZIkSRqo857H3OKkfElqRdaaVyFJ0swrls0OW5WkIqqt6CZJ0qyrls0Wj5JUSKWAkiSpgkrZ7JxHSZIkSdJA9jxKUhHVlgOXJGnWVcvmbovHgJzr9BslaTq1lCNZKKDUkUxYWpp0KyRp8rKdhT0rZbM9j5JUSKV7SUmSVEGlbLZ4lKQisthy4JIkzbpq2eyCOZIkSZKkgex5lKRCKs2rkCSpgkrZ3HnxmPZ1SlJLaq3oJknS7KuVzfY8SlIhla5uSpJUQaVstniUpCKSWpPyJUmaddWy2UGkkiRJkqSB7HmUpCqytfsbS5KkjSiWzd0Xj1sK/e1J0pSpdCNidSOBrPQvG0naoLZ+E1bKZnseJamIXhFQJ6AkSZp11bLZOY+SJEmSpIHseZSkMmrdS0qSpNlXK5stHiWpEKeuSZI0XSpl89DDViNiLiK+GhF/3GyfERG3RcQDEfHpiDiqvWZKkoaRGRt+aPaYzZI0/Spl83p6Ht8L3Ae8sNn+IPChzLwuIj4KXAZ8ZM1PCJxlKUlAGwuvZdaalK+hjJ7NkqTWVMvmoUq5iNgB/F3gY812AK8Hrm8OuQa4uI0GSpKk5zObJUldG7Yf8LeBfwUsNdsvBp7IzIVmex44daU3RsTlEbEvIvYtPvnkSI2VJK1tKWPDj0EiYk9EHIyIu1d5/Z0R8bXm8WcR8fK+174ZEXdFxJ0RsW+Mp7yZjSWbf5CH22+pJG1ilbJ5YPEYEW8FDmbmHf27Vzh0xamgmbk7M3dm5s65Y48dpk2SpA3qDY/Z2GMIVwMXrvH6XwKvycyzgV8Hdi97/XWZeU5m7tzIuemHxpnN22J7K22UJPVUyuZh5jyeD7wtIt4CbKc3r+K3geMiYmtzhXMHcGCYL5QktafNeRWZ+aWIOH2N1/+sb/PL9LJB7TCbJWlGVMrmgT2Pmfm+zNyRmacDlwBfyMx3Al8E3t4ctgu4YZSGSJJGk2x8Nbcm2E48MpSxeVw+QnMuA/7kOc2Dz0fEHSN+rjCbJWlWVMvmUe7z+KvAdRHxG8BXgY+P8FmSpMl7dBxDSiPidfQC6u/07T4/Mw9ExEnAzRHx9cz80qjfpecxmyWplqnK5nUVj5l5K3Br8/xB4Lz1vF+S1K5J34c4Is6mt/rnmzPzsSP7M/NA8+fBiPgsvfyweBwDs1mSplulbPaui5JURU72RsQR8RLgD4Gfz8y/6Nt/TES84Mhz4AJgxVXhJEkqpVg2jzJsVZI0bVq8vBkR1wKvpTf/Yh74ALANIDM/Cryf3u0ifq93y0EWmqE2JwOfbfZtBT6VmZ9rr6WSJE2RQtncbfGYwEJ7qw3Niph037U0YS0uOjY7ZvD3QGZeOuD1XwR+cYX9DwIvf/47JK2m+QedCssh78MgraXrbLbnUZIKaXM5cEmStH6VstniUZIK8UK2JEnTpVI2WzxKUhFJraubkiTNumrZbPEoSVUkTiiVJGmaFMvmTovHWIDtj07R3UEKdSFLGsIU/e6OhUm3QOqJCOKooybdDGnTmaJIUiOe9b/KIPY8SlIhleZVSJJUQaVstniUpEoKBZQkSSUUymaLR0kqI0pNypckafbVymaLR0mqpNDVTUmSSiiUzZ0Wj1sPwwn3LXb5le0p9EMgzZQiF+/mD0+6BVJj6xxbjnvRpFshSZN3eG7SLZh69jxKUhVZ615SkiTNvGLZbPEoSZU4KkKSpOlSKJstHiWplDpXNyVJqqFONm+ZdAMkSZIkSdOv057HuacXeOFdj3X5lZI0leaeXmjngwsNjVE3lo7exuGfPPk5+wpNz5GkVcWyzFx6fFs7X1Qomx22KkmVFAooSZJKKJTNFo+SVEVil5EkSdOkWDZbPEpSIVno6qYkSRVUymYXzJEkSZIkDdRtz+MPFuCRb3f6lZI0lX7ggjmaDgs/Gjx69tGTbsbMKDT6TJvE8kVhtLqFu1r6H7zQfwOHrUpSJf7LVpKk6VIomy0eJakQrzBLkjRdKmWzxaMkVZGUGhojSdLMK5bNLpgjSZIkSRqo257HpSXy8DOdfqUkTaWlpRY+NErNq1A3lrbBU6e28fMoSbNlaVsbn1ormx22KkmVFBoaI0lSCYWy2eJRkiopFFCSJJVQKJud8yhJkiRJGsieR0mqpNDVTUmSSiiUzZ0WjwlkFvrbk6QNauU3YVJqUr46sgWWfsRslqRWxmQWy2Z7HiWpkEo3IpYkqYJK2eycR0mqJEd4DBAReyLiYETcvcrrERH/NiL2R8TXIuIVfa/tiogHmseuUU5RkqSZUiibLR4lScO6GrhwjdffDJzZPC4HPgIQEScAHwBeCZwHfCAijm+1pZIkbQ5X02E2WzxKkoaSmV8CDq1xyEXAJ7Lny8BxEXEK8Cbg5sw8lJmPAzezdtBJkqQhdJ3NznmUpEImPK/iVOChvu35Zt9q+yVJKq9SNls8SlIlo63odmJE7Ovb3p2Zu9fx/pW+PNfYL0lSfYWy2eJRkqoYcnL9Gh7NzJ0jvH8eOK1vewdwoNn/2mX7bx3heyRJmg3Fstk5j5KkcdkL/ONmZbdXAd/JzIeBm4ALIuL4ZjL+Bc0+SZLUrrFmsz2PklRJi4NBI+JaelcpT4yIeXqrtG0DyMyPAjcCbwH2A98H3tW8digifh24vfmoKzNzrcn9kiTVUSibBxaPEbEd+BJwdHP89Zn5gYg4A7gOOAH4CvDzmfns8KcqSRq3NiflZ+alA15P4N2rvLYH2NNGuzYjs1mSZkelbB6m5/EZ4PWZ+WREbAP+NCL+BPhnwIcy87qI+ChwGc19QyRJE+IyNJvFeLPZnxtJak+h37ED5zw29wR5stnc1jwSeD1wfbP/GuDiVlooSZKew2yWJE3CUAvmRMRcRNwJHKR3A8lvAE9k5kJzyKr3BYmIyyNiX0Ts+0EeHkebJUmryREeminjyubFJ59c6RBJ0rgUyuahisfMXMzMc+gt4XoecNZKh63y3t2ZuTMzd26L7RtvqSRpTZGjPTRbxpXNc8ce22YzJWlTq5bN61ptNTOfiIhbgVcBx0XE1uYK55H7hUiSJmm0GxFrBpnNkjTlCmXzMKut/hjwgyacfgT4OeCDwBeBt9Nb1W0XcMPAz+p93kgNlqQKWvtNOIVXKTV+48xmElhqr62SNDPaytBC2TxMz+MpwDURMUdvmOtnMvOPI+Je4LqI+A3gq8DHW2ynJEn6IbNZktS5gcVjZn4NOHeF/Q/Sm2MhSZoS0zg/QuNnNkvS7KiUzeua8yhJmnKFAkqSpBIKZbPFoyRVMaUrs0mStGkVy2aLR0mSNrlYdDE7SdJgFo+SVEmhq5uSJJVQKJstHiWpkkIBJUlSCYWy2eJRkgqpNK9CkqQKKmXzlkk3QJIkSZI0/brteYyALdarkkS4QImmREIsTboRkjQFCvUQtsVhq5JUicEnSdJ0KZTNFo+SVEWxe0lJkjTzimWzxaMkVVIooCRJKqFQNls8SlIlhQJKkqQSCmWzxaMkSZvdkgs4SZIGs3iUpCKCWvMqJEmaddWy2eJRkiopFFCSJJVQKJstHiWpimIrukmSNPOKZfOWSTdAkiRJkjT97HmUpEoKXd1URxJiadKNkKQp0FaGFspmi0dJqqRQQEmSVEKhbLZ4lKRCKs2rkCSpgkrZbPEoSZUUCihJkkoolM0umCNJkiRJGqj7nseIzr9SkjaFpNTVTXWj2g2sJWmjWqlSimWzw1YlqRCLAEmSpkulbHbYqiRVkiM8hhARF0bE/RGxPyKuWOH1D0XEnc3jLyLiib7XFvte2zvKaUqSNDMKZbM9j5JUSJtXNyNiDvgw8EZgHrg9IvZm5r1HjsnMX+k7/peAc/s+4unMPKe9FkqSNH0qZbM9j5KkYZ0H7M/MBzPzWeA64KI1jr8UuLaTlkmStDl1ms32PEpSJaNd3TwxIvb1be/OzN1926cCD/VtzwOvXOmDIuIngDOAL/Tt3t58/gJwVWb+0Uit1fgsTboBklRYoWy2eJSkKkZf0e3RzNy5xusrLUS32jdeAlyfmYt9+16SmQci4qXAFyLirsz8xkYbK0nS1CuWzQ5blaQiYsTHEOaB0/q2dwAHVjn2EpYNi8nMA82fDwK38tw5F5IklVMtmy0eJUnDuh04MyLOiIij6IXQ81Zmi4ifAo4H/kvfvuMj4ujm+YnA+cC9y98rSZLWpdNsdtiqJFXS4opumbkQEe8BbgLmgD2ZeU9EXAnsy8wjYXUpcF1m9rfmLOD3I2KJ3oXLq/pXgpMkqaxC2dx98bjFzk5JakvbNyLOzBuBG5fte/+y7f9jhff9GfC3Wm2cNiYhXDBHklor8iplsz2PklRJywElSZLWqVA2WzxKUiWFAkqSpBIKZbNjSCVJkiRJA9nzKElVZPvzKiRJ0joUy2aLR0mqpFBASZJUQqFstniUpEIqXd1Uh/y5kaTWVMpmi0dJqqRQQEmSVEKhbHbBHEmSJEnSQPY8SlIhlYbGSJJUQaVsHtjzGBGnRcQXI+K+iLgnIt7b7D8hIm6OiAeaP49vv7mSpFXliA/NDLNZkmZEsWwepudxAfjnmfmViHgBcEdE3Az8AnBLZl4VEVcAVwC/2l5TJUkDTWHQqBVjzeZKV8UlaeoU+h07sOcxMx/OzK80z78H3AecClwEXNMcdg1wcVuNlCRJP2Q2S5ImYV1zHiPidOBc4Dbg5Mx8GHohFhEnrfKey4HLAbbHMaO0VZK0hsAepM1o1Gze+iJHtkpSW6pl89CrrUbEscAfAL+cmd8d9n2ZuTszd2bmzqNi+0baKEkaVqF5FRpsHNk8d4wXdiWpVYWyeaiex4jYRi+cPpmZf9jsfiQiTmmubJ4CHGyrkZKk4UROYdKoFWazJM2GStk8zGqrAXwcuC8zf6vvpb3Arub5LuCG8TdPkjS0Yiu6aXVmsyTNiGLZPEzP4/nAzwN3RcSdzb5fA64CPhMRlwHfAt7RThMlSdIyZrMkqXMDi8fM/FN6cz1X8obxNkeSNIpKk/K1OrNZkmZHpWxe12qrkqQpVyigJEkqoVA2WzxKUiGVrm5KklRBpWy2eJSkSgoFlLqTqw2AlSSNrlA2D32fR0mSJEnS5mXPoyRVkbWGxkiSNPOKZbPFoyRVUiigJEkqoVA2WzxKUhFBraubkiTNumrZbPEoSdJm54I5kqQhWDxKUiVZ6PKmJEkVFMpmi0dJKqTS0BhJkiqolM0Wj5JURVJqUr4kSTOvWDZbPEpSIbE06RZIkqR+lbK5++JxqdDfniRJsy5wwRxJAn8XDsGeR0mqpNDQGEmSSiiUzVsm3QBJ0vhEbvwx1OdHXBgR90fE/oi4YoXXfyEivh0RdzaPX+x7bVdEPNA8do3vrCVJml6VstmeR0mqIml1OfCImAM+DLwRmAduj4i9mXnvskM/nZnvWfbeE4APADublt7RvPfx1hosSdKkFctmex4lqZCWr26eB+zPzAcz81ngOuCiIZv2JuDmzDzUhNLNwIUbOUdJkmZJpWy251GSdMSJEbGvb3t3Zu7u2z4VeKhvex545Qqf8/cj4tXAXwC/kpkPrfLeU8fTbI0qvZQsSdNqqrLZ4lGSKhltZMyjmblzjddXWodu+Tf+P8C1mflMRPwT4Brg9UO+V5Kkegpls9caJamIoPWhMfPAaX3bO4AD/Qdk5mOZ+Uyz+e+Anx32vZIkVVMtmy0eJamKzNEeg90OnBkRZ0TEUcAlwN7+AyLilL7NtwH3Nc9vAi6IiOMj4njggmafJEl1Fctmh61KkoaSmQsR8R56wTIH7MnMeyLiSmBfZu4F/reIeBuwABwCfqF576GI+HV6IQdwZWYe6vwkJEkqpOtstniUpEKGvSfURmXmjcCNy/a9v+/5+4D3rfLePcCeVhuodUsgV5r1IkmbTFsRWimbLR4lqRKXoJEkaboUymaLR0kqpO2rm5IkaX0qZbPFoyRVkcBSoYSSJGnWFctmV1uVJEmSJA3Ufc/jcEvOSpI2wjgPLbAAAA4PSURBVF+xWq/AS8mSBL3fh20olM0OW5WkQirNq5AkqYJK2WzxKEmVOLpDkqTpUiibLR4lqZBKVzclSaqgUjY7y0GSJEmSNJA9j5JURVJqUr4kSTOvWDZbPEpSEQFEoXkV6k7O+XMjSW2ols0Wj5JUydKkGyBJkp6jUDY751GSJEmSNJA9j5JUSKWhMZIkVVApmy0eJamKYpPyJUmaecWyudviMROWCg36laSNauUqZJa6EbE6EpBOYpGk3uo2Y1crm+15lKRCKt2IWJKkCipls9caJUmSJEkD2fMoSZUUGhojSVIJhbJ5YM9jROyJiIMRcXffvhMi4uaIeKD58/h2mylJGighljb+0OwwmyVpRhTL5mF6Hq8Gfhf4RN++K4BbMvOqiLii2f7V8TdPkrQuha5uak1XM8ZsdsEcSWpRoWweGBeZ+SXg0LLdFwHXNM+vAS4ec7skSdIqzGZJ0iRsdM7jyZn5MEBmPhwRJ612YERcDlwOsD2O2eDXSZKGUufiptZvQ9k8d7yjWyWpVYWyufUFczJzN7Ab4EVbXlzor06Spk8UGhqj9vRn89EvOc0fGklqUaVs3ugsh0ci4hSA5s+D42uSJGnDMjf+0KwzmyVpGhXK5o32PO4FdgFXNX/eMMybEsgp/EuQpK618pswgSlcmU2d2VA2A7DFbJakVhTL5mFu1XEt8F+An4qI+Yi4jF4wvTEiHgDe2GxLkqQOmM2SpEkY2POYmZeu8tIbxtwWSdIIgiw1r0KrM5slaTZUy+bWF8yRJHWoUEBJklRCoWy2eJSkSgoFlCRJJRTKZotHSaqi2KR8dSTY+NrrklRJtPCZxbLZuJAkSZIkDWTxKEmFROaGH0N9fsSFEXF/ROyPiCtWeP2fRcS9EfG1iLglIn6i77XFiLizeewd42lLkjS1KmWzw1YlqZIW51VExBzwYXq3gZgHbo+IvZl5b99hXwV2Zub3I+KfAv8X8D83rz2dmee01kBJkqZRoWy251GSysheQG30Mdh5wP7MfDAznwWuAy56Tgsyv5iZ3282vwzsGOspSpI0U2pls8WjJOmIEyNiX9/j8mWvnwo81Lc93+xbzWXAn/Rtb28+98sRcfGY2qyRJRnPfUjSZvD8331T+ftvqrLZYauSVEUy6tCYRzNz5xqvr7QO3YpfGBH/CNgJvKZv90sy80BEvBT4QkTclZnf2HhzJUmacsWy2eJRkippdznweeC0vu0dwIHlB0XEzwH/O/CazHzmyP7MPND8+WBE3AqcC1g8SpJqK5TNDluVpEJaXtHtduDMiDgjIo4CLgGeszJbRJwL/D7wtsw82Lf/+Ig4unl+InA+0D+ZX5Kkkiplsz2PklRJiyu6ZeZCRLwHuAmYA/Zk5j0RcSWwLzP3Av8GOBb4jxEB8K3MfBtwFvD7EbFE78LlVctWgpMkqaZC2dxp8RhA02BJ2tRm9TdhZt4I3Lhs3/v7nv/cKu/7M+Bvtds6bUgGW551IJKkzSeWp3HOZjp3mc32PEpSFQksTeVKcZIkbU7FstniUZLKGPqeUJIkqRO1stniUZIqKRRQkiSVUCibLR4lqZJCASVJUgmFsrnb4jECtm3r9CslaSo9M5uT8lVPLMDRjxZZMKfOv8+k2VIk0mJh0i2YfvY8SlIVxSblS5I084pls8WjJJWRkEuTboQkSfprtbLZ4lGSKik0r0KSpBIKZXORSQ6SJEmSpDZ12/O4dY4tJxzX6VdK0lQ6PDf+zyw2r0Ld2HoYTrhvcdLN+Gvhj7C0qeQULbYzf7iFDy2WzQ5blaRKCg2NkSSphELZbPEoSZUUCihJkkoolM0Wj5JURpYKKEmSZl+tbHbBHEmSJEnSQPY8SlIVCSzVuZeUJEkzr1g2d1o8Lh29jafPPKndL9kyRUs2SZpdLa+MtvTYtnY+uNDQGHVj7ukFXnj3Y5Nuxg/5MyxtLjE9/3afe3qhnQ8u9HvNnkdJqqRQQEmSVEKhbLZ4lKQystS9pCRJmn21stkFcyRJkiRJA9nzKElVJGTWmZQvSdLMK5bNnRaPCz8aPHr20V1+pSRNpYWvtbRAQKGhMerIDxbgkW9PuhWSNHk/aGnBnELZbM+jJFVSaFK+JEklFMpm5zxKkiRJkgay51GSqsgsdSNiSZJmXrFstniUpEoKDY2RJKmEQtncafG4tA2+/z/UqbwlaaOWtrXzuVno6qY6srREHn5m0q2QpMlrKUMrZbM9j5JURpa6uilJ0uyrlc0umCNJkiRJGsieR0mqIil1LylJkmZesWweqecxIi6MiPsjYn9EXDGuRkmSNiiXNv5QCWazJE2ZQtm84Z7HiJgDPgy8EZgHbo+IvZl576pv2gJLP1Kn8pakDWth0kACWejqptZvI9mcQBaajyNJG9XGb8Jq2TzKP1/OA/Zn5oOZ+SxwHXDReJolSVq3zNavbg7q1YqIoyPi083rt0XE6X2vva/Zf39EvGls561+ZrMkTZNi2TxK8Xgq8FDf9nyzb3ljL4+IfRGxb/HJJ0f4OknSJPX1ar0ZeBlwaUS8bNlhlwGPZ+b/CHwI+GDz3pcBlwA/DVwI/F7zeRqvdWfzD/JwZ42TJI1X19k8SvEYK+x7Xp9sZu7OzJ2ZuXPu2GNH+DpJ0iC5lBt+DGGYXq2LgGua59cDb4iIaPZfl5nPZOZfAvubz9N4rTubt8X2DpolSZtXpWwepXicB07r294BHBjh8yRJo2p3aMwwvVp/fUxmLgDfAV485Hs1OrNZkqZNoWwe5VYdtwNnRsQZwF/R6/L8h2u94dlvzT/63/7Xf/nfgBOBR0f47mniuUynKudS5TzAc1nuJ8bRkH7f4/Gb/t+8/sQRPmJ7ROzr296dmbv7tofp1VrtmKF6xDSydWfz9/LQozcf/qTZPL08l+lT5TzAc1nObB6QzRsuHjNzISLeA9wEzAF7MvOeAe/5MYCI2JeZOzf63dPEc5lOVc6lynmA59KFzLyw5a8YplfryDHzEbEVeBFwaMj3akRmc4/nMp2qnEuV8wDPpQvVsnmkxeIz88bM/MnM/BuZ+ZujfJYkaer9da9WRBxFr1dr77Jj9gK7mudvB76QvftA7AUuaVZ8OwM4E/jzjtq9qZjNkrSpdJrNowxblSRtIqv1akXElcC+zNwLfBz4DxGxn95VzUua994TEZ8B7gUWgHdn5uJETkSSpCK6zuZJFY+7Bx8yMzyX6VTlXKqcB3guJWTmjcCNy/a9v+/5YeAdq7z3NwF7wqZXpZ9rz2U6VTmXKucBnksJXWZz9HosJUmSJEla3UhzHiVJkiRJm0PnxWNEXBgR90fE/oi4ouvvH0VE7ImIgxFxd9++EyLi5oh4oPnz+Em2cRgRcVpEfDEi7ouIeyLivc3+WTyX7RHx5xHxX5tz+T+b/WdExG3NuXy6mUA89SJiLiK+GhF/3GzP6nl8MyLuiog7jywvPYs/XwARcVxEXB8RX2/+n/nbs3ou0mrM5skzm6eX2Tx9zObJ6bR4jIg54MPAm4GXAZdGxMu6bMOIrgaWL7d7BXBLZp4J3NJsT7sF4J9n5lnAq4B3N/8dZvFcngFen5kvB84BLoyIVwEfBD7UnMvjwGUTbON6vBe4r297Vs8D4HWZeU7fstmz+PMF8DvA5zLzbwIvp/ffZ1bPRXoes3lqmM3Ty2yePmbzhHTd83gesD8zH8zMZ4HrgIs6bsOGZeaX6K1Q1O8i4Jrm+TXAxZ02agMy8+HM/Erz/Hv0/oc7ldk8l8zMJ5vNbc0jgdcD1zf7Z+JcImIH8HeBjzXbwQyexxpm7ucrIl4IvJreKmVk5rOZ+QQzeC7SGszmKWA2TyezefqYzZPVdfF4KvBQ3/Z8s2+WnZyZD0PvFz9w0oTbsy4RcTpwLnAbM3ouzXCSO4GDwM3AN4AnMnOhOWRWfs5+G/hXwFKz/WJm8zyg94+Ez0fEHRFxebNvFn++Xgp8G/j3zZClj0XEMczmuUirMZunjNk8Vczm6WM2T1DXxWOssM/lXickIo4F/gD45cz87qTbs1GZuZiZ5wA76F1BP2ulw7pt1fpExFuBg5l5R//uFQ6d6vPoc35mvoLeMLh3R8SrJ92gDdoKvAL4SGaeCzyFw2BUzyz/rinHbJ4eZvPUMpsnqOvicR44rW97B3Cg4zaM2yMRcQpA8+fBCbdnKBGxjV44fTIz/7DZPZPnckQzZOFWenNFjouII/cxnYWfs/OBt0XEN+kNGXs9vauds3YeAGTmgebPg8Bn6f3DYRZ/vuaB+cy8rdm+nl5gzeK5SKsxm6eE2Tx1zObpZDZPUNfF4+3Amc0qVUcBlwB7O27DuO0FdjXPdwE3TLAtQ2nG638cuC8zf6vvpVk8lx+LiOOa5z8C/By9eSJfBN7eHDb155KZ78vMHZl5Or3/L76Qme9kxs4DICKOiYgXHHkOXADczQz+fGXmfwceioifana9AbiXGTwXaQ1m8xQwm6eP2TydzObJisxue9oj4i30rtrMAXsy8zc7bcAIIuJa4LXAicAjwAeAPwI+A7wE+BbwjsxcPnF/qkTE3wH+E3AXPxzD/2v05lbM2rmcTW9S9By9iyGfycwrI+Kl9K4SngB8FfhHmfnM5Fo6vIh4LfAvMvOts3geTZs/22xuBT6Vmb8ZES9mxn6+ACLiHHoLJRwFPAi8i+ZnjRk7F2k1ZvPkmc3TzWyeLmbz5HRePEqSJEmSZk/Xw1YlSZIkSTPI4lGSJEmSNJDFoyRJkiRpIItHSZIkSdJAFo+SJEmSpIEsHiVJkiRJA1k8SpIkSZIGsniUJEmSJA30/wOBfVy3rDeQEgAAAABJRU5ErkJggg==\n",
Michael Kuron's avatar
Michael Kuron committed
379
      "text/plain": [
380
       "<Figure size 1152x432 with 4 Axes>"
Michael Kuron's avatar
Michael Kuron committed
381
382
      ]
     },
383
384
385
     "metadata": {
      "needs_background": "light"
     },
Michael Kuron's avatar
Michael Kuron committed
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
     "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",
412
   "version": "3.7.2"
Michael Kuron's avatar
Michael Kuron committed
413
414
415
416
417
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}