Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
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
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
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
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
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
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
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from pystencils.session import *\n",
"sp.init_printing()\n",
"frac = sp.Rational"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Phase-field simulation of dentritic solidification in 3D\n",
"\n",
"This notebook tests the model presented in the dentritic growth tutorial in 3D. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"target = 'cpu'\n",
"gpu = target == 'gpu'\n",
"domain_size = (25, 25, 25) if 'is_test_run' in globals() else (300, 300, 300)\n",
"\n",
"dh = ps.create_data_handling(domain_size=domain_size, periodicity=True)\n",
"φ_field = dh.add_array('phi', latex_name='φ')\n",
"φ_delta_field = dh.add_array('phidelta', latex_name='φ_D')\n",
"t_field = dh.add_array('T')"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"ε, m, δ, j, θzero, α, γ, Teq, κ, τ = sp.symbols(\"ε m δ j θ_0 α γ T_eq κ τ\")\n",
"εb = sp.Symbol(\"\\\\bar{\\\\epsilon}\")\n",
"discretize = ps.fd.Discretization2ndOrder(dx=0.03, dt=1e-5)\n",
"\n",
"φ = φ_field.center\n",
"T = t_field.center\n",
"d = ps.fd.Diff\n",
"\n",
"def f(φ, m):\n",
" return φ**4 / 4 - (frac(1, 2) - m/3) * φ**3 + (frac(1,4)-m/2)*φ**2\n",
"\n",
"\n",
"\n",
"bulk_free_energy_density = f(φ, m)\n",
"interface_free_energy_density = ε ** 2 / 2 * (d(φ, 0) ** 2 + d(φ, 1) ** 2 + d(φ, 2) ** 2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here comes the major change, that has to be made for the 3D model: $\\epsilon$ depends on the interface normal, which can not be computed simply as atan() as in the 2D case"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$$\\bar{\\epsilon} \\left(δ \\left(\\frac{{\\partial_{0} {{φ}_{C}}}^{4}}{\\left({\\partial_{0} {{φ}_{C}}}^{2} + {\\partial_{1} {{φ}_{C}}}^{2} + {\\partial_{2} {{φ}_{C}}}^{2}\\right)^{2}} + \\frac{{\\partial_{1} {{φ}_{C}}}^{4}}{\\left({\\partial_{0} {{φ}_{C}}}^{2} + {\\partial_{1} {{φ}_{C}}}^{2} + {\\partial_{2} {{φ}_{C}}}^{2}\\right)^{2}} + \\frac{{\\partial_{2} {{φ}_{C}}}^{4}}{\\left({\\partial_{0} {{φ}_{C}}}^{2} + {\\partial_{1} {{φ}_{C}}}^{2} + {\\partial_{2} {{φ}_{C}}}^{2}\\right)^{2}}\\right) + 1\\right)$$"
],
"text/plain": [
" ⎛ ⎛ 4 \n",
" ⎜ ⎜ D(phi_C) D(phi_C\n",
"\\bar{\\epsilon}⋅⎜δ⋅⎜──────────────────────────────────── + ────────────────────\n",
" ⎜ ⎜ 2 \n",
" ⎜ ⎜⎛ 2 2 2⎞ ⎛ 2 \n",
" ⎝ ⎝⎝D(phi_C) + D(phi_C) + D(phi_C) ⎠ ⎝D(phi_C) + D(phi_C\n",
"\n",
" 4 4 ⎞ ⎞\n",
") D(phi_C) ⎟ ⎟\n",
"──────────────── + ────────────────────────────────────⎟ + 1⎟\n",
" 2 2⎟ ⎟\n",
" 2 2⎞ ⎛ 2 2 2⎞ ⎟ ⎟\n",
") + D(phi_C) ⎠ ⎝D(phi_C) + D(phi_C) + D(phi_C) ⎠ ⎠ ⎠"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = sp.Matrix([d(φ, i) for i in range(3)])\n",
"nLen = sp.sqrt(sum(n_i**2 for n_i in n))\n",
"n = n / nLen\n",
"nVal = sum(n_i**4 for n_i in n)\n",
"σ = δ * nVal\n",
"\n",
"εVal = εb * (1 + σ)\n",
"εVal"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def m_func(temperature):\n",
" return (α / sp.pi) * sp.atan(γ * (Teq - temperature))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"substitutions = {m: m_func(T),\n",
" ε: εVal}\n",
"\n",
"fe_i = interface_free_energy_density.subs(substitutions)\n",
"fe_b = bulk_free_energy_density.subs(substitutions)\n",
"\n",
"μ_if = ps.fd.expand_diff_full(ps.fd.functional_derivative(fe_i, φ), functions=[φ])\n",
"μ_b = ps.fd.expand_diff_full(ps.fd.functional_derivative(fe_b, φ), functions=[φ])"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"dF_dφ = μ_b + sp.Piecewise((μ_if, nLen**2 > 1e-10), (0, True))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7UAAAAXBAMAAAAmS3V+AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZpkQ3Ynvq81UMrtEdiLw+n06AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJeklEQVRoBe1ZW2ycRxX+1r/3fvHmIWlf0mxcTIpayCoE05fKVqFKQYhuojpFrVK7Kre2KCwGYkAYbxWhCATq8kAj8UCWS1VVPNiCquUmskIogievUhVQJcerSK14apK6Td2k7XLmzMzO5Z+1eXBUqWQe/nP7zpkz5/z//PPvAkjsxvXxPqtA7sNVsaLorvfZuq4vhypQeFSUIXExVozS1ANNS5mrSCHzQu5DwL6p7wHnb/0j6aZebIGV7CDNd7x4m8L8enZqqhaN3deV3uZauvz8wtW5J2pG43LR2PaWq/EkDdDUM7P4ZNVok0fuN4LH6RihRF3ovoODo0AbX576g+sVkzRSljBmdhTRyqG2o7AFkzmViwW+qObdIqDJGduB+fPIXmIm+7Agw0uA4FK9Xg3RQ3ilG9Ux3cbt7aGGVLIDm1FBuskY7Or1evXxcukoh7IuyRpSkxgqWyrJyumQaUafjdlshQZIqrxsAPE3GDn6EXilRmNxOlYoUQMTd8oPsL9rNB6njKUKTsXX5WAVUpbQsViCWtMQ8G1L67I6c6bmopp3rwAHervaxOvC9PE9a4L8c0lymfvOAekZpCvpKlIz+C3yS2AlOzCXLaPUYAw+QmAcA2ZFDHukgQsVJNq2TvBqOuwAPufbTs7NmTgawFR7+R5GzlfwJyN5nI51LJCohkazC1UkZpBsaI1PtbH4FkYmfaMja6QsoWMygl7TT4D/GK3H6cyZmotq3mGBDvT2TDta50BF0dvSjfTcCo56Aow0kLxIj136SuFdIbOSHZhLLiGaZAxqwP2gONMCZo9fALuaKNKj4A2eDngeOBU3WlgNkFR5WXaPXe56ClvUsUKJGtxqFfk6snyrG63htLHwFUw0jDrAaSSXMGCXKrWmB4D9AzE6c6bmopo3qLe0Bcs9madIZuzeTiyheCWzRr3Nz4hpuaPswFzhajVZZgzZCp3oKrDaEjhrtAB6DgqWRrFqRW8Cy8241Wg0QNJNe/uc8YxzKlYwUYOm3o7UUXjLaFzOMm6yJ2skl9ANYklqTcu/w0FL67JOFVhQGm5ev7dFejX2KuQ6XJf++2tMeYqb+r3ddraFFD23bwvk2sgPj9CxKS2UdHvVFLdw6WAfkwe+TM9tm+ze+Lwrq2nliqLXqLcdFwB8bK5MKgZqgKKqDo5DtLIkZBl3/amzbccqBSfWoERVZtTbiQoK7wwKY4wbfHYcEN4GSSUMhGMQb5NkTPe2i1XHhgA5VXhB1IwvVDhuHvc23aAD7zPndr4qIoyTQOPlO6tMRdmiju5tvlp6mzZkZMRDvVybeBypJljJDpJL9o73MXTP0Yb8RDy9EqViDzWt7FLpDWCxZpuJp81ODAZqgKKh3g7hoT48Wi/jEfZ2L06sAYnqglBvL9SQo8xig8P0jaV7b44hlCKi+5xOGibMsr9IMkuQ7i2m12UjVAhFGORU4UFRM77UVPN2CeyTLeSaye5LjntSllKULYnMknzfEuJZ4Mf4qdiaDmBiDUN8jiYlpANxJ/92ua0w0UUg08lON53QQkiENzfV29cCvb1gBRF3xmKNDgKShnp7E8wMUa+KD1ZjOUiFjjUgUe0lelsJ95YhlvFwWztpWjxCyephIeUjqg0OVWvKfvce2QjHKAWdOdMHRSn4QlNxLzIV4Df8kZvH1133L3IxxBQ/s3u7t4vEylNXgKGKOBAmxO4MUgLssLebqGP1YYUpNki/86XpFhF3pCluYMgV6c3GASyMmRePBmywJ1flQU+GeJ3Obm0nmhF0rAGJauAGezJDrJ02xfe7dhR0O04a0SCphAOH6u0NKF7mRgSAOnOmzp4se/GNDnCmRY5PY9Jyp+35dFPINEVU6/f2W+rdKTpzll5mDSQugZXswFyqi5x4ZARGfBfTWIhnR4fE0FAroiPBKZ7ewjz2USuKBkiqvCwwsclJI99NvS0b0eV0LNIGEtVY6i0dgrJmM9AGRbWx1EU+BvoMRkzyGskl9KIYUa3pKD2NfiX6IJ05U3NRzcvUCZhu0GU37CepVzW9zc7NLXypwy8Ael/sFUmSW7aCV9Mz4rllJTswN0H2r0kMJkR44ARfncuFiiNqQa2IPkZ3iXnsQXuOGRogabC3YkvS49gGz6348tWTBRLVIVbpG2gJCbtI2sRUG0fW4r2NrohziR4aySXUyhiVaxIvtXQ5ZlQKnTlTc+Fe0KIEStzh0Rp4d1VejwH3dAUvp0gtSa4C3I3EUSx28UvgXI7etzNgJTswR88t/iIxWCTNjm7uXcQGfd7SyXelg30r/7CMqks7EP8dyam7Bkga7O2EmEGNRXrfaj5G+7GCiWo49ZZ+dKCf4QYMbaRSDfvrpQIMmwZpJJdwQDRdeNBzm+8OAvUzF+VigS+qeYeFm/jtojiJL3CI8YYg25B9h7/BZdlGVG/Pi8YWj0ePo3TX1J5JfBPjTbCSHZgrHEeiwxjaWCvAv6u3l+Of86erdD8dj+rRp2Xn5bR8K000MNyMHvV9zlRBPurYqgFMtZdI3IxTAq3g6Vr0iB+vb9Sxwola52R8CuMt9ZuOmUfPoYyJMk51PBA9A8Md8lBHYIXkEnpIAxKFp0r8vYqdgcxlJJ05U3PhXgD93ubL1CgxhuviWhi9tS3YoWPrtL8mT79ZZq40NksVOzTaRoY+hyeRHP0OHVaFkh2k+emx2yQGeKVJzqMkDn9CRO2P5O97JzpIz2+jE1luRqjltDwJsdGRf7V8n8L8TgOUAAWUXu4MwAcEWsXFoT2xeH2jjhVKtA/a89UTZdpnaLmfpAq4Q+bOxgPAX0f/7IPo9k11hA+ZaSikLGE8HIN0JXKzh9p+JUQMBunMmZoL98LqrcBfy/HzQPBUm/59aBbrAROrQj6DsICPvjMG9RExgFBsCsp2g36u0gXR+3ak7QL6kovsq11m06RcuJT6z23IuJW6WiDYSAvViW6e7+iAGSGfEE7qXPSvSvSZ5g0X4Rm1uCkoqZEbUQ/knJNdPw/pGrW0aVIaaFPubeD/WxuzFXzUCUTJN6N6qrXYCpiEKugzAOujS2/QXz/e+J/ibQ66wwsbFD3QWfv71nXwkK5RSZsnFXK7hZXfD5m2VJcIRYvmV6qFefpNMjyCPmEoaT30/O4Y0kPE7KzYHFQOO7paD5Qcu9m1G8lDGoPFbZ6UBdas+o02Gy+ERlx7euO1n+L/cYbcWPW9X7Y8nr/3eVzPYMsrMLTe2vKY1wOqCvwXoyAq8dhqiWoAAAAASUVORK5CYII=\n",
"text/latex": [
"$$\\left \\{ \\pi : 3.14159265358979, \\quad T_{eq} : 1.0, \\quad \\bar{\\epsilon} : 0.01, \\quad j : 6, \\quad α : 0.9, \\quad γ : 10, \\quad δ : 0.3, \\quad θ_{0} : 0.2, \\quad κ : 1.8, \\quad τ : 0.0003\\right \\}$$"
],
"text/plain": [
"{π: 3.14159265358979, T_eq: 1.0, \\bar{\\epsilon}: 0.01, j: 6, α: 0.9, γ: 10, δ:\n",
" 0.3, θ₀: 0.2, κ: 1.8, τ: 0.0003}"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"parameters = {\n",
" τ: 0.0003,\n",
" κ: 1.8,\n",
" εb: 0.01,\n",
" δ: 0.3,\n",
" γ: 10,\n",
" j: 6,\n",
" α: 0.9,\n",
" Teq: 1.0,\n",
" θzero: 0.2,\n",
" sp.pi: sp.pi.evalf()\n",
"}\n",
"parameters"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"dφ_dt = - dF_dφ / τ\n",
"assignments = [\n",
" ps.Assignment(φ_delta_field.center, discretize(dφ_dt.subs(parameters))),\n",
"]\n",
"φEqs = ps.simp.sympy_cse_on_assignment_list(assignments)\n",
"φEqs.append(ps.Assignment(φ, discretize(ps.fd.transient(φ) - φ_delta_field.center)))\n",
"\n",
"\n",
"temperatureEvolution = -ps.fd.transient(T) + ps.fd.diffusion(T, 1) + κ * φ_delta_field.center\n",
"temperatureEqs = [\n",
" ps.Assignment(T, discretize(temperatureEvolution.subs(parameters)))\n",
"]"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$$\\left [ {{T}_{C}} \\leftarrow 0.0111111111111111 {{T}_{B}} + 0.933333333333333 {{T}_{C}} + 0.0111111111111111 {{T}_{E}} + 0.0111111111111111 {{T}_{N}} + 0.0111111111111111 {{T}_{S}} + 0.0111111111111111 {{T}_{T}} + 0.0111111111111111 {{T}_{W}} + 1.8 \\cdot 10^{-5} {{φ_D}_{C}}\\right ]$$"
],
"text/plain": [
"[T_C := 0.0111111111111111⋅T_B + 0.933333333333333⋅T_C + 0.0111111111111111⋅T_\n",
"E + 0.0111111111111111⋅T_N + 0.0111111111111111⋅T_S + 0.0111111111111111⋅T_T +\n",
" 0.0111111111111111⋅T_W + 1.8e-5⋅phidelta_C]"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"temperatureEqs"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"φ_kernel = ps.create_kernel(φEqs, cpu_openmp=4, target=target).compile()\n",
"temperatureKernel = ps.create_kernel(temperatureEqs, cpu_openmp=4, target=target).compile()"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def time_loop(steps):\n",
" φ_sync = dh.synchronization_function(['phi'], target=target)\n",
" temperature_sync = dh.synchronization_function(['T'], target=target)\n",
" dh.all_to_gpu()\n",
" for t in range(steps):\n",
" φ_sync()\n",
" dh.run_kernel(φ_kernel)\n",
" temperature_sync()\n",
" dh.run_kernel(temperatureKernel)\n",
" dh.all_to_cpu()\n",
"\n",
"\n",
"def init(nucleus_size=np.sqrt(5)):\n",
" for b in dh.iterate():\n",
" x, y, z = b.cell_index_arrays\n",
" x, y, z = x - b.shape[0] // 2, y - b.shape[1] // 2, z - b.shape[2] // 2\n",
" b['phi'].fill(0)\n",
" b['phi'][(x ** 2 + y ** 2 + z ** 2) < nucleus_size ** 2] = 1.0\n",
" b['T'].fill(0.0)\n",
"\n",
"\n",
"def plot(slice_obj=ps.make_slice[:, :, 0.5]):\n",
" plt.subplot(1, 3, 1)\n",
" plt.scalar_field(dh.gather_array('phi', slice_obj).squeeze())\n",
" plt.title(\"φ\")\n",
" plt.colorbar()\n",
" plt.subplot(1, 3, 2)\n",
" plt.title(\"T\")\n",
" plt.scalar_field(dh.gather_array('T', slice_obj).squeeze())\n",
" plt.colorbar()\n",
" plt.subplot(1, 3, 3)\n",
" plt.title(\"∂φ\")\n",
" plt.scalar_field(dh.gather_array('phidelta', slice_obj).squeeze())\n",
" plt.colorbar()"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Name| Inner (min/max)| WithGl (min/max)\n",
"----------------------------------------------------\n",
" T| ( 0, 0)| ( 0, 0)\n",
" phi| ( 0, 1)| ( 0, 1)\n",
"phidelta| (inf,inf)| (inf,inf)\n",
"\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7UAAAF1CAYAAAA3Ls2oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X20ZXdd5/n3h0oIymOFJBjz0IkSp404BqkJzDDSSEgILKWwJW0ijcVMMpFepNVRZ3WQETDiWmCrqNNod8WUKRgkYFBT2tEyBBjUpTEViJAHY8o0kkpq5YEKCNiAVfWdP86+957c3Ft17z3n3HP3+b1fa+11z348v3tX1Wed7/n99m+nqpAkSZIkqY+eNO0GSJIkSZK0Vha1kiRJkqTesqiVJEmSJPWWRa0kSZIkqbcsaiVJkiRJvWVRK0mSJEnqLYtaSZIkSVJvWdRKkpqR5MtDy+Ek/31o/XXTbp8kSVo9i1pJUjOq6mlzC/A54PuHtr1/2u2TpFEl+e4k9yX5XJKfnHZ7pPVgUavHSfK/JPmbJC9J8kiS25KcPe12SZIkaUX+DjgH2Aq8Jcn3JDk5yZ8keXWSu5Lcn+RfT7md0thY1Gpekm8Efhd4D/A64P3AdcD7k2SabZMkSdLRVdWXq+ofq+pTDD7LvQL4DeAR4JnA14HXAjuTnDK9lkrjY1GrYecCAa4Gngp8Afh/gO8ETp9iuyRJkrQCSS5NckuSncAh4CnA9wG/ABwHfKGqbgH+Bnj59FoqjY9FrYY9B3igqmpuQ1V9FXgM+KaptUqSJElHleQ7gXcArwH+C3AZgx7aY4D7Fx2+Hz/faUZY1GrYA8Cpw0ONk3wDsBnYN7VWSZIkaSVeDvxhVe0HbgW+wuBWsq/zxFF3p+DnO80Ii1oNu4VB+L2JwTDkTcDPAX9RVQ9Ms2GSJEk6qn9k8PkN4K3AJ6rqHxgUtj8LPBkgyVbgXwI3TqOR0rgdM+0GaOOoqn/uQm4H8N0M7sP4c+D1U22YJEmSVuL9wA8kuQO4F/jfu+0/zmDOlF8EvgE4EfjXVfXYVFopjVmGbp+U5iX5f4G9VfX2abdFkiRJo0tyGfBvq+ql026LNE4OP5YkSZIk9ZZFraSxSrIjycPd0Kel9ifJryfZm+TTSb57vdsoSaNKcmGSe7osu3KJ/S9J8skkB5O8dtG+bUnu7ZZtQ9tfkOQz3TV/3WfES5q2vmSdw48ljVWSlwBfBt5bVc9bYv+rgH8PvAp4IfBrVfXC9W2lJK1dkk3A3wHnM5g99lbgkqq6a+iYM4BnAD8N7Kqq67vtxwN7gC1AAbcBL6iqx5L8NYN7H/+KwQQ+v15Vf7xOv5YkPU6fss6eWkljVVWfAA4c4ZCtDAreqqq/Ap6V5OT1aZ0kjcW5DOaduK+qvs5gZtmtwwdU1Wer6tPA4UXnvgK4qaoOdJP03ARc2OXgM6rqL7vnxb+XwbNGJWlaepN1FrWS1tspPP4B8Pu6bZLUF6Pk2HLnLn5mqNkoadp6k3VTeaTPCSecUGecccY03lqaGbfddtujVXXias55xfc+tT5/4NDa3/PTX7sT+OrQpu1VtX2Vl1nqvomZvA/CrJNGN42sg6Pm3Sg5tty5vc1Gs04anVk3mqkUtWeccQZ79uyZxltLMyPJP6z2nEcPHOKW3aeu+T2PPfnvv1pVW9Z8gYF9wGlD66cCD454zQ3JrJNGN42sg6Pm3Sg5tg946aJzP95tP3XR9l5ko1knjc6sG43DjyWtt13Aj3SzIL8I+GJV7Z92oyRpFW4FzkpyZpInAxczyLaV2A1ckGRzks3ABcDuLge/lORF3UygPwLcMInGS9IK9SbrptJTK2laikO1+D7+8UryAQbfzJ2QZB/wNuBYgKr6zwxmuXsVsBf4J+B/m2iDJDVosllXVQeTXMHgQ9smYEdV3ZnkKmBPVe1K8j8Bvw9sBr4/yc9V1XdU1YEkP8/gwyLAVVU1N7nevwOuBb4B+ONukaRlmHVzLGqlhhRweMK3aFXVJUfZX8CbJtoISU1bp6y7kcGXdMPb3jr0+lYeP8Ru+LgdwI4ltu8BnvAoNElailm3wKJWaszhJ8y4Lkmzx6yT1AKzbsB7aiVJkiRJvWVPrdSQojhUvXhChCStmVknqQVm3YIV99QmeUqSv07yN0nuTPJz3fYzk9yS5N4kH+xmxpK0QR2m1ry0wKyTZsMoWddC3pl10mww6wZWM/z4a8DLquq7gHOAC7vHcbwLeHdVnQU8Blw6/mZKGocCDlFrXhph1kk9N2rWNZJ3Zp3Uc2bdghUXtTXw5W712G4p4GXA9d32ncBrxtpCSWPlt3lHZtZJs8HeiyMz66TZYNYNrGqiqCSbktwOPAzcBPw98IWqOtgdsg84ZZlzL0+yJ8meRx55ZJQ2S9JEmXWSWmDWSZoVqypqq+pQVZ3D4FlE5wLfvtRhy5y7vaq2VNWWE088cfUtlTSyAg5VrXlphVkn9duoWddK3pl1Ur+ZdQvWNPtxVX0hyceBFwHPSnJM963eqcCDY2yfpDHzaWYrZ9ZJ/WXWrZxZJ/WXWTewmtmPT0zyrO71NwAvB+4GPga8tjtsG3DDuBspaTzKyQSOyqyT+m/UrGsh78w6qf/MugWr6ak9GdiZZBODYvhDVfVHSe4CrkvyDuBTwDUTaKekcSg4NDv5NSlmndR3Zt1KmHVS35l181Zc1FbVp4HnL7H9Pgb3YUhS75l1klpg1kmaJWu6p1ZSPxXeeyFp9pl1klpg1i2wqJWaEg6RaTdCkibMrJPUArNujkWt1JACDnvvhaQZZ9ZJaoFZt2BVz6mVJEmSJGkjsadWaozDVCS1wKyT1AKzbsCiVmpIYfhJmn1mnaQWmHULLGqlxhwuw0/S7DPrJLXArBuwqJUa4jd6klpg1klqgVm3wImiJEmSJEm9ZU+t1JAiHPK7LEkzzqyT1AKzboFFrdQY772Q1AKzTlILzLoBi1qpId57IakFZp2kFph1CyxqpaaEQ+UwFUmzzqyT1AKzbo5/BUmSJElSb9lTKzWkgMN+lyVpxpl1klpg1i2wqJUa470Xklpg1klqgVk3YFErNaTKey8kzT6zTlILzLoF/hUkSZIkSb1lT63UmMMOU5HUALNOUgvMugGLWqkhg+eZOUBD0mwz6yS1wKxb4F9Basrg3ou1LpLUD6Nl3UryLsmFSe5JsjfJlUvsPy7JB7v9tyQ5o9v+uiS3Dy2Hk5zT7ft4d825fSeN+Q8jaaaYdXPsqZUa4tTvklow6axLsgl4D3A+sA+4Ncmuqrpr6LBLgceq6rlJLgbeBfxQVb0feH93ne8Ebqiq24fOe11V7ZlY4yXNDLNugZ9uJUmSVudcYG9V3VdVXweuA7YuOmYrsLN7fT1wXpLFN79dAnxgoi2VpLXrTdbZUys15lA5oYCk2TeGrDshyXAvwvaq2t69PgW4f2jfPuCFi86fP6aqDib5IvBs4NGhY36IJ35A/O0kh4APA++oqhrt15A0y8y6AYtaqSFFnFBA0swbU9Y9WlVbltm31KfIxR/IjnhMkhcC/1RVdwztf11VPZDk6Qw+6L0eeO8q2iypIWbdAj/dSo05XE9a8yJJfTFK1q0g7/YBpw2tnwo8uNwxSY4BngkcGNp/MYuG41XVA93PLwG/w2DonyQty6wbsKdWaohTv0tqwTpk3a3AWUnOBB5g8KHthxcdswvYBvwl8Frgo3PD65I8CbgIeMncwd2HwWdV1aNJjgW+D/jIJH8JSf1m1i2wqJUkSVqF7r6xK4DdwCZgR1XdmeQqYE9V7QKuAd6XZC+DXouLhy7xEmBfVd03tO04YHf3IW8Tgw95V6/DryNJS+pT1lnUSg0p4kRRkmbeemRdVd0I3Lho21uHXn+VQQ/FUud+HHjRom1fAV4w9oZKmllm3QKLWqkxPqdWUgvMOkktMOsGLGqlhlTBISd8kjTjzDpJLTDrFvhXkCRJkiT1lj21UlPC4SUfJyZJs8Ssk9QCs27Ointqk5yW5GNJ7k5yZ5If77a/PckDSW7vlldNrrmSRlEMhqmsdWmBWSf136hZ10LemXVS/5l1C1bTU3sQ+Kmq+mSSpwO3Jbmp2/fuqvql8TdP0rj5nNqjMuukGWDWHZVZJ80As25gxUVtVe0H9nevv5TkbuCUSTVM0vgV4bCP9Dkis07qP7Pu6Mw6qf/MugVrKu2TnAE8H7il23RFkk8n2ZFk8zLnXJ5kT5I9jzzyyJoaK0nryayT1AKzTlLfrbqoTfI04MPAT1TVPwK/CXwrcA6Db/x+eanzqmp7VW2pqi0nnnjiCE2WNIpDPGnNS0vMOqnfRsm6lvLOrJP6zawbWNXsx0mOZRB876+q3wOoqoeG9l8N/NFYWyhpbAo4PEOTAkyKWSf1m1m3Mmad1G9m3YIVF7VJAlwD3F1VvzK0/eTuvgyAHwDuGG8TJY1POOTU70dk1kmzwKw7GrNOmgVm3ZzV9NS+GHg98Jkkt3fbfga4JMk5DL4s+Czwo0e70N/ddh/nP+miVTZV0rCns/kFqz3Hb/RWZGxZ95mHH+LMX1ty5J6kFXryaaeadZMxtqyTNB1m3YLVzH7857DkVwE3jq85kjRdZp2kFph1kmbJqu6pldR/DlOR1AKzTlILzLoBi1qpIVVxmIqkmWfWSWqBWbfAolZqzKEJh1+SC4FfAzYBv1VV71y0/3RgJ/Cs7pgrq8rhbpLGatJZJ0kbgVk34F9B0tgk2QS8B3glcDaDCUfOXnTY/w18qKqeD1wM/Mb6tlKSJEmzxJ5aqSEFHJ7svRfnAnur6j6AJNcBW4G7FjXjGd3rZwIPTrJBktqzDlknSVNn1i2wqJWakkkPUzkFuH9ofR/wwkXHvB340yT/Hngq8PJJNkhSiyaedZK0AZh1cyxqpYYMnmc20jd6JyTZM7S+vaq2D60vdfFatH4JcG1V/XKS/xl4X5LnVdXhURomSXPGkHWStOGZdQssaqXGHBrtVvpHq2rLEfbvA04bWj+VJw4vvhS4EKCq/jLJU4ATgIdHaZgkDRsx6ySpF8y6Af8KksbpVuCsJGcmeTKDiaB2LTrmc8B5AEm+HXgK8Mi6tlKSJEkzw55aqSFFJjpMpaoOJrkC2M3gcT07qurOJFcBe6pqF/BTwNVJ/k8GI2feUFWLhyhL0ppNOuskaSMw6xZY1EqNOTzhARrdM2dvXLTtrUOv7wJePNFGSGrepLNOkjYCs27AolZqSBUc8hs9STPOrJPUArNugUWt1BiHqUhqgVknqQVm3YD91ZIkSZKk3rKnVmrIYEIBv8uSNNvMOkktMOsWWNRKjTmEw1QkzT6zTlILzLoBi1qpIYX3XkiafWadpBaYdQvsr5YkSZIk9ZY9tVJTvPdCUgvMOkktMOvm+FeQGnOYrHmRpL4YJetWkndJLkxyT5K9Sa5cYv9xST7Y7b8lyRnd9jOS/Pckt3fLfx465wVJPtOd8+tJDF5JR2TWDdhTKzXEh3RLasGksy7JJuA9wPnAPuDWJLuq6q6hwy4FHquq5ya5GHgX8EPdvr+vqnOWuPRvApcDfwXcCFwI/PGEfg1JPWfWLbCnVmrM4XrSmhdJ6otRsm4FeXcusLeq7quqrwPXAVsXHbMV2Nm9vh4470i9EUlOBp5RVX9ZVQW8F3jNWn53Se0w6wb8lCpJkvREJyTZM7RcPrTvFOD+ofV93TaWOqaqDgJfBJ7d7TszyaeS/H9Jvmfo+H1HuaYkjdtMZJ3Dj6WGDB7S7fBjSbNtTFn3aFVtWWbfUhevFR6zHzi9qj6f5AXAHyT5jhVeU5LmmXULLGqlxjjhk6QWTDjr9gGnDa2fCjy4zDH7khwDPBM40A23+xpAVd2W5O+Bb+uOP/Uo15SkxzHrBhx+LDVk7iHda10kqQ9GzboV5N2twFlJzkzyZOBiYNeiY3YB27rXrwU+WlWV5MRu8hWSfAtwFnBfVe0HvpTkRd39aD8C3DCWP4ikmWTWLbCnVpIkaRWq6mCSK4DdwCZgR1XdmeQqYE9V7QKuAd6XZC9wgMGHQYCXAFclOQgcAt5YVQe6ff8OuBb4BgYzgTrzsaSp6VPWWdRKjXEWY0ktmHTWVdWNDB5FMbztrUOvvwpctMR5HwY+vMw19wDPG29LJc0ys27AolZqicOIJbXArJPUArNunkWt1JDCiaIkzT6zTlILzLoFFrVSY/xGT1ILzDpJLTDrBry5TpIkSZLUW/bUSg2Zm/pdkmaZWSepBWbdghX31CY5LcnHktyd5M4kP95tPz7JTUnu7X5unlxzJY3K59QemVknzYYJP7ux98w6aTaYdQOrGX58EPipqvp24EXAm5KcDVwJ3FxVZwE3d+uSNqDC4FsBs07quVGzrpG8M+uknjPrFqy4qK2q/VX1ye71l4C7gVOArcDO7rCdwGvG3UhJ43OYrHlpgVknzYZRsq6FvDPrpNlg1g2saaKoJGcAzwduAZ5TVfthEJDAScucc3mSPUn2/DNfW1trJWkdjZp1h778lfVqqiSt2ahZ98gjj6xXUyVpSaueKCrJ04APAz9RVf+YrKzCr6rtwHaAZ+T4Wu37ShqDckKBlRpH1h13+mlmnTQNZt2KjSPrtmzZYtZJ02DWzVtVUZvkWAbB9/6q+r1u80NJTq6q/UlOBh4edyMljYez5K2MWSf1m1m3Mmad1G9m3YLVzH4c4Brg7qr6laFdu4Bt3ettwA3ja56kcXMygSMz66TZ4OQpR2bWSbPBrBtYTU/ti4HXA59Jcnu37WeAdwIfSnIp8DngovE2UZLWlVknqQVmnaSZseKitqr+HJadIuu88TRH0iTNTf2u5Zl1Uv+ZdUdn1kn9Z9YtWPVEUZL6rQw/SQ0w6yS1wKwbsKiVGjNLzySTpOWYdZJaYNYNWNRKDSmnfpfUALNOUgvMugUrnv1YkiRJkqSNxp5aqTHeeyGpBWadpBaYdQMWtVJTnCVPUgvMOkktMOvmWNRKjfEbPUktMOsktcCsG7ColRpSOKGApNln1klqgVm3wImiJEmSJEm9ZU+t1JIaTP8uSTPNrJPUArNunkWt1Bgf0i2pBWadpBaYdQMWtVJDCicUkDT7zDpJLTDrFnhPrSRJkiSpt+yplZri88wktcCsk9QCs26ORa3UGCcUkNQCs05SC8y6AYtaqTHeeyGpBWadpBaYdQMWtVJDqgw/SbPPrJPUArNugRNFSZIkSZJ6y6JWaszhypoXSeqLUbJuJXmX5MIk9yTZm+TKJfYfl+SD3f5bkpzRbT8/yW1JPtP9fNnQOR/vrnl7t5w0xj+JpBlk1g04/FhqjBMKSGrBJLMuySbgPcD5wD7g1iS7ququocMuBR6rqucmuRh4F/BDwKPA91fVg0meB+wGThk673VVtWdyrZc0S8y6AXtqpcZUZc2LJPXFKFm3grw7F9hbVfdV1deB64Cti47ZCuzsXl8PnJckVfWpqnqw234n8JQkx43p15bUGLNuwKJWakgx0eCTpA1h1Kzr8u6EJHuGlsuH3uIU4P6h9X08vgficcdU1UHgi8CzFx3zg8CnquprQ9t+uxuO97NJDF5JyzLrFjj8WJIk6Ykeraoty+xb6gPY4kGARzwmyXcwGKZ3wdD+11XVA0meDnwYeD3w3pU3WZJWbSayzp5aqTE1wiJJfTFK1q0g7/YBpw2tnwo8uNwxSY4Bngkc6NZPBX4f+JGq+vv5Nlc90P38EvA7DIb+SdKyzLoBi1qpJeU9tZIaMGLWrSDvbgXOSnJmkicDFwO7Fh2zC9jWvX4t8NGqqiTPAv4r8Oaq+ou5g5Mck+SE7vWxwPcBd4z8t5A0u8y6eQ4/llpjl6ukFkww66rqYJIrGMzmuQnYUVV3JrkK2FNVu4BrgPcl2cug1+Li7vQrgOcCP5vkZ7ttFwBfAXZ3H/I2AR8Brp7cbyFpJph1gEWtpDFLciHwawyC6req6p1LHPNvgLcziOK/qaofXtdGStKIqupG4MZF29469PqrwEVLnPcO4B3LXPYF42yjJI2qL1lnUSs1ZpLDiFfyPLMkZwFvBl5cVY+N44HbkrSYt0xIaoFZN2BRKzVmkg/pZuh5ZgBJ5p5nNvyQ7v8DeE9VPTZoTz080RZJatKEs06SNgSzbsCiVmpIMfI3eick2TO0vr2qtg+tL/U8sxcuusa3AST5CwZDlN9eVX8ySqMkadgYsk6SNjyzboFFrdSSAkYLvyM9ywxW9jyzY4CzgJcymBr+z5I8r6q+MErDJGne6FknSRufWTfPR/pIGqeVPs/shqr656r6b8A9DIpcSZIkadXsqdWK7X7w9iPuf8U3n7NOLdEoJnzvxfzzzIAHGEzrvnhm4z8ALgGu7Z5T9m3AfRNtlaTmeJ+ZpBaYdQMr7qlNsiPJw0nuGNr29iQPJLm9W141mWZKGpsaYTnapasOMngu2W7gbuBDc88zS/Lq7rDdwOeT3AV8DPi/qurz4/r1RmXWSTNilKxr5EOieSfNALMOWF1P7bXAfwLeu2j7u6vql8bWIkkTlIlPKLCC55kV8JPdshFdi1kn9dzks25GXIt5J/WYWTdnxT21VfUJ4MAE2yJpPfht3hGZddKMsPfiqMw7aQaYdcB4Joq6IsmnuyEsm8dwPUnaiMw6Sa0w7yT1yqhF7W8C3wqcA+wHfnm5A5NcnmRPkj3/zNdGfFtJa1KD55mtdWnYmrLu0Je/sl7tkzRsxKwz746ed8NZ98gjj6xn+yTNMevmjVTUVtVDVXWoqg4DVwPnHuHY7VW1paq2HMtxo7ytpFE4RGXV1pp1m5721PVrpKTHc0jemqw074az7sQTT1zfRkpaYNYBIxa1SU4eWv0B4I7ljpW0UWSEpU1mndRHo2Sdedcx76QNz6yDVcx+nOQDwEuBE5LsA94GvDTJOQzq/M8CPzqBNmodHO0ZtOO6hs+y1UZn1klqhXknaVasuKitqkuW2HzNGNsiaT3M0FCTSTDrpBlh1h2VeSfNALMOWN1zaiXNAsNPUgvMOkktMOsAi1qpLQXM0Ex3krQks05SC8y6eRa1UmPKb/QkNcCsk9QCs25g1OfUSpIkSZI0NfbUSq3xGz1JLTDrJLXArAMsaqX2eO+FpBaYdZJaYNYBFrVSc+I3epIaYNZJaoFZN2BR24jdD94+7SYAR2/HK775nHVqSaMKh6lImn1mnaQWmHXznChKkiRJktRb9tRKTYn3XkhqgFknqQVm3RyLWqk1DlOR1AKzTlILzDrAolZqj+EnqQVmnaQWmHWA99RKkiRJknrMnlqpNX6jJ6kFZp2kFph1gEWt1JbCCQUkzT6zTlILzLp5FrWNONrzX9frObY+h3b6fEi3pBaYdZJaYNYNeE+t1JoaYZGkvhgl61aQd0kuTHJPkr1Jrlxi/3FJPtjtvyXJGUP73txtvyfJK1Z6TUl6ArMOsKiVJElalSSbgPcArwTOBi5Jcvaiwy4FHquq5wLvBt7VnXs2cDHwHcCFwG8k2bTCa0rSuulT1lnUSpIkrc65wN6quq+qvg5cB2xddMxWYGf3+nrgvCTptl9XVV+rqv8G7O2ut5JrStJ66k3WWdRKjUmtfZGkvhgl61aQd6cA9w+t7+u2LXlMVR0Evgg8+wjnruSakvQ4Zt2AE0VJrXGWPEktGD3rTkiyZ2h9e1Vt714vdfHFHw+XO2a57Ut1NPh1oqQjM+sAi1qpLU74JKkF48m6R6tqyzL79gGnDa2fCjy4zDH7khwDPBM4cJRzj3ZNSVpg1s1z+LEkSdLq3AqcleTMJE9mMBnKrkXH7AK2da9fC3y0qqrbfnE3Y+iZwFnAX6/wmpK0nnqTdfbUSq2xp1ZSCyaYdVV1MMkVwG5gE7Cjqu5MchWwp6p2AdcA70uyl0GvxcXduXcm+RBwF3AQeFNVHQJY6pqT+y0kzQSzDrCoVecV33zOUY/Z/eDtI19D0+eET5JaMOmsq6obgRsXbXvr0OuvAhctc+4vAL+wkmtK0pGYdQMWtVJrLGoltcCsk9QCsw7wnlpJkiRJUo/ZUyu1xm/0JLXArJPUArMOsKiVmrLCB21LUq+ZdZJaYNYtsKiVWjP6Q7olaeMz6yS1wKwDLGql9viNnqQWmHWSWmDWAU4UJUmSJEnqMXtqtWI+h3Y2eO+FpBaYdZJaYNYNrLinNsmOJA8nuWNo2/FJbkpyb/dz82SaKWlsaoSlEeadNANGybpG8s6sk2aAWQesbvjxtcCFi7ZdCdxcVWcBN3frkjaqWpgpby1LQ67FvJP6a8SsayjvrsWsk/rLrJu34qK2qj4BHFi0eSuws3u9E3jNmNolSVNj3klqgVknaVaMek/tc6pqP0BV7U9y0hjaJGmSZuhbuXVm3kl9YtatlVkn9YlZB6zjRFFJLgcuB3gK37hebytpMcNvooazbtNmb0WTpsasm6jhrDv99NOn3BqpYWYdMPojfR5KcjJA9/Ph5Q6squ1VtaWqthzLcSO+raS18r6LNVtR3g1n3aanPXVdGyhpgfeZrdmqs+7EE09c1wZKWmDWDYxa1O4CtnWvtwE3jHg9SdqozDtJLTDrJPXOah7p8wHgL4H/Icm+JJcC7wTOT3IvcH63Lkm9Zt5JaoFZJ2lWrPie2qq6ZJld542pLZLWwwwNNZkU806aAWbdUZl10gww64B1nChK0gYwY/dPSNKSzDpJLTDr5lnUSq0x/CS1wKyT1AKzDrColdpj+ElqgVknqQVmHTD67MeSJEmSJE2NPbVSQ4L3XkiafWadpBaYdQssaqXWGH6SWmDWSWqBWQdY1EptcZY8SS0w6yS1wKyb5z21ksYqyYVJ7kmyN8mVRzjutUkqyZb1bJ8kSZJmiz21Umsm+I1ekk3Ae4DzgX3ArUl2VdVdi457OvBjwC2Ta42kptl7IakFZh1gT63UnhphObpzgb1VdV9VfR24Dti6xHE/D/wi8NURfhNJWt4oWeeHREl9YdYBFrVSc1JrX4ATkuwZWi5fdPlTgPuH1vd12xbeP3k+cFpV/dEkf09JbRsl67xHTVJfmHUDDj+WWjNagD1aVUe6BzZHesckTwLeDbxhpFZI0tHM0Ic1SVqWWQfYUytpvPYBpw2tnwo8OLT+dOB5wMeTfBbSJW2dAAAWTElEQVR4EbDLyaIkSZK0VvbUSi2Z/P0TtwJnJTkTeAC4GPjh+bev+iJwwtx6ko8DP11VeybaKkltmbF7xSRpSWbdPItaqTGTvH+iqg4muQLYDWwCdlTVnUmuAvZU1a7JvbskLZile8UkaTlm3YDDj6XWTHiGvKq6saq+raq+tap+odv21qUK2qp6qb20kiZiijOCJjk+yU1J7u1+bl7muG3dMfcm2dZt+8Yk/zXJ3ya5M8k7h45/Q5JHktzeLZeN1lJJvWfWARa1UnOcIU9SC6Y8I+iVwM1VdRZwc7f++PYlxwNvA17I4HFobxv6QPhLVfUvgecDL07yyqFTP1hV53TLb43cUkm9ZtYNWNRKkiSN11ZgZ/d6J/CaJY55BXBTVR2oqseAm4ALq+qfqupjAN3zvj/JYNI9SdpoNkzWWdRKrZnSEBVJWlejD8k72nO5j+Q5VbUfoPt50hLHrOS53s8Cvp9BD8icH0zy6STXJxmebV5Si8w6wImipLZYnEpqwXiy7ojP5U7yEeCbltj1lhVe/2jP9T4G+ADw61V1X7f5D4EPVNXXkryRQc/Iy1b4fpJmjVk3z6JWakhYOlkkaZasR9ZV1cuXff/koSQnV9X+JCcDDy9x2D7gpUPrpwIfH1rfDtxbVb869J6fH9p/NfCuNTRd0oww6xY4/FiSJGm8dgHbutfbgBuWOGY3cEGSzd2kKRd020jyDuCZwE8Mn9B9aJzzauDuMbdbklZjw2SdPbVSaxx+LKkF0826dwIfSnIp8DngIoAkW4A3VtVlVXUgyc8Dt3bnXNVtO5XBsL6/BT6ZBOA/dbN//liSVwMHgQPAG9bzl5K0AZl1gEWt1BwfzSOpBdPMum7o3HlLbN8DXDa0vgPYseiYfSwzorCq3gy8eayNldRrZt2ARa3UGotaSS0w6yS1wKwDLGql9hh+klpg1klqgVkHOFGUJEmSJKnH7KmVWlLeUyupAWadpBaYdfMsaqXWGH6SWmDWSWqBWQdY1ErN8Rs9SS0w6yS1wKwbsKiVWmP4SWqBWSepBWYd4ERRkiRJkqQes6dWaozDVCS1wKyT1AKzbmAsRW2SzwJfAg4BB6tqyziuK2nMCoepjMCsk3rCrBuZeSf1gFk3b5w9td9bVY+O8XqSJsHwG5VZJ/WBWTcO5p200Zl1gPfUSpIkSZJ6bFxFbQF/muS2JJeP6ZqSxiwM7r1Y6yKzTuqDUbPOvAPMO2nDM+sWjGv48Yur6sEkJwE3JfnbqvrE8AFdIF4O8BS+cUxvK2nVZijApmBVWbdp8+ZptFESmHWjO2LeDWfd6aefPq02SjLrgDH11FbVg93Ph4HfB85d4pjtVbWlqrYcy3HjeFtJa5CqNS+tW23WbXraU9e7iZI6o2SdeXf0vBvOuhNPPHEaTZSEWTdn5KI2yVOTPH3uNXABcMeo15U0ATXi0jCzTuqRUbPOvDPvpD4w6+aNY/jxc4DfTzJ3vd+pqj8Zw3UlaSMx6yS1wryT1CsjF7VVdR/wXWNoi6R1MEuTAqwns07qF7Nu7cw7qT/MuoFxPqdWUh8YfpJaYNZJaoFZB1jUSs3xGz1JLTDrJLXArBuwqJVaY/hJaoFZJ6kFZh0wpkf6SJIkSZI0DfbUSi0ph6lIaoBZJ6kFZt08i1qpNYafpBaYdZJaYNYBFrVSU4Lf6EmafWadpBaYdQu8p1aSJEmS1Fv21EqtKb/Sk9QAs05SC8w6wKJWao7DVCS1wKyT1AKzbsCiVmpJ4YQCkmafWSepBWbdPItaqTE5PO0WSNLkmXWSWmDWDThRlCRJkiSpt+yplVrjMBVJLTDrJLXArAPsqZWak1r7Ikl9MUrWjZp3SY5PclOSe7ufm5c5blt3zL1Jtg1t/3iSe5Lc3i0ndduPS/LBJHuT3JLkjNFaKqnvzLoBi1qpJcVg6ve1LpLUB6Nm3eh5dyVwc1WdBdzcrT9OkuOBtwEvBM4F3rboA+Hrquqcbnm423Yp8FhVPRd4N/CuURsqqcfMunkWtVJj7KmV1IJp9l4AW4Gd3eudwGuWOOYVwE1VdaCqHgNuAi5cxXWvB85LkpFbK6m3zLoBi1pJkqQnOiHJnqHl8lWc+5yq2g/Q/TxpiWNOAe4fWt/XbZvz291wvJ8d+jA3f05VHQS+CDx7Fe2SpMVmIuucKEpqjT2uklowetY9WlVbltuZ5CPANy2x6y0rvP5SvQ5zrX5dVT2Q5OnAh4HXA+89yjmSWmTWARa1UlOCw4glzb71yLqqevmy7588lOTkqtqf5GTg4SUO2we8dGj9VODj3bUf6H5+KcnvMLgP7b3dOacB+5IcAzwTODD6byOpj8y6BQ4/lloy3ckEJGl9jJp1o+fdLmBb93obcMMSx+wGLkiyuZs05QJgd5JjkpwAkORY4PuAO5a47muBj1YZzlKzzLp59tRKkiSN1zuBDyW5FPgccBFAki3AG6vqsqo6kOTngVu7c67qtj2VwQe+Y4FNwEeAq7tjrgHel2Qvg16Li9fvV5KkJ9gwWWdRKzXG4ceSWjDNrKuqzwPnLbF9D3DZ0PoOYMeiY74CvGCZ636V7kOjJIFZN8eiVmqNRa2kFph1klpg1gEWtVJz7KmV1AKzTlILzLoBi1qpJQUcNv0kzTizTlILzLp5zn4sSZIkSeote2ql1viFnqQWmHWSWmDWARa1UnO890JSC8w6SS0w6wYsaqXWjP6gbUna+Mw6SS0w6wDvqZWak1r7sqLrJxcmuSfJ3iRXLrH/J5PcleTTSW5O8i/G/TtK0ihZZ8+HpL4w6wYsaiWNTZJNwHuAVwJnA5ckOXvRYZ8CtlTV/whcD/zi+rZSkiRJs8SiVmpJjbgc3bnA3qq6r6q+DlwHbH1cE6o+VlX/1K3+FXDqSL+TJC02atbNUO+FpBlm1s0bS1F7tOGGkjaGAKla87ICpwD3D63v67Yt51Lgj9f+G60vs07qh1GzboV5N9PMO2njM+sWjDxR1NBww/MZfIC9Ncmuqrpr1GtLmoDDI519QpI9Q+vbq2r70HqWOGfJxEzyb4EtwL8aqUXrxKyTema0rGuaeSf1iFkHjGf24/nhhgBJ5oYbGnzS7Hm0qrYcYf8+4LSh9VOBBxcflOTlwFuAf1VVXxtvEyfGrJPUCvNOUq+Mo6hdarjhC8dwXUkTMOGhJrcCZyU5E3gAuBj44ce9f/J84L8AF1bVw5NszJiZdVKPzNKwuikw76SeMOsGxlHUrmi4YZLLgcsBnsI3juFtJa3ahCcFqKqDSa4AdgObgB1VdWeSq4A9VbUL+I/A04DfTQLwuap69eRaNTarzrpNmzdPuk2SljJjE6BMwVHzbjjrTj/99PVok6TFzLp54yhqVzTcsLvvbjvAM3K8f35pKmriD+muqhuBGxdte+vQ65dPtAGTs+qsO+7008w6aSomn3Uz7qh5N5x1W7Zs8Y8tTYVZN2ccsx/PDzdM8mQGww13jeG6kibAB3SvmVkn9cgoWWfemXdSX5h1AyP31C433HDklknSBmLWSWqFeSepb8Yx/HjJ4YZH8m0v+BZu2vO743hrqVlJblvTiQ5TWbPVZt13nvQc9vz4T02wRdLsy0/8tFk3BavNO0lTYtYBYypqJfVEQXyemaRZZ9ZJaoFZN8+iVmqN3+hJaoFZJ6kFZh0wnomiJEmSJEmaCntqpdb4hZ6kFph1klpg1gEWtVJz4jAVSQ0w6yS1wKwbsKiVWmP4SWqBWSepBWYdYFErtaUAZ8mTNOvMOkktMOvmOVGUJEmSJKm37KmVGhLKey8kzTyzTlILzLoFFrVSaww/SS0w6yS1wKwDLGql9hh+klpg1klqgVkHWNRKbXFCAUktMOsktcCsm+dEUZIkSZKk3rKolRqTqjUvktQXo2TdqHmX5PgkNyW5t/u5eZnjtnXH3JtkW7ft6UluH1oeTfKr3b43JHlkaN9lIzVUUu+ZdQMOP5ZaY3EqqQXTzborgZur6p1JruzW/8PwAUmOB94GbGEwiPC2JLuq6jHgnKHjbgN+b+jUD1bVFZP+BST1hFkH2FMrNaYG4bfWRZJ6YcSsGz3vtgI7u9c7gdcsccwrgJuq6kD34e4m4MLhA5KcBZwE/NmoDZI0i8y6ORa1kiRJT3RCkj1Dy+WrOPc5VbUfoPt50hLHnALcP7S+r9s27BIGvRXDnzx/MMmnk1yf5LRVtEmSljITWefwY6klhT2ukmbfeLLu0arastzOJB8BvmmJXW9Z4fWzxLbFjb4YeP3Q+h8CH6iqryV5I4OekZet8P0kzRqzbp5FrdQap36X1IIJZ11VvXy5fUkeSnJyVe1PcjLw8BKH7QNeOrR+KvDxoWt8F3BMVd029J6fHzr+auBda2u9pJlh1gEOP5aa4+zHklowzRlBgV3Atu71NuCGJY7ZDVyQZHM3Y+gF3bY5lwAfeNzvNPjQOOfVwN2jNlRSv5l1A/bUSq2xOJXUgulm3TuBDyW5FPgccBFAki3AG6vqsqo6kOTngVu7c66qqgND1/g3wKsWXffHkrwaOAgcAN4wwd9BUh+YdYBFrSRJ0lh1Q+fOW2L7HuCyofUdwI5lrvEtS2x7M/Dm8bVUktZuI2WdRa3UkgIO21MracaZdZJaYNbNs6iVmuLzZiW1wKyT1AKzbo5FrdQaw09SC8w6SS0w6wCLWqk9hp+kFph1klpg1gE+0keSJEmS1GP21EotcUIBSS0w6yS1wKybZ1ErNaWgDk+7EZI0YWadpBaYdXMsaqXWeO+FpBaYdZJaYNYB3lMrSZIkSeoxe2qllnjvhaQWmHWSWmDWzbOolVrjMBVJLTDrJLXArANGHH6c5O1JHkhye7e8alwNkzQhVWtfGmXWST00StaZd+ad1BdmHTCentp3V9UvjeE6kiZutgJsnZl1Um+YdSMy76ReMOvmOFGUJEmSJKm3xlHUXpHk00l2JNk8hutJmpQCDh9e+9I2s07qi1Gzzrwz76Q+MOvmHbWoTfKRJHcssWwFfhP4VuAcYD/wy0e4zuVJ9iTZ88gjj4ztF5C0St53sSSzTpox3me2rHHknVknbRBmHbCCe2qr6uUruVCSq4E/OsJ1tgPbAbZs2TI7f0Gpb2YowMbJrJNmjFm3rHHknVknbRBmHTDiRFFJTq6q/d3qDwB3jN4kSZNTPs9sDcw6qW/MurUy76Q+MevmjDr78S8mOYfBiO7PAj86coskaeMx6yS1wryT1DsjFbVV9fpxNUTSOiiomp1JAdaLWSf1jFm3Zuad1CNm3bxxPKdWUp84TEVSC8w6SS0w6wCLWqk9TiggqQVmnaQWmHXAeJ5TK0mSJEnSVNhTK7WkaqYetC1JSzLrJLXArJtnUSu1xmEqklpg1klqgVkHWNRKzSm/0ZPUALNOUgvMugGLWqkp5Td6khpg1klqgVk3x4miJEmSJEm9ZU+t1JLC55lJmn1mnaQWmHXzLGql1pT3XkhqgFknqQVmHWBRKzWlgPIbPUkzzqyT1AKzboH31EotqRp8o7fWZQWSXJjkniR7k1y5xP7jknyw239LkjPG/FtKat2oWTdiz0eS45PclOTe7ufmZY77kyRfSPJHi7af2eXjvV1ePrnbbn5KWmDWzbOolTQ2STYB7wFeCZwNXJLk7EWHXQo8VlXPBd4NvGt9WylJE3clcHNVnQXc3K0v5T8Cr19i+7uAd3fnP8YgN8H8lLSxbJiss6iVGlOHa83LCpwL7K2q+6rq68B1wNZFx2wFdnavrwfOS5Kx/YKSxGhZN4bhfMM5txN4zZJtrLoZ+NLwti4PX8YgHxefb35KehyzbsB7aqXWTHZCgVOA+4fW9wEvXO6YqjqY5IvAs4FHJ9kwSY2Z7uQpz6mq/QBVtT/JSas499nAF6rqYLe+j0FugvkpaTGzDphSUXvbbbc9muQfhjadQD8C2XaOl+0czb9Y7Qlf4rHdH6nrTxjhPZ+SZM/Q+vaq2j60vtS3aIu/BlzJMTPBrJs42zleG7Wd08g6OEreJfkI8E1LnPeWEd/3SBm5IfPTrJs42zleG7WdZt0IWTeVoraqThxeT7KnqrZMoy2rYTvHy3auv6q6cMJvsQ84bWj9VODBZY7Zl+QY4JnAgQm3ayrMusmynePVl3auxDpkHVX18uX2JXkoycldz8XJwMOruPSjwLOSHNP1YAzn6IbMT7NusmznePWlnSth1i3wnlpJ43QrcFY3m92TgYuBXYuO2QVs616/FvhoVU29p0GSxmg457YBN6z0xC4PP8YgHxefb35K2kg2TNZZ1Eoam+6btiuA3cDdwIeq6s4kVyV5dXfYNcCzk+wFfpLlZ8qTpL56J3B+knuB87t1kmxJ8ltzByX5M+B3GUyCsi/JK7pd/wH4yS4nn80gN8H8lLSxbJis2ygTRW0/+iEbgu0cL9s5g6rqRuDGRdveOvT6q8BF692uDaIv/5Zs53jZzsZU1eeB85bYvge4bGj9e5Y5/z4Gs8kv3t6X/OzLvyXbOV62szEbKeviqBVJkiRJUl85/FiSJEmS1FtTL2qTXJjkniR7k2zYe0OSfDbJZ5Lcvmja66lKsiPJw0nuGNp2fJKbktzb/dw8zTZ2bVqqnW9P8kD3N709yaum2cauTacl+ViSu5PcmeTHu+0b7m+qfjHrRmPWjZdZp0kx60Zj1o2XWdeOqRa1STYB7wFeCZwNXJLk7Gm26Si+t6rO2WDTgF8LLJ7O+0rg5qo6C7iZjTGRxLU8sZ0A7+7+pud092JO20Hgp6rq24EXAW/q/k1uxL+pesKsG4trMevGyazT2Jl1Y3EtZt04mXWNmHZP7bnA3qq6r6q+DlwHbJ1ym3qlqj7BE5/btBXY2b3eCbxmXRu1hGXaueFU1f6q+mT3+ksMZvA9hQ34N1WvmHUjMuvGy6zThJh1IzLrxsusa8e0i9pTgPuH1vd12zaiAv40yW1JLp92Y47iOVW1Hwb/mYGTptyeI7kiyae7YSwbauhHkjOA5wO30K+/qTYes24y+vT/0qxTC8y6yejT/0uzTlMx7aI2S2zbqNMxv7iqvpvBkJo3JXnJtBs0A34T+FbgHGA/8MvTbc6CJE8DPgz8RFX947Tbo94z69pm1qkVZl3bzDpNzbSL2n3AaUPrpwIPTqktR1RVD3Y/HwZ+nyWeqbSBPJTkZIDu58NTbs+SquqhqjpUVYeBq9kgf9MkxzIIvvdX1e91m3vxN9WGZdZNRi/+X5p1aohZNxm9+H9p1mmapl3U3gqcleTMJE8GLgZ2TblNT5DkqUmePvcauAC448hnTdUuYFv3ehtwwxTbsqy5MOn8ABvgb5okwDXA3VX1K0O7evE31YZl1k1GL/5fmnVqiFk3Gb34f2nWaZpSNd1RId10378KbAJ2VNUvTLVBS0jyLQy+xQM4BvidjdLOJB8AXgqcADwEvA34A+BDwOnA54CLqmqqN/Mv086XMhiiUsBngR+du79hWpL8r8CfAZ8BDnebf4bB/Rcb6m+qfjHrRmPWjZdZp0kx60Zj1o2XWdeOqRe1kiRJkiSt1bSHH0uSJEmStGYWtZIkSZKk3rKolSRJkiT1lkWtJEmSJKm3LGolSZIkSb1lUStJkiRJ6i2LWkmSJElSb1nUSpIkSZJ66/8HT0NsxKftztUAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1152x432 with 6 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"init()\n",
"plot()\n",
"print(dh)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"if 'is_test_run' in globals():\n",
" time_loop(2)\n",
" assert np.isfinite(dh.max('phi'))\n",
" assert np.isfinite(dh.max('T'))\n",
" assert np.isfinite(dh.max('phidelta'))\n",
"else:\n",
" vtk_writer = dh.create_vtk_writer('dentritic_growth_large', ['phi'])\n",
" last = perf_counter()\n",
" for i in range(300):\n",
" time_loop(100)\n",
" vtk_writer(i)\n",
" print(\"Step \", i, perf_counter() - last, dh.max('phi'))\n",
" last = perf_counter()"
]
}
],
"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.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}