Newer
Older
{
"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": [
"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, default_target=target)\n",
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
"φ_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": [
"$\\displaystyle \\bar{\\epsilon} \\left(δ \\left(\\frac{{\\partial_{0} {{φ}_{(0,0,0)}}}^{4}}{\\left({\\partial_{0} {{φ}_{(0,0,0)}}}^{2} + {\\partial_{1} {{φ}_{(0,0,0)}}}^{2} + {\\partial_{2} {{φ}_{(0,0,0)}}}^{2}\\right)^{2}} + \\frac{{\\partial_{1} {{φ}_{(0,0,0)}}}^{4}}{\\left({\\partial_{0} {{φ}_{(0,0,0)}}}^{2} + {\\partial_{1} {{φ}_{(0,0,0)}}}^{2} + {\\partial_{2} {{φ}_{(0,0,0)}}}^{2}\\right)^{2}} + \\frac{{\\partial_{2} {{φ}_{(0,0,0)}}}^{4}}{\\left({\\partial_{0} {{φ}_{(0,0,0)}}}^{2} + {\\partial_{1} {{φ}_{(0,0,0)}}}^{2} + {\\partial_{2} {{φ}_{(0,0,0)}}}^{2}\\right)^{2}}\\right) + 1\\right)$"
" ⎛ ⎛ 4 \n",
" ⎜ ⎜ D(φ[0,0,0]) \n",
"\\bar{\\epsilon}⋅⎜δ⋅⎜───────────────────────────────────────────── + ───────────\n",
" ⎜ ⎜ 2 \n",
" ⎜ ⎜⎛ 2 2 2⎞ ⎛ \n",
" ⎝ ⎝⎝D(φ[0,0,0]) + D(φ[0,0,0]) + D(φ[0,0,0]) ⎠ ⎝D(φ[0,0,0]\n",
" 4 4 \n",
" D(φ[0,0,0]) D(φ[0,0,0]) \n",
"────────────────────────────────── + ─────────────────────────────────────────\n",
" 2 \n",
" 2 2 2⎞ ⎛ 2 2 \n",
") + D(φ[0,0,0]) + D(φ[0,0,0]) ⎠ ⎝D(φ[0,0,0]) + D(φ[0,0,0]) + D(φ[0,0,0]\n",
"\n",
" ⎞ ⎞\n",
" ⎟ ⎟\n",
"────⎟ + 1⎟\n",
" 2⎟ ⎟\n",
" 2⎞ ⎟ ⎟\n",
") ⎠ ⎠ ⎠"
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
]
},
"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": "iVBORw0KGgoAAAANSUhEUgAAA8wAAAAaCAYAAABip9XBAAAABHNCSVQICAgIfAhkiAAAFZFJREFUeJztnXncXUV5x78JCRBWQYoRgoQEMCxCCAiUUrggSyOFRsWmtWxFalk+hRgVClrzghWksqSAAqHSC0oLBSogSAGRsglBliiRsgh5I2EJW8KWEAi8/eM3wz3vuefce2bOnHvywnw/n/dzk7PMM/Ocmec8M/PMnGFEIpFIpA6mAR/JON7X43xEIpFIJBKJRDz5EnAv8CrwLjCx3uxEIpHIB4Z+YCDjL1KOQ5Ae/6HujEQikSHDKODbwKPAW8DTwKnAyDozFYlEKqcPeAd4HrgG2Nw1gS2B94DXgYuB7wCjw+XPmdOBW5ERWwq8AjwEzAA+6pjWgcC5wJ3Aa8i5+oljGgfTcnCPCCSjn2wHegA9yCyGAYejgY3XgSVIL8cCK6Wu/ajJ60+B3yM9vgrcBXwZGN4lf38KXA08BywzvzcDnw1UFtdn7CMDYD+T7wVGzlPAlcAf51zvomOAwzrky/69W1JGCF4skM/k35EV5cOVEO03yRhk455F9bofmAmsUyqX9RO6XC7phX5GLpxj5O0SME0X2zcUCVVX6rBjvcS3XtdlY/ZG77b56FksA+YBF6EOYihClS+kn+fCx4G5wNvAFSYfv0PP98IK5XYihE7L+nxDgdDvGlf/cKgR0ha5puUj2/UeHxvSAL4H/A+qP/dlXTSsQyaPAc4DjkMOSN28DTwIPAK8AKwO7AzsgBS5M1JQEeYA2wJvoEYxAbgMOKjg/RsBDyMHYA3g74B/CyCjH4Vozsw49wZwRsbxS1Hn/QXgZ8CbwF5owONq4Iu0Zq2OBM5Hzt5twB+AjwGfB9bOuD7Jt9CgyUvA9SaN9YDtTFrHByiL6zP2kXG6yevLaCTpJWBT4ABgBJqdShtbFx2DIjGmZMgGOd57AjcAf15CRllWB76ROjYC+CZ6Dqdl3NNEOq+bsu03yXjgV8D6wLVodmFHYA/gMeBPUF0ZaoQul2t6IZ+RKxsCayIHcXmA9Fxt31AjZF3ptR3rNT71ui4bcwiy2c8AP0cDpCOBscCuwCdoH7j1IWT5Qvp5RVkZ5X8CsC9wtzm+Buo0j0E2pdMgfGhC6bSMzzdUCPmu8fEPhxIh26prWj6yfe4pa0MeQO/2tVCdKsQ/oUa0d9EbKmbVnOPfRfn8oUNaewCboQGDBm4jUsOAXwBPAt8nf4bZR0Y/bh2SKSbdp5ADZxmJRhQH0GynZU9gf9pHFUcjQzoAfCFDjjWotyBHNE1WyFI/7p0r12fsKmM0chCeRw0wyR60dJnEVcfduMfcc0CFMnzZ1si6vweyylCm/aa5iezw3bPM8Qs8062b0OVyTS/kM6oTH9s31AhVV1YUO1YlPvW6Lhsz16S/Xsa5VQLKCVm+kH5eUU4yaR+TcW6mOffFCuR2IpROfX2+oUSod42PfzjUCNlWXdPyke1zT1kbcpm5bqMu1w2iz9zUyDm/H8VDOatc+2yd/Fs872/g1sCOQ6Hqu9HSUVaH2UdGP24dwEvJN/Rbm3MPFEzLvjTOTR0fjozEm8AfOeStn3CzkXnP2FXGTiada3POv4bCCJOE1LG9fgGDwxNDyijDl6k3BM2HBv4vyHHm3nm0OxRrotHFN9EIZWgmIMP9BAqVfA2Nnl5BeWc2dLnKptegdx3m7YysWQHS8rV9IRiJNoWbg8LKFgBno9mw1YCF6MVelpB1pS479nm0/mw2sHHONScb+V8JKLdB93pdp4251cj+JZr1H1+BjF6Vr6yfl8coYBGaecoa/LITR73cD6FXOs3z+UKwhUn7v1PHd0IRD/Oppl/QwP9d4+MflmU91J/o1odaRn5HsCgh65VrWj6yQ7eDojakaa4bmz5RZv3CqyjUAxS2cnLi7wpz/C4UN/7bLhk7rEQ+9je/eTJCsgWKc/9X4I6KZKyCQklOQp3zPchf+2XXlGeNetljk8jeiTfNO+Y3HcK4C7AJCulahAZKTjB567amw6Usnej0jF1kPIHCNXakfdR9N9QIf5E6HlLHf29+f8TgULiQMsqwvfmtY4a5Se9noPY0vzejl1aS11Fo3moohCckDbSm5nDgN8ieNNGz3ha9HPNo0l1PoctVl558mGR+H+pwTZNida2M7SvDukinZ6PQsXPQ+3YaWgt8pLlmRoc0mhQrY8hnW5cdW4icoB3JDo8fZ47/mvalU1Ct7amz7UwHHkeh8dug5QRPAX9b4N4mva8/najKz/scqo//QcsHSmI7KW8HkNVkxdJpns8Xgiw7fChwO1oqswMaDMyjSe/9AR//sCxrAKfQ6jtdYo4/wOA+1XS0EV0WTXpfr1zT8pEduh2UtiEjOpyzL7WlOefvQqPuR6FY/77EueOBqehB/sg3czl8HVWytVGj2xUp4HuB5aQZAfwYhbGcVKGc0UZOknnoJXd76vhL5neTjHTGJf49AW3CkoddmwFa9J7k0+Z3IVoX8KnU+TvQpgsvZqTrUpYkLs/YRcYryOE9C61vuAatfxiPQqRvodWptYTS8SjUsX+PdqctlIyy2A5zL2azVwQ+aX4fzzn/BLAP2jHx1oByv4tmM3ZEbSo0octVl558sI5aCL2WsX1luNzITu4f8n00yzwZrQduIsezLCGfbV127G7U8XkZvS/SnINm5o+m3fGqmrrazoFo3f1uqP4CnIiiAC5G9idEFEZV5euVn7ef+d2Q7M/57WV+Q6+b7kQv6kwnny8ESTu8EtpPZhqaZDuO7MGJuvHxD8vSz+B6dwQaWLic7D14yhCyXrmm5SO7bH59bYgdmFi7y3XvMww5ze/ReWdsG+Z0YOr45eb49m13DObj6OVZOGNofUEyVOFGtImBLw2KhXCcgmYFkzMLfYQNyZ6BRlU+hkZOtkYx+u+hsM1tU9d/yaT7ezTjYBmBNnOwOprcRe4Z5robMs6dZs4tRxX0M6gSbkVrR7n/DVCWJEWfsa+MKcg4JmU8gfSZJpSODzXXXV+hjDKMQINjy5Bz2Wt8bAGUC8GaRef2a9e8nOiRdiceQy9inxCrInoKXa6y6TXoXUj2PchWddoBuGhd87V9ZdjLpHsH7ZtyPmrOvUX39VVFyxiyrtRtx+5DYZNJDjAyz2+//H2qtD112JhdUIdk/4xzY4y833RJo476kyS0n5fH/JScvL+xAWTVrdMknXy+ENxm0t8SdTSX0d1PTlKHP2Bx8Q9Dc56Rt1e3CxPUUa9c0/KRXTa/vjbkH831x6VPpEOy/ww4E714JpkMd9oZ0K5BSBvfScjJmNslY88hJ+DVLtclGY0cidFo7dI4FPYxqdNNJdkRzSqfiRyyqjgZrTdaiDp8c1H43VnIAexLXX85qgTj0YjYLLRJxRz0uZMnzHWddsI8Fvgaeg4HZ5y34c3D0MDIrWjtwO/QqP4CYHfaQxRdy5Kk6DP2kXE8cBWaoRmP1j9sj0LVLgP+JXV9CB1Da+1c1vrgUDLKsBXqwD1MmPAzV3xsQdXYDstA4HSnI/v4ILIpfWj39CKE0FPoclWlJ1dWQqGn/0d+ZBQU16Gv7SuDtcF2s6EkduT7QrrPeIVqTy7Ptm479igKmxxj/j/KyH+JzlFhddqeKtrOD9GAXFZnyEZDdHMc66g/SXrh562Odgqfa2Sl/9ZEAw8LaN8r5WgUyfYWmlwqYr/r1qmlm88XgonIj7gOdf7mkL0cIo+62qSrfxgaO8nTbUAryYpSr8qk5SO72z2+NuTHaCDtDOC/UMd8LGR3mKej6eu5dN8MYCJyIJKhYWugbdgfpfNavLIsRDtv7oO+rXVpRXJsKPbjaAOIOrA7we2WOv4eGkH/OhrYOBiti1yAwg/sdusv5KR7DFo/+Qha+/tKxjWLzO9TtDfipWgXO9CgQhHyypKF7zPOk9FAnw24DtXzp1BH+0HkAD+DXiTJ8MGyOgaNsu5i7vl5xvkQMspiQxk/LOHY0HrB5I3MrpW6LgTDkLM6H+2JMB1FSmwQUEboctWhJx8moEiTUGHuoW1fEXZHTnpemOQS4NSA8kI+27rt2KPmd0vzeyIKDz+B1rPsNb1uO9sgp/smssPPbbj8s4HkVV2+Kv28Dc1vni72RaHr6YGHqchvOhVtMng3Gij6RKB8VanTIj5fWcahJZ0ro7b+MLKRWREPKxIN3P3DkAxD7fdZwi/zgbD1yjUtH9mh8utqQ55B0WXvoN3xTyKnwzzNZG4K6vTeRP4GSusgA/FbBvfwt0MPvtOC/pDMR41/K7I/n1CWNVCM/BZoNDE5xW83XbnI/D/re8AhsA5G1m5wy9Es1UQ0or4WGvh4xBxbimZE0kxD4R9zkeHMiyR4zPwuzjlvHZFOIZBJOpUlD9dnnCfDfvf4tox7lqDIiuGoDifx1bElb7OvkDLKUnTDrzFopPUVVCeupv0TDJ9EM/9voc+vTUaDZ58JldlA2Lq9ec75zcxv3hoaH85B9eB+tIRgVWQvr+h0kyOhy1WHnnwIuX4Zwtu+boxC79T5yB4lGYcGBGbTWpMagtDPtk47luwwj0ezRfcC/16RvCL0uu1sY37/kHN+ivm9K5C8XpWvCj/PLj3Km9ixm6NdnDo+Hc1AXoSiWY5Fs3xHBcpXVTot6vOVxdrhn6JO6DfN/0+l3EbDVePrH4ZiE2QvXWaXXQhZr1zT8pEduh0UtSH7oUm3B5GPtjIFll79BHUCt845b79Llv6m1XHm+Ne6CQjIQiNzHY97G3Re8zAKhZJk/T1o7r3T/H+qp4xu7Gvuf8Thnq+Ye5oZ504w5x6i+8tnPTTSspjsta03mrT+qmC+fMoCbs84T8a55vgpOffdac4XHQntpGPLqqhz+S6O33VzkBGC2UZOp3CVceg5nIYGkCaijdWuSlyzGRrxOwsNuu2HnIkBqhnQauDftsabe+eR/9mCJYT75Mv6qB5UsclKktDlKpteg96sYT7byCka3t6N0LavG+ua9B7LOHetOXdzIFmWXrWBXtixLY2MWWiviOVU59xCsXrdaxtj9ZwVHbg2rW/vhtJLL8tXxs/LYj2T3n0Z53ZGM/TpiLCVUb1Kf5f5B3TeyNSFKnTq4vOV5VQjK/l953vMsUMy7whHA/93TWj/0JXPmvSrCvsOWa9c0/KRXUU7KGJDzsR9Hfn7C6rzQme/as6nv2toZ1v3KSCj6GL1CWRvPjY8kc+7M86PN/dmfV/P0sC/gfXReVG6i4ytGLxRimVjtPZrgOx1WGtlHPs06qS9Tnv4iP2u4P058rKwgyf/nDq+N3qpLGbwp0J8yuL6jH1k/KU5/jytcCzLZFOWpShsI4mrjpMcbGT+rMM1vjKaJu3DuqTdjaIbft1MezjovgwOd7yJ9pCXS1BYZieq3OSjkx24ydyf/s7mWeb4BanjTfx1br9F+yuyI3eKzFQW1ZNruSCsnpI0KGZjm5Srz7ejNrxml+tc6pqr7QP/cgxDbf1dWjOFoJkrG9WU5dxn4VLGkHWlTjtmOzOLyO80ZlH1BkO9tDG7mnufY/D7cVW06+8AxdaS1lF/fP28Jv76eoT29rYxWrO8iPbNvjYg2y/+NtkDXUmq1Gkn2+3q8zUp1x7thojJtm43M5xHsU1F6/AHfPzDJmFsF2jgdQD4luN9ddUr17R87ITrPb42JMnF5rpNu1w3iD5zUyPn/CXm/E6p43Y3zK/S2nwjjybFKts0NNJ/Kxo9Pg0V6klaL4ctM+7rN+fHpo5PMbKbtBr3k4ljRbdz7yO/w+wqow+FsN6IZu1PRzN3S829N5BtaGajcIHzkF6uQ07Dm6gzk+RQk9ZyNBvTl/F3WIaM9Wl1Qu8web/SpGPj/MuWxfUZ+8gYjnZsHEC7qV5Ca82K/Xh82854uOk4TdFRSR8Zl5q0D+qSdjfsB907hWNvbK5Zgkb27N9SWiGiG5E9e3ER3XfjbFL8xePatvrJtgOgF4IddbwG6f6XtGb60oMnZXQ+0qRp0z4PGfELUYhkOvQviybF9ORaLgirJx8bW0a3w1BkQzenFdzqmqvtg3LlsLMcL6BZq6uRbbqG1s6zF9D65FUeTYqXMWRdqdOOQetZLaT4956bVGd7oLc2Blr1ZD6aLZlJ63ldRbEOS5Pe1x9fP6+MvuzO7i8in2gWGgRbhPYdSWM7zOkolhm0lgTk0aQ6nfaT3R59fL6y9W8h0mF6l39bL7N8rDRNeu8P+PiHIW3XJJPWS6hD+NcF72vS+3rlk5aPnXC9x9eGJGl2KHMuM+jcYZ6DRuZWSx2fih74ErT5R5GMHdbluq2R8zDHpL0cOUe/Rg0+b9Ssn+yC95njeX/9XfKTTierw+wqY3fgP5HRXYwe+ouoAR9Cu/GxfANt1LQYzRDOQw7VWI88DZAfq78uasTz0O6HL6MwwayPhvuUxfUZ++prJGpU9yKjuBw5p9eTHxXhouMkWyCdPk3+XgBlZDxkylA2RO1wk8+sHbwtB5i8bZrxZ0PN/wLpM+2QzUYdw040KW70+3BrW/10NoAboXWOz6G6PR9tjJJlV8rqfAwy5LYdvYmM+ZUUCyVuUlxPLuWCsHrqw93GltHtpibdywtc26S4DsHN9kG5cqyKHLynkU17AXV67He77ael9u6SThO3MoaqK3XaMVAkj0u5oVrbY+mljVkL1Zl+I+sVFB2U/vxnJ5r0vv74+nll9XUoWlf/FoqEmkX7DKOlTEh2k+p02o+fr5vl85XR54Y5aYIGIAaQTesWBdSkHn/A1T8MabtA672fQx30Mwve06T39conLZ/rXe/xtSFJmqjMGxe49n3st6g+53JTJBLpCR9BA1ZVf+bAMhkZn04vuv2RoU+uKdkd2ZGsmbihRq91/mGirG7/BtWz6cFy5EesI26E1tddJr1Qa4J7Taw/btShr9moU53kcTSbNdSJ9a84UVcfTOx+IU7fff+Cuel61MNfkXe3i0Q+bOyPRsSz1mtUwTpoBv9aFHI9Hs1y/YDWbP4GJk/novVLU5AjMUD+TodDiV7r/MNEWd3a9ZkTg+XIj1hH3AipL7sG3HVDyRWJWH/cqENfU9FM1xEoimwmWp7kNCO1ghLrX3Girj5YrIKiuF41f0593lHo+2nJ0Ia6nZFIJFIfO6D1I4uRYzoHRaIkOQjtxPoG+qzENPPvOOAWCc1w4Du0wnDvqDc7kZrZHNWDy+rOSOQDz9EoXHUZWoKQtzluJBJZ8eljcF93hk8iq6BQzOloF8A4khKJRFw4Be0KHYmEZgJaJvA02mQq1DqyyNBkKnJ2uu2dEolEIpGIpYF2Jz+KODEciURq4hrg/LozEYlEIpFIJBKJ+BDDJCORSJV8CoVuRyKRSCQSiUQikUgkEolEIpFIJBKJRCKRSCQSiUQikUgkEolEIpFIJBKJRCKRSCQSiUQikUgkEolEIpFIJPKh5v8Bsb8cqcFJbMQAAAAASUVORK5CYII=\n",
"$\\displaystyle \\left\\{ \\pi : 3.14159265358979, \\ T_{eq} : 1.0, \\ \\bar{\\epsilon} : 0.01, \\ j : 6, \\ α : 0.9, \\ γ : 10, \\ δ : 0.3, \\ θ_{0} : 0.2, \\ κ : 1.8, \\ τ : 0.0003\\right\\}$"
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
],
"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": [
"$\\displaystyle \\left[ {{T}_{(0,0,0)}} \\leftarrow 0.0111111111111111 {{T}_{(-1,0,0)}} + 0.0111111111111111 {{T}_{(0,-1,0)}} + 0.0111111111111111 {{T}_{(0,0,-1)}} + 0.933333333333333 {{T}_{(0,0,0)}} + 0.0111111111111111 {{T}_{(0,0,1)}} + 0.0111111111111111 {{T}_{(0,1,0)}} + 0.0111111111111111 {{T}_{(1,0,0)}} + 1.8 \\cdot 10^{-5} {{φ_D}_{(0,0,0)}}\\right]$"
],
"text/plain": [
"[T_C := 0.0111111111111111⋅T_W + 0.0111111111111111⋅T_S + 0.0111111111111111⋅T\n",
"_B + 0.933333333333333⋅T_C + 0.0111111111111111⋅T_T + 0.0111111111111111⋅T_N +\n",
" 0.0111111111111111⋅T_E + 1.8e-5⋅phidelta_C]"
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
]
},
"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",
"\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7kAAAF1CAYAAAAtEi0mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3df7RlZX3n+ffHQtCoUX6ooSnSkLacDjHT/qhBZpzYtgiWrsSyO9pijJY9MCx7ZJKM6VnBcaIJmrUk3Wk7ztB2oxLRNqLBpK2kK00jyiSdpaSKiCjQhJKk5QJLxEKippFU1Xf+2Pvee7icW/cW55577t3P+7XWXvecvZ9zznMO5cfzPc+zn52qQpIkSZKkIXjcrDsgSZIkSdJasciVJEmSJA2GRa4kSZIkaTAsciVJkiRJg2GRK0mSJEkaDItcSZIkSdJgWORKkiRJkgbDIleS1LQk3x3ZDif5byP33zDr/kmSpKNjkStJalpVPXl+A74O/NTIvo/Pun+SNIkkz09yZ5KvJ3nbrPsjrQeLXB1Rkv8pyZeTvDjJN5PcmOSMWfdLkiRJq/LnwHOBncA7kvxEkpOT/Mckr0pya5K7kvyjGfdTWjMWuVpWkh8Afge4DHgD8HHgKuDjSTLLvkmSJGllVfXdqvqrqvoS3Xe5lwP/Gvgm8FTgYeA1wJVJTpldT6W1Y5GrIzkTCPBB4EnAt4H/B/hx4Idn2C9JkiStQpLzk9yQ5ErgEPAE4CeBXwOOA75dVTcAXwZeNrueSmvHIldH8kzg7qqq+R1V9RDwAPBDM+uVJEmSVpTkx4H3AK8G/i1wAd0I7jHAXUua34vf7zQQFrk6kruBraNTk5M8ETgemJtZryRJkrQaLwN+v6ruBfYC36M79exhHj0r7xT8fqeBsMjVkdxAF4ZvpZu2vAX4VeBPquruWXZMkiRJK/oruu9vAO8E/qiq/itdofvLwLEASXYCfxfYM4tOSmvtmFl3QBtXVf1NH3pXAM+nO4/jPwNvnGnHJEmStBofB/5hkq8CdwD/S7//5+nWXPl14InA04F/VFUPzKSX0hrLyOmW0rKS/Dtgf1X9yqz7IkmSpMkluQD42ap6yaz7Iq0lR3K1KlX1s7PugyRJkiStxHNyJU1VkiuS3NdPlRp3PEnen2R/kpuTPH+9+yhJk0qyI8ntfZZdPOb4i5P8WZKDSV6z5NiuJHf0266R/S9I8pX+Od/vNeolzdpmyTqLXEnT9hFgxxGOvwLY1m8XAh9Yhz5J0ppJsgW4jC7PzgBen+SMJc2+DrwZ+O0ljz0BeBfwQrrr078ryfH94Q/Q5eJ8Rh4pS6WjVlUfcqqyVmszZZ1FrqSpqqo/Ag4coclO4KPV+SLwtCQnr0/vJGlNnEm3bsWdVfUw3cq1O0cbVNVfVtXNwOElj305cG1VHegX/bkW2NHn4A9W1Rf669V/lO5ap5I0K5sm6yxyJc3aKTzygvRz/T5J2iwmybHlHrv0mqVmo6RZ2zRZtyEWnjrppJPqtNNOm3U3pE3txhtvvL+qnn40j3n5P3hSfevAocle9+bv3wI8NLLr8qq6/CieYtx5F4Nc9t2skya3QbNukhxb7rGbNhvNOmlyZt1kNkSRe9ppp7Fv375Zd0Pa1JL816N9zP0HDnHDNVsnet3Hn/y1h6pq+wRPMQecOnJ/K3DPRJ3aoMw6aXIbNOsmybE54CVLHnt9v3/rkv2bIhvNOmlyZt1knK4sadZ2A2/qV1k+C3iwqu6ddack6SjsBbYlOT3JscB5dNm2GtcA5yY5vl+E5Vzgmj4Hv5PkrH6l0TcBn5lG5yVplTZN1m2IkVxJs1IcqqXrAqytJJ+g++XupCRzdCvrPR6gqv4NsAd4JbAf+Gvgn0y1Q5IaNN2sq6qDSS6i+xK3Bbiiqm5Jcgmwr6p2J/kfgN8Djgd+KsmvVtWPVdWBJO+m+/IIcElVzS/W90/pVqh/IvCH/SZJyzDr5lnkSg0r4PCUT/GqqtevcLyAt061E5Katk5Zt4fuR7vRfe8cub2XR07JG213BXDFmP37gOesbU8lDZVZt8giV2rc4Uet8C5Jw2PWSWqBWdfxnFxJkiRJ0mA4kis1rCgO1aa4IoUkPWZmnaQWmHWLLHKlxk373A1J2gjMOkktMOs6FrlSwwo4ZBhKGjizTlILzLpFFrlS4/zFT1ILzDpJLTDrOi48JUmSJEkaDEdypYYVuECBpMEz6yS1wKxbZJErNc6rqUlqgVknqQVmXcciV2pYUS5QIGnwzDpJLTDrFlnkSi0rOGQWSho6s05SC8y6BS48JUmSJEkaDEdypYYVnrshafjMOkktMOsWWeRKTQuHyKw7IUlTZtZJaoFZN88iV2pYAYc9d0PSwJl1klpg1i3ynFxJkiRJ0mA4kis1zmktklpg1klqgVnXsciVGlYYhpKGz6yT1AKzbpFFrtS4w2UYSho+s05SC8y6jkWu1DB/8ZPUArNOUgvMukUuPCVJkiRJGgxHcqWGFeGQv3VJGjizTlILzLpFFrlS4zx3Q1ILzDpJLTDrOha5UsM8d0NSC8w6SS0w6xZZ5EpNC4fKaS2Shs6sk9QCs26en4IkSZIkaTAcyZUaVsBhf+uSNHBmnaQWmHWLLHKlxnnuhqQWmHWSWmDWdSxypYZVee6GpOEz6yS1wKxb5KcgSZIkSRoMR3Klxh12WoukBph1klpg1nUscqWGdddTc0KHpGEz6yS1wKxb5KcgNa07d2OSTZI2vulnXZIdSW5Psj/JxWOOH5fkk/3xG5Kc1u9/Q5KbRrbDSZ7bH7u+f875Y89Y4w9G0qCYdfMcyZUa5lLzklow7axLsgW4DDgHmAP2JtldVbeONDsfeKCqnpXkPOBS4HVV9XHg4/3z/Djwmaq6aeRxb6iqfVPrvKTBMOsW+e1WkiRpMmcC+6vqzqp6GLgK2LmkzU7gyv721cDZSZaePPd64BNT7akkPXabJuscyZUad6hcoEDS8K1B1p2UZHSU4fKqury/fQpw18ixOeCFSx6/0KaqDiZ5EDgRuH+kzet49BfG30pyCPg08J6qqsnehqQhM+s6FrlSw4q4QIGkwVujrLu/qrYvc2zct8qlX9CO2CbJC4G/rqqvjhx/Q1XdneQpdF/83gh89Cj6LKkhZt0iv91KjTtcj5tok6TNYMpZNwecOnJ/K3DPcm2SHAM8FTgwcvw8lkzfq6q7+7/fAX6bbqqgJC3LrOs4kis1zKXmJbVgHbJuL7AtyenA3XRf4n5mSZvdwC7gC8BrgM/NT8dL8jjgtcCL5xv3Xw6fVlX3J3k88JPAZ6f5JiRtbmbdIotcSZKkCfTnnV0EXANsAa6oqluSXALsq6rdwIeBjyXZTzeqcd7IU7wYmKuqO0f2HQdc03/p20L3pe+D6/B2JGmszZR1FrlSw4q48JSkwVuPrKuqPcCeJfveOXL7IboRjHGPvR44a8m+7wEvWPOOShoss26RRa7UOK+TK6kFZp2kFph1HYtcqWFVcMjFoyQNnFknqQVm3SI/BUmSJEnSYDiSKzUtHB57OTNJGhKzTlILzLp5FrlSwwqntUgaPrNOUgvMukUWuVLjvE6upBaYdZJaYNZ1LHKlhhXhsJcQkjRwZp2kFph1iyz1JUmSJEmDsWKRm+QJSf40yZeT3JLkV/v9pye5IckdST6Z5Nh+/3H9/f398dOm+xYkTeIQj5toGwqzTho2s65j1knDZtZ1VvNOvg+8tKr+HvBcYEeSs4BLgfdV1TbgAeD8vv35wANV9SzgfX07SRtQAYfrcRNtA2LWSQNl1j2CWScNlFm3aMV3Up3v9ncf328FvBS4ut9/JfDq/vbO/j798bOTODlc2pDCoQm3oTDrpCEz6+aZddKQmXXzVrXwVJItwI3As4DLgK8B366qg32TOeCU/vYpwF0AVXUwyYPAicD9S57zQuBCgCfwA5zzuNdO9k6kxj2F419wtI+Z/8VPnWln3Zbjj+f03/yNab8NadCOPXWrWTchs07a+My6yazqU6iqQ1X1XGArcCbwo+Oa9X/H/QRQj9pRdXlVba+q7Y/nuNX2V5KmZtpZt+XJT1q7zkrSY2TWSRq6o7qEUFV9O8n1wFnA05Ic0//qtxW4p282B5wKzCU5BngqcGDtuixpLQ1paspaMeuk4THrHs2sk4bHrOusZnXlpyd5Wn/7icDLgNuAzwOv6ZvtAj7T397d36c//rmqetQvfpJmryouUNAz66ThMusWmXXScJl1i1YzknsycGV//sbjgE9V1R8kuRW4Ksl7gC8BH+7bfxj4WJL9dL/0nTeFfktaI4emHGhJdgC/CWwBPlRV711y/IfpFjV5Wt/m4qraM9VOjWfWSQM27azbRMw6acDMus6KRW5V3Qw8b8z+O+nO41i6/yHAVaQkzS9uchlwDt2Ut71JdlfVrSPN/m+6L1kfSHIGsAc4bb37atZJaoFZJ6kFR3VOrqRhKeDwdM/dOBPY3395IslVdJejGC1yC/jB/vZTWTwPTJLWxDpknSTNnFm3yCJXalqmPa1l4dITvTnghUva/Arwn5L878CT6M4Pk6Q1NPWsk6QNwKybZ5ErNay7ntrEv/idlGTfyP3Lq+ry/vZqLj3xeuAjVfUbSf5HunO/nlNVhyftmCTBmmWdJG1oZt0ii1ypcYdWd7nsI7m/qrYvc2z+0hPzRi9LMe98YAdAVX0hyROAk4D7Ju2YJM1bg6yTpA3PrOv4KUiapr3AtiSnJzmWblXO3UvafB04GyDJjwJPAL65rr2UJEnSYDiSKzWsyFSntVTVwSQXAdfQXR7oiqq6JcklwL6q2g38IvDBJP8H3UybN3sNRklradpZJ0kbgVm3yCJXatzhKU/o6K95u2fJvneO3L4VeNFUOyGpedPOOknaCMy6jkWu1LAqOOQvfpIGzqyT1AKzbpFFrtQ4p7VIaoFZJ6kFZl3H8WxJkiRJ0mA4kis1rFugwN+6JA2bWSepBWbdIotcqXGHcFqLpOEz6yS1wKzrWORKDSs8d0PS8Jl1klpg1i1yPFuSJEmSNBiO5EpN89wNSS0w6yS1wKyb56cgNe4wmWiTpM1g2lmXZEeS25PsT3LxmOPHJflkf/yGJKf1+09L8t+S3NRv/2bkMS9I8pX+Me9PYuhKOiKzruNIrtQwLxouqQXTzrokW4DLgHOAOWBvkt1VdetIs/OBB6rqWUnOAy4FXtcf+1pVPXfMU38AuBD4IrAH2AH84ZTehqRNzqxb5Eiu1LjD9biJNknaDKacdWcC+6vqzqp6GLgK2LmkzU7gyv721cDZRxqtSHIy8INV9YWqKuCjwKsfy3uX1A6zruM3VEmSpJWdlGTfyHbhyLFTgLtG7s/1+xjXpqoOAg8CJ/bHTk/ypST/X5KfGGk/t8JzStJaG0TWOV1Zalh30XCnK0satjXKuvuravsyx8Y9ea2yzb3AD1fVt5K8APj3SX5slc8pSQvMukUWuVLjXDxKUgumnHVzwKkj97cC9yzTZi7JMcBTgQP99LzvA1TVjUm+Bjy7b791heeUpEcw6zpOV5YaNn/R8Ek2Sdro1iHr9gLbkpye5FjgPGD3kja7gV397dcAn6uqSvL0fjEXkvwIsA24s6ruBb6T5Kz+fLY3AZ9Zkw9E0iCZdYscyZUkSZpAVR1MchFwDbAFuKKqbklyCbCvqnYDHwY+lmQ/cIDuyyHAi4FLkhwEDgFvqaoD/bF/CnwEeCLdSqOurCxpZjZT1lnkSo1zhWRJLZh21lXVHrpLX4zue+fI7YeA14553KeBTy/znPuA56xtTyUNmVnXsciVWuaUY0ktMOsktcCsW2CRKzWscOEpScNn1klqgVm3yCJXapy/+ElqgVknqQVmXceT8SRJkiRJg+FIrtSw+aXmJWnIzDpJLTDrFlnkSo0zDCW1wKyT1AKzrmORKzWscBU+ScNn1klqgVm3yCJXapyr8ElqgVknqQVmXceFpyRJkiRJg+FIrtSy8twNSQ0w6yS1wKxbYJErNcxV+CS1wKyT1AKzbpFFrtQ4w1BSC8w6SS0w6zqekytJkiRJGgxHcqWGudS8pBaYdZJaYNYtssiVGleGoaQGmHWSWmDWdSxypcZ5PTVJLTDrJLXArOtY5EoNK5eal9QAs05SC8y6RS48JUmSJEkaDEdypcZ57oakFph1klpg1nUscqWmuQqfpBaYdZJaYNbNs8iVGucvfpJaYNZJaoFZ17HIlRpWuECBpOEz6yS1wKxb5MJTkiRJkqTBcCRXall1y81L0qCZdZJaYNYtsMiVGudFwyW1wKyT1AKzrmORKzWscIECScNn1klqgVm3yHNyJUmSJEmD4Uiu1DSvpyapBWadpBaYdfMscqXGuUCBpBaYdZJaYNZ1LHKlxnnuhqQWmHWSWmDWdSxypYZVGYaShs+sk9QCs26RC09JkiRJkgbDIldq3OHKRJskbQbTzrokO5LcnmR/kovHHD8uySf74zckOa3ff06SG5N8pf/70pHHXN8/50399ow1/EgkDZBZ13G6stQ4FyiQ1IJpZl2SLcBlwDnAHLA3ye6qunWk2fnAA1X1rCTnAZcCrwPuB36qqu5J8hzgGuCUkce9oar2Ta/3kobErOs4kis1rioTbZK0GUw5684E9lfVnVX1MHAVsHNJm53Alf3tq4Gzk6SqvlRV9/T7bwGekOS4NXrbkhpj1nUscqWGFZMFoUWupM1gjbLupCT7RrYLR17iFOCukftzPHKE4hFtquog8CBw4pI2Pw18qaq+P7Lvt/rpe7+cxNCVtCyzbpHTlSVJklZ2f1VtX+bYuC9kSycNHrFNkh+jm9Z37sjxN1TV3UmeAnwaeCPw0dV3WZKO2iCybsWR3CSnJvl8ktuS3JLk5/v9JyS5Nskd/d/j+/1J8v7+ZOObkzx/kg5Kmq6acBsKs04atiln3Rxw6sj9rcA9y7VJcgzwVOBAf38r8HvAm6rqawt9rrq7//sd4LfppgpOxKyThs2s66xmuvJB4Ber6keBs4C3JjkDuBi4rqq2Adf19wFeAWzrtwuBD0zaSUlTUp6TO8Ksk4Zq+lm3F9iW5PQkxwLnAbuXtNkN7Opvvwb4XFVVkqcB/wF4e1X9yXzjJMckOam//XjgJ4GvTvxZmHXScJl1C1Yscqvq3qr6s/72d4Db6OZaj55UfCXw6v72TuCj1fki8LQkJ0/aUUlT4lAuYNZJgzfFrOvPO7uIbrXQ24BPVdUtSS5J8qq+2YeBE5PsB97GYhF5EfAs4JeXXD7jOOCaJDcDNwF3Ax+c8FMw66ShM+uAozwnN911jp4H3AA8s6ruhS4wR65ntNwJyfcuea4L6X4R5An8wGPouqTNIMkO4DeBLcCHquq9Y9r8Y+BX6OL1y1X1M+vayUf35zSmkHVbjj9+qv2WNDtVtQfYs2TfO0duPwS8dszj3gO8Z5mnfcFa9nEps07S0dosWbfqIjfJk+lOBP6FqvqrIyx6tZoTkqmqy4HLAX4wJwxoPEjaXKY55TiruJ5akm3A24EXVdUDWYMLgE9imll33A+fatZJMzKw0ysmZtZJw2TWdVZ1CaF+fvSngY9X1e/2u78xP12l/3tfv381JyRL2iCqJttWsJrrqf2vwGVV9UDXn7qPGTHrpOGactZtKmadNFxmXWc1qyuHbm71bVX1L0cOjZ5UvAv4zMj+N/Wr8Z0FPDg//UXSxlKsyQIFk15P7dnAs5P8SZIv9tOb151ZJw3XGmXdIJh10nCZdYtWM135RXTXKvpKkpv6ff8X8F7gU0nOB77O4tzrPcArgf3AXwP/ZE17LGntFDB5oE16PbVj6FbtfAndCMEfJ3lOVX170o4dJbNOGqq1ybqhMOukoTLrFqxY5FbVf2b8F1WAs8e0L+CtE/ZL0jCs9npqX6yqvwH+IsntdEXv3vXpYsesk9QCs05SC1Z1Tq6k4ZryuRuruZ7avwf+AUB/nbRnA3eu7buU1DrPU5PUArOuc1SXEJI0QFMMtKo6mGT+empbgCvmr6cG7Kuq3f2xc5PcChwC/s+q+tb0eiWpSQP68iZJyzLrAItcqXHTX2RgFddTK7qLhb9tqh2R1LBhLagiSeOZdfMscqXW+YufpBaYdZJaYNYBnpMrSZIkSRoQR3KllhVOa5E0fGadpBaYdQsscqXWOa1FUgvMOkktMOsAi1xJy14uUZKGxKyT1AKzDjwnV5IkSZI0II7kSq1zWoukFph1klpg1gEWuZIMQ0ktMOsktcCsAyxypbYV4Cp8kobOrJPUArNugUWu1LjyFz9JDTDrJLXArOu48JQkSZIkaTAcyZVa5y9+klpg1klqgVkHWORK8twNSS0w6yS1wKwDLHKl5sVf/CQ1wKyT1AKzrmORK7WscFqLpOEz6yS1wKxb4MJTkiRJkqTBcCRXalo8d0NSA8w6SS0w6+ZZ5Eqtc1qLpBaYdZJaYNYBFrmSDENJLTDrJLXArAM8J1eSJEmSNCCO5Eqt8xc/SS0w6yS1wKwDLHKlthUuUCBp+Mw6SS0w6xZY5EqN86Lhklpg1klqgVnX8ZxcqXU14SZJm8GUsy7JjiS3J9mf5OIxx49L8sn++A1JThs59vZ+/+1JXr7a55SkRzHrAEdytUrX3HPTwu2X/63nzrAnkiRtLEm2AJcB5wBzwN4ku6vq1pFm5wMPVNWzkpwHXAq8LskZwHnAjwF/C/hskmf3j1npOSVp3WymrHMkV5IkaTJnAvur6s6qehi4Cti5pM1O4Mr+9tXA2UnS77+qqr5fVX8B7O+fbzXPKUnradNknUWu1LjUZJskbQZTzrpTgLtG7s/1+8a2qaqDwIPAiUd47GqeU5IewazrOF1Zq+IU5QFzFT5JLZg8605Ksm/k/uVVdXl/e9yTL/26uFyb5faPG4jwp0VJR2bWARa5UttcPEpSC9Ym6+6vqu3LHJsDTh25vxW4Z5k2c0mOAZ4KHFjhsSs9pyQtMusWOF1ZkiRpMnuBbUlOT3Is3eIqu5e02Q3s6m+/BvhcVVW//7x+RdLTgW3An67yOSVpPW2arHMkV2qdI7mSWjDFrKuqg0kuAq4BtgBXVNUtSS4B9lXVbuDDwMeS7Kcb1Tivf+wtST4F3AocBN5aVYcAxj3n9N6FpEEw6wCLXKl5Lh4lqQXTzrqq2gPsWbLvnSO3HwJeu8xjfw34tdU8pyQdiVnXsciVWmeRK6kFZp2kFph1gOfkSpIkSZIGxJFcqXX+4iepBWadpBaYdYBFrtS0VV74W5I2NbNOUgvMukUWuVLrJr9ouCRtfGadpBaYdYBFriR/8ZPUArNOUgvMOsCFpyRJkiRJA+JIrtQ4z92Q1AKzTlILzLqORa7UOsNQUgvMOkktMOsAi1ypba7CJ6kFZp2kFph1CzwnV5IkSZI0GI7kSq3zFz9JLTDrJLXArAMsciUZhpJaYNZJaoFZB1jkSs3z3A1JLTDrJLXArOt4Tq4kSZIkaTAsciVJkiRJg+F0Zal1TmuR1AKzTlILzDrAIldqm9dTk9QCs05SC8y6BRa5UusMQ0ktMOsktcCsAyxyJRmGklpg1klqgVkHuPCUJEmSJGlAHMmVGhY8d0PS8Jl1klpg1i2yyJVaZxhKaoFZJ6kFZh1gkSu1zVX4JLXArJPUArNugefkSpqqJDuS3J5kf5KLj9DuNUkqyfb17J8kSZKGxZFcqXVT/MUvyRbgMuAcYA7Ym2R3Vd26pN1TgJ8DbphebyQ1zdENSS0w64BVjOQmuSLJfUm+OrLvhCTXJrmj/3t8vz9J3t+P2Nyc5PnT7LykNVATbkd2JrC/qu6sqoeBq4CdY9q9G/h14KHJ3sxjZ9ZJAzfdrNtUzDtpwMw6YHXTlT8C7Fiy72LguqraBlzX3wd4BbCt3y4EPrA23ZQ0LanJNuCkJPtGtgtHnv4U4K6R+3P9vsXXT54HnFpVfzDt97qCj2DWSYO1Blk3JB/BvJMGyazrrFjkVtUfAQeW7N4JXNnfvhJ49cj+j1bni8DTkpy8Vp2VNAWT/+J3f1VtH9kuH3n2LPOK3cHkccD7gF9c2zd19Mw6aeAc3Vhg3kkDZtYBj33hqWdW1b0A/d9n9PtXHLWZl+TC+ZGfv+H7j7Ebkja4OeDUkftbgXtG7j8FeA5wfZK/BM4Cdm+gxafWNOsOffd7U+2sJE1gorwz6yRtJGu9uvIRR20esbPq8vmRn8dz3Bp3Q9KqTPpr38q/+O0FtiU5PcmxwHnA7oWXr3qwqk6qqtOq6jTgi8CrqmrfGr3DaXlMWbflyU+acrckjTX9rBuyVeWdWSdtAGbdgsda5H5jfqpK//e+fv9KozaSNphpnrtRVQeBi4BrgNuAT1XVLUkuSfKq6b+7iZl10kB4ntqKzDtpAMy6zmMtcncDu/rbu4DPjOx/U78S31nAg/NTXyRtUFP+xa+q9lTVs6vq71TVr/X73llVu8e0fckGG8U166ShmOHoxnIrF49pt6tvc0eSXf2+H0jyH5L8lyS3JHnvSPs3J/lmkpv67YIJumneSUNg1gGru4TQJ4AvAP9dkrkk5wPvBc5Jcgfd9S/nO7EHuBPYD3wQ+N9Wen5Js+Uvfh2zThq2GWfdcisXL/YvOQF4F/BCusuvvWvkC+K/qKq/CzwPeFGSV4w89JNV9dx++9BqOmPeScNl1nWOWalBVb1+mUNnj2lbwFtXek5J2mjMOklTtBN4SX/7SuB64JeWtHk5cG1VHQBIci2wo6o+AXweoKoeTvJndFOGHzPzTtKUbJisW+uFpyRtNjOc1iJJ62byrDvSNcFXstzKxaNWc13xpwE/RTdCMu+nk9yc5Ooko+fOSmqRWQesYiRX0oBZqEpqwdpk3f1VtezlzZJ8FvihMYfescrnX+m64scAnwDeX1V39rt/H/hEVX0/yVvoRk5eusrXkzQ0Zt0Ci1ypYWF80kjSkKxH1lXVy5Z9/eQbSU6uqnuXrFw8ao7FaX7QTdO7fuT+5cAdVfWvRl7zWyPHPwhc+hi6LmkgzLpFTleWJEmaruVWLh51DXBukuP7RVjO7feR5D3AU4FfGH3A/CV/eq+iu1SbJM3Khsk6R3Kl1jldWVILZpt17wU+1a9i/HXgtQBJtgNvqaoLqq1a5RcAABLzSURBVOpAkncDe/vHXNLv20o3DfC/AH+WBOD/7VcX/bn+muMHgQPAm9fzTUnagMw6wCJXat6QLgMkScuZZdb1U+3GrVy8D7hg5P4VwBVL2syxzAzEqno78PY17aykTc2s61jkSq2zyJXUArNOUgvMOsAiV5JhKKkFZp2kFph1gAtPSZIkSZIGxJFcqWXlObmSGmDWSWqBWbfAIldqnWEoqQVmnaQWmHWARa7UPH/xk9QCs05SC8y6jkWu1DrDUFILzDpJLTDrABeekiRJkiQNiCO5UuOc1iKpBWadpBaYdR2LXKllhdNaJA2fWSepBWbdAotcqXWGoaQWmHWSWmDWAZ6TK0mSJEkaEEdypYYFz92QNHxmnaQWmHWLLHKl1hmGklpg1klqgVkHWORKzUuZhpKGz6yT1AKzrmORK7XMVfgktcCsk9QCs26BC09JkiRJkgbDkVypcS5QIKkFZp2kFph1HYtcqXWGoaQWmHWSWmDWARa5UvP8xU9SC8w6SS0w6zoWuVLrDENJLTDrJLXArANceEqSJEmSNCCO5EotK6e1SGqAWSepBWbdAotcqXWGoaQWmHWSWmDWARa5UtOCv/hJGj6zTlILzLpFnpMrSZIkSRoMR3Kl1pU/+UlqgFknqQVmHWCRKzXPaS2SWmDWSWqBWdexyJVaVrhAgaThM+sktcCsW2CRKzUuh2fdA0maPrNOUgvMuo4LT0mSJEmSBsORXKl1TmuR1AKzTlILzDrAkVypeanJNknaDGaZdUlOSHJtkjv6v8cv025X3+aOJLtG9l+f5PYkN/XbM/r9xyX5ZJL9SW5IctpkPZW02Zl1HYtcqWVFt9T8JJskbXSzz7qLgeuqahtwXX//EZKcALwLeCFwJvCuJV8Q31BVz+23+/p95wMPVNWzgPcBl07aUUmbmFm3wCJXapwjuZJaMOOs2wlc2d++Enj1mDYvB66tqgNV9QBwLbDjKJ73auDsJJm4t5I2LbOuY5ErSZK0spOS7BvZLjyKxz6zqu4F6P8+Y0ybU4C7Ru7P9fvm/VY/fe+XR77cLTymqg4CDwInHkW/JGmpQWSdC09JrXM0VlILJs+6+6tq+3IHk3wW+KExh96xyucfNyox3+s3VNXdSZ4CfBp4I/DRFR4jqUVmHWCRKzUtOOVY0vCtR9ZV1cuWff3kG0lOrqp7k5wM3Dem2RzwkpH7W4Hr++e+u//7nSS/TXce20f7x5wKzCU5BngqcGDydyNpMzLrFjldWWrZpIsTuPCUpM1g9lm3G9jV394FfGZMm2uAc5Mc3y/Cci5wTZJjkpwEkOTxwE8CXx3zvK8BPldlMEvNMusWOJIrSZI0Xe8FPpXkfODrwGsBkmwH3lJVF1TVgSTvBvb2j7mk3/ckui+Ajwe2AJ8FPti3+TDwsST76UY1zlu/tyRJj7Jhss4iV2qc05UltWCWWVdV3wLOHrN/H3DByP0rgCuWtPke8IJlnvch+i+RkgRm3TyLXKl1FrmSWmDWSWqBWQdY5ErNcyRXUgvMOkktMOs6FrlSywo4bBpKGjizTlILzLoFrq4sSZIkSRoMR3Kl1vmDn6QWmHWSWmDWARa5UvM8d0NSC8w6SS0w6zoWuVLrJr/wtyRtfGadpBaYdYDn5ErNS022rfj8yY4ktyfZn+TiMcffluTWJDcnuS7J357G+5TUtmlnnSRtBGZdxyJX0tQk2QJcBrwCOAN4fZIzljT7ErC9qv574Grg19e3l5IkSRoSi1ypZbUG25GdCeyvqjur6mHgKmDnI7pQ9fmq+uv+7heBrZO+LUl6hOlnnSTNnlm3wHNypYYFyHTP3TgFuGvk/hzwwiO0Px/4w2l2SFJ71iHrJGnmzLpFFrlS6w5P/AwnJdk3cv/yqrq8v50x7cemb5KfBbYDf3/iHknSUpNnnSRtfGYdYJEraXL3V9X2ZY7NAaeO3N8K3LO0UZKXAe8A/n5VfX/tuyhJkqRWWORKjZvytJa9wLYkpwN3A+cBP/OI10+eB/xbYEdV3TfNzkhql1P4JLXArOtMZeGplS4ZImmDmPICBVV1ELgIuAa4DfhUVd2S5JIkr+qb/XPgycDvJLkpye41fIdTZ95Jm4CLsUzMrJM2AbNuwZqP5I5cMuQcuqmKe5Psrqpb1/q1JE2qpn7R8KraA+xZsu+dI7dfNtUOTJF5J20W08+6ITPrpM3CrJs3jZHcFS8ZImnj8KLhEzHvpE3CrJuIWSdtEmZdZxpF7rhLhpyytFGSC5PsS7Lvb3CdGUmb0op5N5p1h777vXXtnCStEbNO0qYyjYWnVnXJkP4SI5cDbN++va7d9ztT6IrUjiQ3PqYHOq1lEivm3dKs2/fzv7ge/ZIGK7/wz8y69WfWSevMrJvMNIrcVV0yRNIGUBCvpzYJ807aDMy6SZl10mZg1i2YxnTlhUuGJDmW7pIhm2q1VKkpVZNtbTPvpM3CrJuEWSdtFmYdMIWR3Ko6mGT+kiFbgCuq6pa1fh1JmjXzTlILzDpJm800piuPvWSIpA1qOD/azYR5J20SZt1EzDppkzDrgCkVuZI2jwxoaookLcesk9QCs65jkSu1zjCU1AKzTlILzDrAIldqWwGuwidp6Mw6SS0w6xZMY3VlSZIkSZJmwpFcqWGhPHdD0uCZdZJaYNYtssiVWmcYSmqBWSepBWYdYJEryTCU1AKzTlILzDrAIldqmwsUSGqBWSepBWbdAheekiRJkiQNhkWu1LhUTbRJ0mYwy6xLckKSa5Pc0f89fpl2u/o2dyTZ1e97SpKbRrb7k/yr/tibk3xz5NgFE3VU0qZn1nWcriy1zkJVUgtmm3UXA9dV1XuTXNzf/6XRBklOAN4FbKebdHhjkt1V9QDw3JF2NwK/O/LQT1bVRdN+A5I2CbMOcCRXalx1YTjJJkkb3syzbidwZX/7SuDVY9q8HLi2qg70X/auBXaMNkiyDXgG8MeTdkjSEJl18yxyJUmSVnZSkn0j24VH8dhnVtW9AP3fZ4xpcwpw18j9uX7fqNfTjWaMfhP96SQ3J7k6yalH0SdJGmcQWed0ZallhaOxkoZvbbLu/qravtzBJJ8FfmjMoXes8vkzZt/STp8HvHHk/u8Dn6iq7yd5C93IyUtX+XqShsasW2CRK7XOpeYltWDKWVdVL1vuWJJvJDm5qu5NcjJw35hmc8BLRu5vBa4feY6/BxxTVTeOvOa3Rtp/ELj0sfVe0mCYdYDTlaXmubqypBbMOOt2A7v627uAz4xpcw1wbpLj+xVJz+33zXs98IlHvKfuS+S8VwG3TdpRSZubWddxJFdqnYWqpBbMNuveC3wqyfnA14HXAiTZDrylqi6oqgNJ3g3s7R9zSVUdGHmOfwy8csnz/lySVwEHgQPAm6f4HiRtBmYdYJErSZI0Vf1Uu7PH7N8HXDBy/wrgimWe40fG7Hs78Pa166kkPXYbKesscqWWFXDYkVxJA2fWSWqBWbfAIldqmte6ldQCs05SC8y6eRa5UusMQ0ktMOsktcCsAyxyJRmGklpg1klqgVkHeAkhSZIkSdKAOJIrtcwFCiS1wKyT1AKzboFFrtS0gjo8605I0pSZdZJaYNbNs8iVWue5G5JaYNZJaoFZB3hOriRJkiRpQBzJlVrmuRuSWmDWSWqBWbfAIldqndNaJLXArJPUArMOsMiVZBhKaoFZJ6kFZh1gkSs1rgxDSQ0w6yS1wKyb58JTkiRJkqTBcCRXalkBh72emqSBM+sktcCsW2CRK7XOaS2SWmDWSWqBWQdY5EoyDCW1wKyT1AKzDrDIlRpXXk9NUgPMOkktMOvmufCUJEmSJGkwHMmVWlZQ5QIFkgbOrJPUArNugUWu1DqntUhqgVknqQVmHWCRK8kFCiS1wKyT1AKzDvCcXEmSJEnSgDiSK7WsyouGSxo+s05SC8y6BRa5Uuuc1iKpBWadpBaYdYBFrtS88hc/SQ0w6yS1wKzrWORKTSt/8ZPUALNOUgvMunkuPCVJkiRJGgxHcqWWFV5PTdLwmXWSWmDWLbDIlVpXnrshqQFmnaQWmHWARa7UtALKX/wkDZxZJ6kFZt0iz8mVWlbV/eI3ybaCJDuS3J5kf5KLxxw/Lskn++M3JDltCu9UUsvWIeuOJMkJSa5Nckf/9/hl2v3HJN9O8gdL9p/e5+MdfV4e2+83PyUtMusWWORKmpokW4DLgFcAZwCvT3LGkmbnAw9U1bOA9wGXrm8vJWnqLgauq6ptwHX9/XH+OfDGMfsvBd7XP/4ButwE81PSxrJhss4iV2pcHa6JthWcCeyvqjur6mHgKmDnkjY7gSv721cDZyfJmr5JSc2bctatZDTnrgRePbaPVdcB3xnd1+fhS+nycenjzU9Jj2DWdTwnV2rddBcoOAW4a+T+HPDC5dpU1cEkDwInAvdPs2OSGjPbxVieWVX3AlTVvUmecRSPPRH4dlUd7O/P0eUmmJ+SljLrgA1S5N54443fTXL7rPvRO4mN838O9mU8+zLe3z7aB3yHB675bF190oSv+4Qk+0buX15Vl/e3x/3KtvRnwtW0GYQbb7zx/iTfY+P8m9lI/37ty3j25dE2YtaR5LPAD4153DsmfN0jZeSGzE+z7ojsy3j25dHMugmybkMUucDtVbV91p0ASLLPvjyafRlvI/XlsaiqHVN+iTng1JH7W4F7lmkzl+QY4KnAgSn3ayaq6ukb6d+MfRnPvoy3kfpytNYh66iqly13LMk3kpzcj2ycDNx3FE99P/C0JMf0IxyjOboh89OsW559Gc++rA2zbpHn5Eqapr3Atn61vGOB84DdS9rsBnb1t18DfK6qZj4SIUlraDTndgGfWe0D+zz8PF0+Ln28+SlpI9kwWWeRK2lq+l/iLgKuAW4DPlVVtyS5JMmr+mYfBk5Msh94G8uvxCdJm9V7gXOS3AGc098nyfYkH5pvlOSPgd+hW1RlLsnL+0O/BLytz8kT6XITzE9JG8uGybqNMl358pWbrBv7Mp59GW8j9WVDqqo9wJ4l+945cvsh4LXr3a8Z2kj/ZuzLePZlvI3Ul02lqr4FnD1m/z7ggpH7P7HM4++kW61+6f6NnJ8b6d+LfRnPvoy3kfqyqWykrIuzWiRJkiRJQ+F0ZUmSJEnSYMy8yE2yI8ntSfYnWfdzSZL8ZZKvJLlpfrnsJCckuTbJHf3f46f02lckuS/JV0f2jX3tdN7ff043J3n+OvTlV5Lc3X82NyV55cixt/d9uX1kHv1a9eXUJJ9PcluSW5L8fL9/3T+bI/RlJp+NNi+zzqwb0xezToNj1pl1Y/pi1mn9VdXMNmAL8DXgR4BjgS8DZ6xzH/4SOGnJvl8HLu5vXwxcOqXXfjHwfOCrK7028ErgD+muE3UWcMM69OVXgH82pu0Z/X+r44DT+/+GW9awLycDz+9vPwX48/411/2zOUJfZvLZuG3Ozawz65bpi1nnNqjNrDPrlumLWee27tusR3LPBPZX1Z1V9TBwFbBzxn2Crg9X9revBF49jRepqj/i0dd4Wu61dwIfrc4X6a4jdfKU+7KcncBVVfX9qvoLYD9jThKfoC/3VtWf9be/Q7cq7ynM4LM5Ql+WM9XPRpuWWWfWjeuLWaehMevMunF9Meu07mZd5J4C3DVyf44j/0ObhgL+U5Ibk1zY73tmVd0L3f8YgGesY3+We+1ZfVYX9VNFrhiZ3rNufUlyGvA84AZm/Nks6QvM+LPRprIR/l2YdUdm1o3vC5h1Wr2N8O/CrDsys258X8CsG5RZF7kZs2+9l3t+UVU9H3gF8NYkL17n11+tWXxWHwD+DvBc4F7gN9azL0meDHwa+IWq+qsjNZ12f8b0ZaafjTadjfDvwqxbnlm3fF/MOh2NjfDvwqxbnlm3fF/MuoGZdZE7B5w6cn8rcM96dqCq7un/3gf8Ht0UhG/MT4vo/963jl1a7rXX/bOqqm9U1aGqOgx8kMXpGVPvS5LH04XPx6vqd/vdM/lsxvVllp+NNqWZ/7sw65Zn1i3fF7NOR2nm/y7MuuWZdcv3xawbnlkXuXuBbUlOT3IscB6we71ePMmTkjxl/jZwLvDVvg+7+ma7gM+sV5+O8Nq7gTf1K86dBTw4P8VjWpac//AP6T6b+b6cl+S4JKcD24A/XcPXDfBh4Laq+pcjh9b9s1muL7P6bLRpmXWPZtaZdRoes+7RzDqzTrNQM175im4FtT+nW63sHev82j9Ct2Lal4Fb5l8fOBG4Drij/3vClF7/E3RTIv6G7pei85d7bbrpEpf1n9NXgO3r0JeP9a91M93/yE8eaf+Ovi+3A69Y4778z3RTQW4Gbuq3V87iszlCX2by2bht3s2sM+vG9MWscxvcZtaZdWP6Yta5rfuW/j+eJEmSJEmb3qynK0uSJEmStGYsciVJkiRJg2GRK0mSJEkaDItcSZIkSdJgWORKkiRJkgbDIleSJEmSNBgWuZIkSZKkwbDIlSRJkiQNxv8PtjZY9H+HVNYAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1152x432 with 6 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"init()\n",
"plot()\n",
"print(dh)"
]
},
{
"cell_type": "code",
"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",
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
" 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",
}
},
"nbformat": 4,
"nbformat_minor": 2
}