{ "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": [ "# Tutorial 05: Phase-field simulation of dentritic solidification\n", "\n", "This is the second tutorial on phase field methods with pystencils. Make sure to read the previous tutorial first. \n", "\n", "In this tutorial we again implement a model described in **Programming Phase-Field Modelling** by S. Bulent Biner.\n", "This time we implement the model from chapter 4.7 that describes dentritic growth. So get ready for some beautiful snowflake pictures.\n", "\n", "We start again by adding all required arrays fields. This time we explicitly store the change of the phase variable φ in time, since the dynamics is calculated using an Allen-Cahn formulation where a term $\\partial_t \\phi$ occurs." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "dh = ps.create_data_handling(domain_size=(300, 300), periodicity=True, \n", " default_target='cpu')\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": "markdown", "metadata": {}, "source": [ "This model has a lot of parameters that are created here in a symbolic fashion. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$\\frac{{{φ}_{C}}^{4}}{4} - {{φ}_{C}}^{3} \\left(- \\frac{m}{3} + \\frac{1}{2}\\right) + {{φ}_{C}}^{2} \\left(- \\frac{m}{2} + \\frac{1}{4}\\right) + \\frac{ε^{2}}{2} \\left({\\partial_{0} {{φ}_{C}}}^{2} + {\\partial_{1} {{φ}_{C}}}^{2}\\right)$$" ], "text/plain": [ " 4 2 ⎛ 2 2⎞\n", "φ_C 3 ⎛ m 1⎞ 2 ⎛ m 1⎞ ε ⋅⎝D(phi_C) + D(phi_C) ⎠\n", "──── - φ_C ⋅⎜- ─ + ─⎟ + φ_C ⋅⎜- ─ + ─⎟ + ──────────────────────────\n", " 4 ⎝ 3 2⎠ ⎝ 2 4⎠ 2 " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ε, m, δ, j, θzero, α, γ, Teq, κ, τ = sp.symbols(\"ε m δ j θ_0 α γ T_eq κ τ\")\n", "εb = sp.Symbol(\"\\\\bar{\\\\epsilon}\")\n", "\n", "φ = φ_field.center\n", "T = t_field.center\n", "\n", "def f(φ, m):\n", " return φ**4 / 4 - (frac(1, 2) - m/3) * φ**3 + (frac(1,4)-m/2)*φ**2\n", "\n", "free_energy_density = ε**2 / 2 * (ps.fd.Diff(φ,0)**2 + ps.fd.Diff(φ,1)**2 ) + f(φ, m)\n", "free_energy_density" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The free energy is again composed of a bulk and interface part." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEWCAYAAAAZwvJqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8lOW9///XJ3sC2UkCCWvYERAwsmgV64pVQa22aLVotaCn1nNOz2lr23Nqj63ftqc/62mtraK1Ym3r2ioqapEqagUk7PsWIiRhyUISCNlz/f7IYNMYJCGTuWcy7+fjMY+ZuZe5P7mzvHPf93VflznnEBERCRcRXhcgIiISSAo+EREJKwo+EREJKwo+EREJKwo+EREJKwo+EREJKwo+kdNkZs7MRvheP2lmP+rkellm9q6ZHTWzB3q2ShFpL8rrAkS8YmaFQBbQDDQCHwB3OOf29/Cm5wNlQJLTjbQiAacjPgl3Vznn+gIDgEPAQwHY5hBg68lCz8xC8h/SUK1bwo+CTwRwztUBLwDjTkwzs3fM7PY2728xs/dP9Vlmlmhmb5vZL83M2s17EpgHfMvMjpnZxWb2AzN7wcyeNrNq4BYzizCze8xsj5mVm9lzZpbW5nOmm9kHZlZpZhvM7IJPqSfbzF40s1Iz22tmd7eZ9wPfZz/lO/W6xczyurBu+7rjzWyRmR0xs21m9i0zK/It/00ze7FdbQ+Z2f+dap+K+JOCTwQwswTgi8DKbn5OOrAM+Ltz7u72R3XOuVuAPwD/65zr65x7yzdrDq3Bm+KbfzdwNTATyAaOAA/7tpEDvAb8CEgD/hN40cwyOqgnAngF2ADkABcB/2Zml7VZbDbwjG/bi4FfdWHd9nXfCwwFcoFLgJvaLPs0MMvMUnyfH0XrPv/9SXanSI9Q8Em4e8nMKoFqWv9Q/6wbn5UNLAeed879VxfXXeGce8k51+KcqwUWAN9zzhU55+qBHwDX+cLiJmCJc26Jb/mlQD7wuQ4+92wgwzl3n3OuwTlXADwGzG2zzPu+z2qmNYTO7MK67ev+AvD/nHNHnHNFwC9PLOicOwC8C1zvmzQLKHPOrenivhLpFp2Tl3B3tXPuLTOLpPXoZbmZjXPOHTyNz7oCOAY8chrrtm9QMwT4i5m1tJnWTGtjnCHA9WZ2VZt50cDbHXzuECDbF+4nRALvtXnf9ms9DsT5ArYz67avO7vdtPbzFwF30hqgN6GjPfGAjvhEAOdcs3Puz7SGy2d8k2uAhDaL9T/FxzwGvAEsMbM+XS2h3fv9wOXOuZQ2jzjnXLFv3u/bzevjnPtJB5+7H9jbbtlE51xHR4ens277ug8AA9u8H9Ru/kvARDMbD1xJ6+lRkYBS8IkA1moOkAps801eD1xrZgm++/Vu68RH3QXsAF41s/hulPQIcL+ZDfHVl+GrD1qvlV1lZpeZWaSZxZnZBWY2sIPP+RCoNrNv+xqeRJrZeDM7uxM1nM66zwHfMbNU37XIu9rObNOI6I/Ah865fZ2oQ8SvFHwS7l4xs2O0XuO7H5jnnNvim/cg0EDrbQ6L6MTRia8xy3xaj5ZeNrO406zrF7Q2NPmrmR2ltdHNNN829tN6Wva7QKlvW9+kg99n33W7q4BJwF5a7x98HEjuxNdyOuveBxT5ln+L1pCrb7fMImACOs0pHjHdPysiPcXM7gTmOudmtpk2GNgO9HfOVXtWnIQtHfGJiN+Y2QAzO9d3H+Jo4D+Av7SZHwF8A3hGoSdeUatOEfGnGOBRYBhQSev9gb8G8DX4OQR8ROutDCKe0KlOEREJKzrVKSIiYSUkT3X269fPDR061OsyREQkSKxZs6bMOfeJbvs6EpLBN3ToUPLz870uQ0REgoSZfdTZZXWqU0REwoqCT0REwoqCT0REwoqCT0REwoqCT0REwoqCT0REwoqCT0REwkpYBt/Bqjoe+OsOPiqv8boUEREJsLAMvoamFh76226W7yz1uhQREQmwsAy+QWnxZCfHsbKg3OtSREQkwMIy+MyM6bnprCyoQKNTiIiEl7AMPoDpw9OpqGlg56FjXpciIiIBFLbBNyM3HUCnO0VEwkzYBt/A1HhyUuIVfCIiYSZsg+8f1/nKaWnRdT4RkXARtsEHMD03jSPHG9l5+KjXpYiISID4JfjMbJaZ7TCz3WZ2Twfzv2FmW81so5ktM7MhbebNM7Ndvsc8f9TTWdNPXOfbo9OdIiLhotvBZ2aRwMPA5cA44AYzG9dusXVAnnNuIvAC8L++ddOAe4FpwFTgXjNL7W5NnTUoLYGBqfGs0HU+EZGw4Y8jvqnAbudcgXOuAXgGmNN2Aefc28654763K4GBvteXAUudcxXOuSPAUmCWH2rqtOm56azaW6HrfCIiYcIfwZcD7G/zvsg37WRuA17v6rpmNt/M8s0sv7TUf12NTc9Np/J4IzsO6TqfiEg48EfwWQfTOjx8MrObgDzgZ11d1zm30DmX55zLy8jIOK1COzI9Nw2AFbrOJyISFvwRfEXAoDbvBwIl7Rcys4uB7wGznXP1XVm3Jw1MTWBQmu7nExEJF/4IvtXASDMbZmYxwFxgcdsFzGwy8CitoXe4zaw3gUvNLNXXqOVS37SAmj5M1/lERMJFt4PPOdcE3EVrYG0DnnPObTGz+8xstm+xnwF9gefNbL2ZLfatWwH8kNbwXA3c55sWUDOGp1NV28j2g7rOJyLS20X540Occ0uAJe2mfb/N64s/Zd0ngCf8Ucfpmua7n29FQTnjspO8LEVERHpYWPfcckJOSjyD0xJ0nU9EJAwo+Hym56bxoa7ziYj0ego+nxPX+bYeqPa6FBER6UEKPp9pwzQ+n4hIOFDw+WSnxDMkPYGVBQFvVCoiIgGk4GtjRm46q/aW06zrfCIivZaCr43puekcrWtim67ziYj0Wgq+Nj4en0/X+UREei0FXxv9k+MY1q+POqwWEenFFHztnLifT9f5RER6JwVfO9Nz0zla38TWEl3nExHpjRR87Uz/uN/OMo8rERGRnqDgaycrKY4RmX1ZvtN/o7yLiEjwUPB14KKxmawqqKC6rtHrUkRExM8UfB24eGwWTS2Od3XUJyLS6yj4OjBlcCqpCdG8tfWQ16WIiIif+SX4zGyWme0ws91mdk8H8883s7Vm1mRm17Wb1+wblf3jkdm9FhlhfHZ0Jm/vKKWpucXrckRExI+6HXxmFgk8DFwOjANuMLNx7RbbB9wC/LGDj6h1zk3yPWZ3tx5/uXhcFlW1jaz56IjXpYiIiB/544hvKrDbOVfgnGsAngHmtF3AOVfonNsIhMzh03kj+xEdaSzbftjrUkRExI/8EXw5wP4274t80zorzszyzWylmV3th3r8IjEumum56by1Tdf5RER6E38En3UwrSv9fQ12zuUBNwL/Z2bDO9yI2XxfQOaXlgamteVFYzIpKK2hoPRYQLYnIiI9zx/BVwQMavN+IFDS2ZWdcyW+5wLgHWDySZZb6JzLc87lZWRknH61XXDR2CwAlm3T6U4Rkd7CH8G3GhhpZsPMLAaYC3SqdaaZpZpZrO91P+BcYKsfavKLQWkJjOmfqNOdIiK9SLeDzznXBNwFvAlsA55zzm0xs/vMbDaAmZ1tZkXA9cCjZrbFt/pYIN/MNgBvAz9xzgVN8EFrLy75Hx2h6rh6cRER8SfnHI+9W8DOQ0cDut0of3yIc24JsKTdtO+3eb2a1lOg7df7AJjgjxp6ykVjs3j47T28s/MwcyZ1pc2OiIh8mqIjtdy/ZBsJsZGMykoM2HbVc8spTBqYQr++Mbyl63wiIn61qbgKgAk5yQHdroLvFCIijAvHZPLOjsM0qhcXERG/2VxcRVSEMbp/4I72QMHXKReNzeJoXROr91Z4XYqISK+xqbiKUVmJxEZFBnS7Cr5OOG9kP2KiInS6U0TET5xzbCmpDvhpTlDwdUpCTBTnDE9n2fZDONeVe/NFRKQjJVV1VNQ0MH6ggi9oXTw2i4/Kj7NHvbiIiHTbpqLWhi3js5MCvm0FXyddNDYTgKVbdbpTRKS7tpRUERlhjB2g4AtaA5LjOSM7iWXqxUVEpNs2FVcxMrMvcdGBbdgCCr4uuWhsFmv3HaGipsHrUkREQpZzjs3FVYz3oGELKPi65JKxWbQ4eFtj9ImInLZD1fWUHWvwpEUnKPi6ZHxOEllJseq0WkSkG0702DI+J/DX90DB1yVmxqXj+vP2jsMcrVOn1SIip2NTcRURBuMG6IgvJFwzJYe6xhZe33zQ61JERELSluIqRmT2JT4m8A1bQMHXZZMHpTCsXx/+vLbI61JERELSpuIqxmd7c7QHCr4uMzOunZzDyoIKio4c97ocEZGQcri6jsNH6z1r0QkKvtNy9eTWcfleWlfscSUiIqFlc4lvKCIPuio7QcF3GgalJTBtWBovri1W350iIl2wqagaMxjnQY8tJ/gl+MxslpntMLPdZnZPB/PPN7O1ZtZkZte1mzfPzHb5HvP8UU8gfH7KQPaW1bBuf6XXpYiIhIxNxVXk9utDn9goz2rodvCZWSTwMHA5MA64wczGtVtsH3AL8Md266YB9wLTgKnAvWaW2t2aAuHyCf2Ji45QIxcRkS7YUlLl2Y3rJ/jjiG8qsNs5V+CcawCeAea0XcA5V+ic2wi0H8L8MmCpc67COXcEWArM8kNNPS4xLprLzujPKxsOUN/U7HU5IiJBr+xYPQeq6jxt2AL+Cb4cYH+b90W+aX5d18zmm1m+meWXlpaeVqH+du2UgVTVNqoLMxGRTvhHjy2hH3zWwbTOtvjo9LrOuYXOuTznXF5GRkani+tJ5w5PJzMxlhfXqnWniMipbPEF3xkejMHXlj+CrwgY1Ob9QKAkAOt6Lioygqsn5/D29sOUH6v3uhwRkaC2qbiKYf36kBgX7Wkd/gi+1cBIMxtmZjHAXGBxJ9d9E7jUzFJ9jVou9U0LGddOyaGpxfHKhpDJaxERT2wurvb8NCf4Ificc03AXbQG1jbgOefcFjO7z8xmA5jZ2WZWBFwPPGpmW3zrVgA/pDU8VwP3+aaFjDH9kxg3IIk/62Z2EZGTqqhpoLiylvEen+YE8MuNFM65JcCSdtO+3+b1alpPY3a07hPAE/6owyvXTsnhR69tY9eho4zMSvS6HBGRoLPZd33P61sZQD23+MWcSTlERpiO+kRETuJEi84zFHy9Q0ZiLDNHZfDSumKaW9SFmYhIe1tKqhiclkByvLcNW0DB5zfXTsnhQFUdK/aUe12KiEjQ2VTsfY8tJyj4/OTisVkkxkWpCzMRkXaqjjeyv6I2KFp0goLPb+KiI7ly4gBe33yQ6rpGr8sREQkaJ4YiGp/jfYtOUPD51Y1Th1Db2MwzH+7zuhQRkaDxcVdlHo663paCz48mDExmRm46T7xfSENT+/64RUTC0+biKgamxpPaJ8brUgAFn9/Nn5nLweo6FqsnFxERoDX4guVoDxR8fnfBqAxGZyXy2LsFGp1dRMJedV0jheXHmTBQwddrmRlfPT+XHYeO8s7O4Bg+SUTEK5uDZCiithR8PWD2mdn0T4pj4fICr0sREfHUun2VAExU8PVuMVER3HruUFYUlLOpqMrrckREPLO6sIKRmX2DpmELKPh6zA3TBtM3NopH393jdSkiIp5obnGsKTzC2cPSvC7lnyj4ekhSXDQ3ThvMkk0H2F9x3OtyREQCbvvBao7WNzF1qIIvbNx67lAizPjt+3u9LkVEJODyC48AkDc01eNK/pmCrwcNSI5n9qRsnl29nyM1DV6XIyISUB8WVpCdHMfA1ASvS/knfgk+M5tlZjvMbLeZ3dPB/Fgze9Y3f5WZDfVNH2pmtWa23vd4xB/1BJP55+dS29jM0ys/8roUEZGAcc6RX1hBXpCd5gQ/BJ+ZRQIPA5cD44AbzGxcu8VuA44450YADwI/bTNvj3Nuku9xR3frCTZj+icxc1QGi1YUUtfY7HU5IiIBsb+ilkPV9UHXsAX8c8Q3FdjtnCtwzjUAzwBz2i0zB1jke/0CcJGZmR+2HRIWnJ9L2bEG/rxWI7SLSHhYXVgBwNlBdn0P/BN8OcD+Nu+LfNM6XMY51wRUAem+ecPMbJ2ZLTez8/xQT9CZMTyd8TlJPP5egUZoF5GwsLqwguT4aEZlJnpdyif4I/g6OnJr/9f9ZMscAAY75yYD3wD+aGYdDthkZvPNLN/M8ktLQ6srMDNj/vnDKSir4bVNB7wuR0Skx60urCBvSCoREcF3cs8fwVcEDGrzfiDQfmiCj5cxsyggGahwztU758oBnHNrgD3AqI424pxb6JzLc87lZWRk+KHswLpiwgDGDUjiJ0u2Uduga30i0nuVH6tnT2lNUDZsAf8E32pgpJkNM7MYYC6wuN0yi4F5vtfXAX9zzjkzy/A1jsHMcoGRQK/s4DIywvjB7DMoqarjkeXqzUVEeq/Vvvv3pg4Lvut74Ifg812zuwt4E9gGPOec22Jm95nZbN9ivwXSzWw3rac0T9zycD6w0cw20Nro5Q7nXEV3awpWU4elceXEATyyfA9FR9Sbi4j0TvmFFcRERQTViAxtWSiOGZeXl+fy8/O9LuO0lFTWcuED73DhmEx+/aWzvC5HRMTv5vzqfWKjI3luwYyAbdPM1jjn8jqzrHpuCbDslHjunDmCJZsO8sGeMq/LERHxq+MNTWwuqQ7K2xhOUPB5YMHMXHJS4rnvla00Nbd4XY6IiN+s21dJc4vj7CBt2AIKPk/ERUfyX1eMZfvBo/zpw31elyMi4jerCyuIMDhriI74pJ1Z4/szIzedB5buVAfWItJrrC6sYEz/JBLjor0u5aQUfB4xM+6dPY7q2kZ+vnSn1+WIiHRbY3ML6/ZVMjUI++dsS8HnoTH9k7hp+hD+sOojth2o9rocEZFu2VpSzfGG5qAbf689BZ/HvnHJKJLio/mfV7YQireWiIic8I+OqXXEJ58iJSGG/7hkFCsLKnh1o/rxFJHQtbqwgsFpCWQlxXldyqdS8AWBG6YOZkJOMt/7yyb2V6hHFxEJPa0Dzx4J+qM9UPAFhajICB6+cQoAd/5hjQasFZGQU1BWQ3lNQ9D2z9mWgi9IDE5P4IEvTGJzcTX3vbrV63JERLpk9d7W63vBOiJDWwq+IHLJuCzumDmcP67ax5/XFnldjohIp60uPEJ6nxhy+/XxupRTUvAFmf+8dBTThqXxvb9sZsfBo16XIyLSKasLK8gbmopZ8A08256CL8hERUbw0I2T6RsXxZ1/WMOx+iavSxIR+VSHquvYV3E8JBq2gIIvKGUmxvHQDZMpLKvh2y9u1P19IhLUQuX+vRMUfEFqem4637xsDK9tPMCiDwq9LkdE5KRWFVSQEBPJuOwkr0vpFAVfEFtwfi4Xj83i/iXbWLvviNfliIh8gnOO5TtLOWd4OtGRoREpfqnSzGaZ2Q4z221m93QwP9bMnvXNX2VmQ9vM+45v+g4zu8wf9fQWERHGA9efSf/kOOY/lc/WEvXnKSLBpaCshn0Vx7lgdKbXpXRat4PPzCKBh4HLgXHADWY2rt1itwFHnHMjgAeBn/rWHQfMBc4AZgG/9n2e+CQnRPPkrVOJjoxg7sIVOvITkaDy9vbDAFwwOsPjSjovyg+fMRXY7ZwrADCzZ4A5QNu7sOcAP/C9fgH4lbW2eZ0DPOOcqwf2mtlu3+et8ENdvcbwjL48t2AGN/12FTc9vorH5+VxzvB+XpfVqzW3OEoqaymurKW6tpHquibfcyNHfa9rGpqIioggNiqCuOhIYqMiiI2OIDaq9XVmUiyD0/owOC2Bfn1jQqKZt0hXLd9ZysjMvgxMTfC6lE7zR/DlAPvbvC8Cpp1sGedck5lVAem+6SvbrZvT0UbMbD4wH2Dw4MF+KDu0DEpL4Hlf+N3yu9U8ctMULhyT5XVZIe9gVR07Dx2lsLyGwrLjrc/lNeyvOE5jc8etafvERJIUH02f2CiaWxx1jc3UN7VQ73tuavnkegkxkQxOS2BQWgKD0xIYnZXI1GFpDElPUCBKyKqpb2JVQQW3nDvU61K6xB/B19Fvbfvf/JMt05l1Wyc6txBYCJCXlxeW7fszk+J4Zv4M5j3xIfOfWsP/zZ3ElROzvS4rZFQdb2RjcSUb9leyoaiKDfsrOXy0/uP58dGRDO3Xh9FZiVx2Rn+GpicwMDWB5PhokuKiSYqPom9sFFGnuIDf1NxCXVMLB6vq2F9xnI/Ka9hXUcs+3+v3dpVS19gCQGZiLFOHpTFtWBpTh6UzMrMvEREKQgkNH+wpp6G5hQtGhc5pTvBP8BUBg9q8HwiUnGSZIjOLApKBik6uK22k9Ynhj1+dxm1P5nP3n9ZxvL6ZL5w96NQrhqH9FcdZsaecFQXlrNt3hMLyf4x8kZvRh3NH9GPiwGTGDkhiWL8+ZCbG+uXoKyoygr6REYzI7MuIzL6fmO+cY0/pMVbtreDDvRWsajMkVUpCNDNy07lyYjYXjc0kLlqXvCV4vb3jMH1iIkOif862/BF8q4GRZjYMKKa1scqN7ZZZDMyj9drddcDfnHPOzBYDfzSznwPZwEjgQz/U1KslxkWz6CtTWfD0Gr714kaq6xq57TPDwv6U2aHqOlbsKeeDPWV8sKecoiO1APTrG8NZQ1K5Pm8QkwalMD4nmeT4aM/qNDNGZCYyIjORL00bgnOOoiO1viAs550dpby++SB9Y6O47Iz+zJmUzTnD0095pCkSSM45lu8o5TMj+xETFVo/m90OPt81u7uAN4FI4Ann3BYzuw/Id84tBn4L/N7XeKWC1nDEt9xztDaEaQK+5pzTmDydEB8TyWNfPot//dN6fvTaNv6+u4wfXTOBnJR4r0sLmLrGZlYWlLN8ZynLd5ZSUFoDQHJ8NNNz0/jqebnMGN56+jCY/ykwMwb5rv9dd9ZAmlscqwrKeWl9Ma9vPsiLa4vo1zeWKycOYM6kbCYNSgnqr0fCw67DxyiurOXrF47wupQus1DsDisvL8/l5+d7XUZQaG5xLPqgkJ+9uYMIg2/NGsPN04f0yutEzjkKy4/zzo7DvLOjlJUF5dQ3tRAbFcG03HTOG9GPGcPTGTsgiche8vXXNTbzzo5SXl5fzLLth2loauHMQSncOTOXS8b17zVfp4SeR5fv4cevb2fFdy5kQLL3/3Cb2RrnXF6nllXw9Q77K47z3b9s4r1dZUwZnMJPPz+RkVmJXpfVbbUNzawoKOOdHa1HdR/5rtPl9uvDzNEZzByVwfTc9LC4FlZd18jL60t4/L0CPio/Tm6/Pnz1/FyumZwTFl+/BJe5C1dQebyRN/7tfK9LARR8Ycs5x1/WFXPfq1upqW/ia58dwb9cMCKkzr875ygoq+GdHaW8s+Mwq/ZW0NDUQnx0JOcMT+eC0RnMHJXJ4PTQuWfI35pbHG9sPsgjy/ewqbiKjMRYbj13KF+aNsTTa5cSPo7WNTL5vqV89fxcvj1rjNflAAq+sFd2rJ77XtnK4g0ljMzsyx0zh/O5CQOIjwnOo4IjNQ2sKGhtlPLuzjL2VbQe1Q3P6MMFozO5YHQGZw9N01FNO845Vuwp55F3C3h3Zyl9Y6O47TPDWDAzl4QYf7RbE+nYG5sPcMfTa3l2/nSm5aZ7XQ6g4BOfv20/xA9f3cbeshoSY6OYPSmbL549iAk5yZ42jqipb+LDwgo+2F3G33eXs+1gNc613hg+Pbf1qO6C0ZkMSgvfo7qu2lpSzcNv7+a1TQfonxTHty8fzZwzc3rltV7x3rdf2MiSzQdY+9+XBE3H1Ao++ZhzjlV7K3hu9X6WbD5AXWMLY/on8sWzB3HN5BxSEmJ6dPstLa2nLjcWVbKxqIoNRZVsKqqiqcURExnB5MEpnDuiH+eOSGfiwJSg+SUKVfmFFfzw1a1sKKrizIHJfP+qcZw1JLTusZLg5pxj+o+XkTckjYe/NMXrcj6m4JMOVdc1snh9Cc/l72djURUxkRFMGpTC6P6JjBmQyJj+iYzKSiQxruvXiZpbHGXH6imurGV/xXG2lFSzYX8lW0qqPx5FPj46kjOyk8gbmsa5I9LJG5IWtKdfQ1lLi+Ol9cX89I3tHKqu56ozs/n2rNEh1ZeiBK+tJdV87pfv8b/XTeQLecHTeYaCT05p24FqXlxTxPr9lew4eJSjvnACGJgaz5j+iQxIjicywogwIzICIiMiWp/NaGh2HKyqpaSyjpKqWg5W1f1TH5UxURGMHZDEmQOTmZCTzMSBKQzP6KObsAPoeEMTjywv4NHlewC484LhIdfYSYLPw2/v5mdv7uDD715EZlKc1+V8TMEnXeKco7iylh0Hj7L9xONANWXH6mlucbS41iO65hZHs2t9jo40+ifHkZ0cT3ZKPAOS48hOiSc7pfU5t19f/YENEsWVtfx4yTZe3XiAsQOSeOD6M0NmpGwJPtc/8gG1jc28+vXzvC7ln3Ql+NT0SzAzBqa2dsh80djOjfjgnFPvISEiJyWeX904hdlnHuS7f9nM7F+9z9cvHMm/fHa4rqlKl1Qdb2TtvkrunDnc61K6RT/1cloUeqHn0jP6s/Tfz+eKiQN48K2dXP3w39l+sNrrsiSEvLe7lOYWx2fHhNZoDO0p+ETCSGqfGH4xdzKP3HQWh6rruOqh93lo2S4am1u8Lk1CwNvbS0lJiGbSoFSvS+kWBZ9IGJo1vj9//feZXHZGfx5YupPP/+YD9rUZtkmkvZYWx/KdhzlvZEbI9xGr4BMJU2l9YvjVjVN4+MYpFJbVcMVD7/HmloNelyVBaktJNWXHGvjs6NA+zQkKPpGwd8XEAbx293kMTe/Dgt+v4f7XturUp3zC2zsOYwbnh9ho6x1R8IkIg9ISeOHOGdw8fQiPvbeXGxau5EBVrddlSRBZtv0wE3OS6dc31utSuk3BJyIAxEZF8sOrx/OLuZPYeqCaK375Pu/tKvW6LAkCH5XXsGF/JbPGD/C6FL/oVvCZWZqZLTWzXb7nDpv6mNk83zK7zGxem+nvmNkOM1vve2R2px4R6b45k3JYfNdn6Nc3hi8/8SEPLt1Jc0vodXQh/vPy+hIA5kzK9rgS/+juEd89wDLn3Ehgme/9PzGzNOBeYBowFbi3XUB+yTk3yffhE9MAAAAUJUlEQVQ43M16RMQPRmT25aWvncs1k3P4xbJdLPj9GmradGsn4cM5x0vripmem0Z2ivcjrftDd4NvDrDI93oRcHUHy1wGLHXOVTjnjgBLgVnd3K6I9LCEmCgeuP5M/mf2Gfxt+yGue2QFxZW67hduNhZVUVBWw9WTcrwuxW+6G3xZzrkDAL7njk5V5gD727wv8k074Xe+05z/beoORCSomBnzzhnK726dSlHFceb86u+s3XfE67IkgF5aX0xMZASXT+gd1/egE8FnZm+Z2eYOHnM6uY2OwuzEBYMvOecmAOf5Hjd/Sh3zzSzfzPJLS3XBXSSQZo7K4M//cg4JMZHMXbiSl9cXe12SBEBTcwuvbCjhorGZJMd3fbiyYHXK4HPOXeycG9/B42XgkJkNAPA9d3SNrghoO2jTQKDE99nFvuejwB9pvQZ4sjoWOufynHN5GRmhfx+JSKgZmZXIS187l0kDU/jXZ9bz86U7aVGjl17t/d1llB1rYE4vOs0J3T/VuRg40UpzHvByB8u8CVxqZqm+Ri2XAm+aWZSZ9QMws2jgSmBzN+sRkR6U1ieG398+levOGsgvl+3i68+so7ah2euypIe8vL6EpLiokO+Uur3uBt9PgEvMbBdwie89ZpZnZo8DOOcqgB8Cq32P+3zTYmkNwI3AeqAYeKyb9YhID4uNiuRn103kO5ePYcmmA3zp8ZUcqWnwuizxs5r6Jt7YfJArJmYTGxXpdTl+1a3x+Jxz5cBFHUzPB25v8/4J4Il2y9QAZ3Vn+yLiDTNjwczhDE5L4F+fXc91j3zAU7dNI6eXNHcXWLr1ELWNzVzdS+7da0s9t4jIabt8wgCe+spUDh+t59pfa3y/3uQv64rJSYnn7KFpXpfidwo+EemW6bnpPH/HDACuf2QFqwrKPa5Iuqv0aD3v7y5jzqRsIkJ8CKKOKPhEpNvG9E/ixTvPISMxlpuf+JA3Nh/wuiTphlc3ltDc4rhmcu9qzXmCgk9E/GJgagIv3nEOZ2Qncecf1vL7lR95XZKcppfWFTNuQBIjsxK9LqVHKPhExG9S+8Twx9unc+HoTP77pc38fOlOnNO9fqGkoPQYG4qqeu3RHij4RMTP4mMiefTms7jed6/f/7yyVTe6h5CX1pdgBrN7YWvOE7p1O4OISEeiIiP46ecn0jcuit/9vZCa+iZ+8vmJRPbChhK9yYmRGM4Znk5WUpzX5fQYBZ+I9IiICOP7V44jMS6aXy7bxfGGZh784iRionSiKVit21/JvorjfP3CEV6X0qMUfCLSY8yMb1wyisTYKO5fso2ahiZ+86WziI/pXT2B9BYvrSsmNiqCWeP7e11Kj9K/XiLS4756fi7/75oJLN9ZyrzffcjRukavS5J2GppaeHXjAS4el0ViXO8ZiaEjCj4RCYgbpw3m/744ibUfHeFLj69S/55B5uX1xVTUNPCFvEGnXjjEKfhEJGDmTMrh0ZvPYvvBo3xx4QoOH63zuiQBWlocj71XwJj+iZw/sp/X5fQ4BZ+IBNRFY7N48pazKTpSy9xHV3KgqtbrksLeOzsPs/PQMRbMzMWs97e8VfCJSMCdM6Lfx51bf+HRFeyvOO51SWHt0eUFZCfHceXE3nvvXlsKPhHxRN7QNP5w+zSqa5v4wqMrKCg95nVJYWn9/kpW7a3gK58ZRnRkeERCeHyVIhKUzhyUwp++Op2Gpha+uHAlOw8d9bqksLPw3T0kxkUxd+pgr0sJGAWfiHhqXHYSz8yfjgFzF65kS0mV1yWFjcKyGl7ffJCbpw+hb2z43NbdreAzszQzW2pmu3zPqSdZ7g0zqzSzV9tNH2Zmq3zrP2tmMd2pR0RC08isRJ5bMIO4qAhuWLiS9fsrvS4pLDz+fgHRERHccs5Qr0sJqO4e8d0DLHPOjQSW+d535GfAzR1M/ynwoG/9I8Bt3axHRELU0H59eHbBDFISYrjp8VV8uLfC65J6tbJj9TyfX8S1U3LI7MX9cnaku8E3B1jke70IuLqjhZxzy4B/OnlvrW1mLwReONX6IhIeBqUl8NyCGWQmxTLviQ95f1eZ1yX1Wk+t+Ij6phZuPy/X61ICrrvBl+WcOwDge87swrrpQKVzrsn3vgg46QBQZjbfzPLNLL+0tPS0CxaR4NY/OY5n589gSHoCX1m0mr9tP+R1Sb3O8YYmnlpRyCXjshiR2dfrcgLulMFnZm+Z2eYOHnO6ue2O7pI86aBdzrmFzrk851xeRkZGNzctIsEsIzGWP311OmP6JzL/qTUs2XTA65J6lefzi6g83siC88PvaA86EXzOuYudc+M7eLwMHDKzAQC+58Nd2HYZkGJmJ5oSDQRKuvoFiEjvlNonhqdvn8aZg1K4649r+cu6Iq9L6hWamlt47L0CzhqSSt7QNK/L8UR3T3UuBub5Xs8DXu7sis45B7wNXHc664tI75cUF81TX5nK9Nx0vvHcBv704T6vSwp5r28+SNGRWuaH6dEedD/4fgJcYma7gEt87zGzPDN7/MRCZvYe8DxwkZkVmdllvlnfBr5hZrtpveb3227WIyK9TJ/YKJ645WxmjsrgO3/exBPv7/W6pJDlnOPRd/eQ268Pl4zN8rocz3TrjkXnXDlwUQfT84Hb27w/7yTrFwBTu1ODiPR+cdGRPHrzWdz9p3Xc9+pWjjc08bXPjgiLDpX9afnOUjYXV/PjaycQERG++049t4hISIiNiuThG6dwzeQc/r+/7uT+17bResVEOqOusZkfLN7CsH59uHbKSRvQh4Xw6aNGREJeVGQED1x/Jsnx0Tz+/l6qahv58bUTiAqTzpW745HleygsP87vb5tKbFSk1+V4SsEnIiElIsK496pxJMdH84tlu6iua+QXcycTFx3ef8w/TWFZDb9+Zw9XThzAeSN1O5j+TRKRkGNm/Pslo7j3qnG8ueUQX3lyNcfqm069YhhyzvHfL28mJjKC/75ynNflBAUFn4iErFvPHcbPv3Amq/ZW8KXHVnKkpsHrkoLOkk0HeW9XGf9x6SiywqxPzpNR8IlISLt2ykAeveksth08yvWPruBAVa3XJQWNo3WN3PfqFs7ITuLm6UO8LidoKPhEJORdPC6LRbdO5WBVHZ//9QdsP1jtdUlB4cGluzh8tJ77r1EDoLa0J0SkV5gxPJ1nF0yn2Tmu+80K3t0Z3p3Zbymp4skP9nLj1MFMGpTidTlBRcEnIr3GGdnJvPS1cxmYGs+tT64O2y7OWloc//XSZtL6xPCty8Z4XU7QUfCJSK8yIDme5++YwWdG9OM7f97ET9/YTktLeN3o/szq/azbV8l3PzeW5IRor8sJOgo+Eel1EuOi+e28PG6cNpjfvLOHu59ZR11js9dlBUTZsXp++sZ2pg1L45rJ4d1Dy8noBnYR6ZWiIiO4/+rxDElL4Mevb+dAVR2PfTmPtD4xXpfWY5qaW/jWCxupqW/iR1ePV1+mJ6EjPhHptcyMBTOH8/CNU9hUXMXVD/+dLSVVXpfVI5xz/OCVLfxt+2HuvWocI7MSvS4paCn4RKTXu2LiAJ6ZP536pmau/fUHPLd6v9cl+d3Cdwt4euU+Fpyfy80zhnpdTlBT8IlIWJgyOJXX7j6PvKGpfOvFjXzz+Q3UNvSO636vbizhx69v54qJA/j2LLXiPBUFn4iEjX59Y3nqK9O4+8IRPL+miGt+/Xf2ltV4XVa35BdW8I3nNpA3JJUHrj8zrMfZ66xuBZ+ZpZnZUjPb5XtOPclyb5hZpZm92m76k2a218zW+x6TulOPiMipREYY37h0NL+79WwOVtcx+6H3eWPzAa/LOi0Fpce4/al8BqbE89iX8zRCRSd194jvHmCZc24ksMz3viM/A24+ybxvOucm+R7ru1mPiEinfHZ0Jq/dfR65mX254+m1/OjVrdQ3hc6pz7Jj9dzyu9VEmvG7W88mtRe3VvW37gbfHGCR7/Ui4OqOFnLOLQOOdnNbIiJ+lZMSz/MLZnDLOUN5/P29fO4X77GqoNzrsk6ptqGZ2xflc/hoHY/Py2NIeh+vSwop3Q2+LOfcAQDfc+ZpfMb9ZrbRzB40s9hu1iMi0iUxURH8YPYZLPrKVBqaW/jiwpV8+4WNVB4PziGOKmoaWPD0GjYUVfKLuZOZPLjDK0zyKU4ZfGb2lplt7uAxxw/b/w4wBjgbSAO+/Sl1zDezfDPLLy0N785nRcT/Zo7K4K//NpMFM3N5YW0RF/98OS+vL8a54Onu7IPdZVz+i3dZuaecH18zgcvO6O91SSHJuvNNNbMdwAXOuQNmNgB4xzk3+iTLXgD8p3PuytOZ31ZeXp7Lz88/7bpFRD7NlpIqvvvnTWwoquL8URncf/V4BqUleFZPY3MLP1+6k0eW72FYvz78cu5kxucke1ZPMDKzNc65vM4s291TnYuBeb7X84CXu7KyLyyx1n51rgY2d7MeEZFuOyM7mT//y7nce9U41hRWcMmDy/nx69s4WFUX8FoKy2q47jcf8Jt39jD37MG8+vXPKPS6qbtHfOnAc8BgYB9wvXOuwszygDucc7f7lnuP1lOafYFy4Dbn3Jtm9jcgAzBgvW+dY6faro74RCRQSipr+fHr23ltYwmREcbsM3P46vnDGNM/qUe365zjxbXF3PvyZqIiI/jJtRO4fMKAHt1mKOvKEV+3gs8rCj4RCbR95cd54u97eXb1fmobmzl/VAbzz8vl3BHpfu0MuqXFsbqwgkUrClmy6SDThqXx4BcnkZ0S77dt9EYKPhGRHlJ5vIGnV37Ekx98RNmxesYOSOILeQM5e2gaYwckEXkaPac459hSUs3iDSW8sqGEA1V1xEVHcNdnR3DnBSNO6zPDjYJPRKSH1TU28/L6Yh5/by+7DrdeoekbG8XkwSmcPTSNvKGpTBqUQkLMJ0d/a25x1DU2c6Cqjtc2HmDxhmL2lNYQFWHMHJXB7EnZXDw2iz6xGjmusxR8IiIBVFxZS35hBasLK8gvPMKOQ0dxDqIijIGp8TQ2twZdfVMLdY3NNLUbEX7asDRmT8rmc+MHqAeW09SV4NO/EyIi3ZSTEk/OpBzmTGod8bzqeCNr9x1hdWEF+4/UEhsVQWxUBHHRkf/0nBgXzQWjM3T9LsAUfCIifpacEM1nx2Ty2TGn05mV9DQNSyQiImFFwSciImFFwSciImFFwSciImFFwSciImFFwSciImFFwSciImFFwSciImElJLssM7NS4KMe3kw/oKyHtxFqtE86pv3ySdonHdN+6Zg/9ssQ51xGZxYMyeALBDPL72y/b+FC+6Rj2i+fpH3SMe2XjgV6v+hUp4iIhBUFn4iIhBUF38kt9LqAIKR90jHtl0/SPumY9kvHArpfdI1PRETCio74REQkrCj4REQkrCj4fMzsejPbYmYtZnbSZrVmNsvMdpjZbjO7J5A1BpqZpZnZUjPb5XtOPclyzWa23vdYHOg6A+VU33szizWzZ33zV5nZ0MBXGVid2Ce3mFlpm5+P272oM5DM7AkzO2xmm08y38zsl759ttHMpgS6Ri90Yr9cYGZVbX5Wvt9TtSj4/mEzcC3w7skWMLNI4GHgcmAccIOZjQtMeZ64B1jmnBsJLPO970itc26S7zE7cOUFTie/97cBR5xzI4AHgZ8GtsrA6sLvw7Ntfj4eD2iR3ngSmPUp8y8HRvoe84HfBKCmYPAkn75fAN5r87NyX08VouDzcc5tc87tOMViU4HdzrkC51wD8Awwp+er88wcYJHv9SLgag9r8Vpnvvdt99cLwEVmZgGsMdDC7fehU5xz7wIVn7LIHOAp12olkGJmAwJTnXc6sV8CRsHXNTnA/jbvi3zTeqss59wBAN9z5kmWizOzfDNbaWa9NRw7873/eBnnXBNQBaQHpDpvdPb34fO+U3ovmNmgwJQW1MLt70hXzDCzDWb2upmd0VMbieqpDw5GZvYW0L+DWd9zzr3cmY/oYFpI3w/yafukCx8z2DlXYma5wN/MbJNzbo9/Kgwanfne97qfj1PozNf7CvAn51y9md1B6xHxhT1eWXALt5+TzlpLa3+bx8zsc8BLtJ4O9ruwCj7n3MXd/IgioO1/rAOBkm5+pqc+bZ+Y2SEzG+CcO+A7FXP4JJ9R4nsuMLN3gMlAbwu+znzvTyxTZGZRQDJBcmqnh5xynzjnytu8fYxeft2zk3rd3xF/cM5Vt3m9xMx+bWb9nHN+79Rbpzq7ZjUw0syGmVkMMBfota0Yaf3a5vlezwM+cVRsZqlmFut73Q84F9gasAoDpzPf+7b76zrgb6539xBxyn3S7trVbGBbAOsLVouBL/tad04Hqk5cUghnZtb/xDVxM5tKaz6Vf/papyesjvg+jZldAzwEZACvmdl659xlZpYNPO6c+5xzrsnM7gLeBCKBJ5xzWzwsu6f9BHjOzG4D9gHXA/hu97jDOXc7MBZ41MxaaP1B/YlzrtcF38m+92Z2H5DvnFsM/Bb4vZntpvVIb653Ffe8Tu6Tu81sNtBE6z65xbOCA8TM/gRcAPQzsyLgXiAawDn3CLAE+BywGzgO3OpNpYHVif1yHXCnmTUBtcDcnvrHUV2WiYhIWNGpThERCSsKPhERCSsKPhERCSsKPhERCSsKPhERCSsKPhERCSsKPhERCSsKPpEQZGY/NbN/NbPNZrbazMZ6XZNIqFDwiYQYMzsHuBRYDxQD/wP8wtOiREKIgk8k9EylddQDAxqBN4CzPK1IJIQo+ERCT0f9DDYHvAqREKXgEwk97wFX4OvgF/i8b5qIdIJGZxAJMc65tWb2Aq3j26XROu7fTd5WJRI6NDqDSIgyswuA/3TOXel1LSKhRKc6RUQkrOiIT0REwoqO+EREJKwo+EREJKwo+EREJKwo+EREJKwo+EREJKz8/23ES97CcMA7AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(7,4))\n", "plt.sympy_function(f(φ, m=1), x_values=(-1.05, 1.5) )\n", "plt.xlabel(\"φ\")\n", "plt.title(\"Bulk free energy\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compared to last tutorial, this bulk free energy has also two minima, but at different values. \n", "\n", "Another complexity of the model is its anisotropy. The gradient parameter $\\epsilon$ depends on the interface normal." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$\\bar{\\epsilon} \\left(δ \\cos{\\left (j \\left(- θ_{0} + \\operatorname{atan_{2}}{\\left ({\\partial_{1} {{φ}_{C}}},{\\partial_{0} {{φ}_{C}}} \\right )}\\right) \\right )} + 1\\right)$$" ], "text/plain": [ "\\bar{\\epsilon}⋅(δ⋅cos(j⋅(-θ₀ + atan2(D(phi_C), D(phi_C)))) + 1)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def σ(θ):\n", " return 1 + δ * sp.cos(j * (θ - θzero))\n", "\n", "θ = sp.atan2(ps.fd.Diff(φ, 1), ps.fd.Diff(φ, 0))\n", "\n", "ε_val = εb * σ(θ)\n", "ε_val" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def m_func(T):\n", " return (α / sp.pi) * sp.atan(γ * (Teq - T))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, we just insert these parameters into the free energy formulation before doing the functional derivative, to make the dependence of $\\epsilon$ on $\\nabla \\phi$ explicit." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/latex": [ "$${{φ}_{C}}^{3} - \\frac{{{φ}_{C}}^{2} α}{\\pi} \\operatorname{atan}{\\left ({{T}_{C}} γ - T_{eq} γ \\right )} - \\frac{3 {{φ}_{C}}^{2}}{2} + \\frac{{{φ}_{C}} α}{\\pi} \\operatorname{atan}{\\left ({{T}_{C}} γ - T_{eq} γ \\right )} + \\frac{{{φ}_{C}}}{2} - \\bar{\\epsilon}^{2} δ^{2} \\cos^{2}{\\left (j θ_{0} - j \\operatorname{atan_{2}}{\\left ({\\partial_{1} {{φ}_{C}}},{\\partial_{0} {{φ}_{C}}} \\right )} \\right )} {\\partial_{0} {\\partial_{0} {{φ}_{C}}}} - \\bar{\\epsilon}^{2} δ^{2} \\cos^{2}{\\left (j θ_{0} - j \\operatorname{atan_{2}}{\\left ({\\partial_{1} {{φ}_{C}}},{\\partial_{0} {{φ}_{C}}} \\right )} \\right )} {\\partial_{1} {\\partial_{1} {{φ}_{C}}}} - 2 \\bar{\\epsilon}^{2} δ \\cos{\\left (j θ_{0} - j \\operatorname{atan_{2}}{\\left ({\\partial_{1} {{φ}_{C}}},{\\partial_{0} {{φ}_{C}}} \\right )} \\right )} {\\partial_{0} {\\partial_{0} {{φ}_{C}}}} - 2 \\bar{\\epsilon}^{2} δ \\cos{\\left (j θ_{0} - j \\operatorname{atan_{2}}{\\left ({\\partial_{1} {{φ}_{C}}},{\\partial_{0} {{φ}_{C}}} \\right )} \\right )} {\\partial_{1} {\\partial_{1} {{φ}_{C}}}} - \\bar{\\epsilon}^{2} {\\partial_{0} {\\partial_{0} {{φ}_{C}}}} - \\bar{\\epsilon}^{2} {\\partial_{1} {\\partial_{1} {{φ}_{C}}}}$$" ], "text/plain": [ " 2 2 \n", " 3 φ_C ⋅α⋅atan(T_C⋅γ - T_eq⋅γ) 3⋅φ_C φ_C⋅α⋅atan(T_C⋅γ - T_eq⋅γ) φ_C\n", "φ_C - ─────────────────────────── - ────── + ────────────────────────── + ───\n", " π 2 π 2 \n", "\n", " \n", " 2 2 2 \n", " - \\bar{\\epsilon} ⋅δ ⋅cos (j⋅θ₀ - j⋅atan2(D(phi_C), D(phi_C)))⋅D(D(phi_C)) - \\\n", " \n", "\n", " \n", " 2 2 2 \n", "bar{\\epsilon} ⋅δ ⋅cos (j⋅θ₀ - j⋅atan2(D(phi_C), D(phi_C)))⋅D(D(phi_C)) - 2⋅\\ba\n", " \n", "\n", " \n", " 2 \n", "r{\\epsilon} ⋅δ⋅cos(j⋅θ₀ - j⋅atan2(D(phi_C), D(phi_C)))⋅D(D(phi_C)) - 2⋅\\bar{\\e\n", " \n", "\n", " \n", " 2 \n", "psilon} ⋅δ⋅cos(j⋅θ₀ - j⋅atan2(D(phi_C), D(phi_C)))⋅D(D(phi_C)) - \\bar{\\epsilon\n", " \n", "\n", " \n", " 2 2 \n", "} ⋅D(D(phi_C)) - \\bar{\\epsilon} ⋅D(D(phi_C))\n", " " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fe = free_energy_density.subs({\n", " m: m_func(T),\n", " ε: εb * σ(θ),\n", "})\n", "\n", "dF_dφ = ps.fd.functional_derivative(fe, φ)\n", "dF_dφ = ps.fd.expand_diff_full(dF_dφ, functions=[φ])\n", "dF_dφ" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we insert all the numeric parameters and discretize:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA78AAAAXBAMAAAAxaeW3AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZpkQ3Ynvq81UMrtEdiLw+n06AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJa0lEQVRoBe1ZX4ybRxH/+b6zfT77fM5DQ1/SOClHigrECiX0pbpToUpBCNyql6KiNFeVf21ROQLkUKVwriIUgUD1C43EAzF/qgr14SyEKP9ELIQi8RQrVQFVSmJFah+TpkfbNGl7zMzufPvtet2LBEelknnYmZ35zex45vvWuzaA3DZco/diBSY/3OCPldzxXvxw1z4TVaDyMJch9/JQMarzX2hnlJN1M5l4bvKDwK757wLnbv4D6eaf70CU4mDMtz3/IYv55YH5+WYyc+/AeLux+tqzy1eWnmg6jS8lM9d1fE0wU4DywCzTJxtOm993n5sEksaIJepDd901OgrUKNxUx/f2ZwpfH0l+UlLfP5259KlmMjGaF+d/T5ibGJdf4NGjcyhdFEXpQWbjXYClwtpaE8kDeGmQLGJ/D7f2xlpGKQ5iRh3FtmCwdW1tbXF3rfqoF5sm+SYKcxirhXpZhJQT7eSzQ7asQgGGmySzdpbf5xTJD80ndZqMpLFiiToYPy7fxy0DpwkkNTI31QkA/tTC10Hah5RL6ru7maYv3A3VOo5Rfe9hYKTBZ9v4F5s+vmOV2d+7Rpq49wxQXECxXmygsIBfo9yFKMVBpFIN1ZZg8BEC4yBwgGNkqQhcqCPXy+pYtsthM/D50HZ0acnFUYBw9Qo93Lxcxx/dLJA01sFIogpNDiw3kFtAvqWakKtRuKlOCMnOFf6OSFmVvKSkWe+srOkLd8PUG5ieA/YyNNLgE73kkkSZ4gZXr+8CLFFjgOkW8i/TC1i8XHmb56IUB5HyXSRzgkETuA8UZz/DsvQzYGsbU42sTmRZDngWODZszKAVYLj1ytgD8eQgUGSnGiuWqMOdbaC8iJI8707rJDUKl+o4W0RS+DpIWpVIShoJIipNX7gbKl/DbGt0g+mpMVu01C4/0c00eLaLqcsTq9Tg8gIvIW0VB5EqVxr5mmDIVuknV4CzHcZlqAPQG1HJaKxoW/U6cLI9bHUaBRi+boN/6zyHJRsrmqhDU6mnF1F5w2l8SY3CpTq+PZgpfB2kabCUNAiQTr1SyEQ1vEWnb/AUfV2u1clrfNG43tIULrW7IW3wplMdFOgNfpORq9M/2EfnqSIrAXIw0vLFu1JMGfgqvcE9sgf0RX9ulzWtSl6hBvd9APCxJUrY5KcAy2MNTk53OYCJe+mpUz2eBSRGjTUqUZsZlXq2jspbQQieCkKNysdHvut72Edh7B5DCohejQaDwSWNEIM0feHPceFkoOrJ5UgaXGzRcfhXZ7ac5xi7aUL04u0mNNcu6WuDy43qm7Q/Y4Jf75PN2cdRaEOU4mCk/NrhFEN50f78RE1iZocqpZIlu6xpVfVVYKWZNZNMew6TABVgeazBY3gghSeXanhIvP3BizUiUS0IlfpCE5OU2RBJGDUqPxnmr24JPfF0BHGxYkgDShvMJR0mAXmluJ8LJ0MT1XtuJJet7PZkB5Pt/OAFL0Te1JNrl8dE12zRhPgN8CP8mHeqPZhdxRidrUUJ40Dmo399rWcxyct0Iu6X9rcZ5FEuvtfZBr8SafCFTBB+PFaadDgwPNbgG+BWSNYa+EDDW99NNNaIRBXIDa7HGywQNSq3r6C6E5/aRxkrKYzmw0gFpQ2WkjptVtL0hd/P9ZCBl9rbo+rXgWfkMlzGN7N+wJelIly7n2QbvHOA3OmnLgNjdT6m5XizBimtw85BbhFnH7SYqRbpt7ywv0PMpyLFjZBplW47HmB5xm1TCrA81uCGOQGaEHQj2NrzormJxhqRqAKpwZltVbWOq9Fyqk5I1+GoUylc6ujUoUSrEpmShjYz1/SFe1s0CvzyfasPnOiQ8DTmaFQ6Dxxv84RqlzTTBn/Hfp9ye07Rt0cLuYsQpTiIVBhgkl8exvD9mWhZ0hRRBzo9xsi2io4Jx2T5DOaRj2aiKMDwWIORz6xwJzW4lonliRqLlJFEFUqlpoNRyW0LarBcjZZTdUL6DKbdJ1C41DFEurlpsJbU6TOSpi/cDdUBypTsxCJBiy0atiH7Tq01XINLS0vLX+nLNYm+PnZykuRWquN8cYHfYFGKg0izZP+GwWCWwwNHZPQG2qJiZFtFl9atrhgGSLuPIwUYHm0wb05KB9/hDeYbsi4WSVRDnKVrUhe5bJHUJFyNhnN1PDOdZC7zgUVJ4VJHVQ5zWpVISzpsJ42mL9wN06vSYPkO5mc9WYVstjbGI8DnBiyb2hW6RqoDdyL3KFYG+DlwZpK+gxcgSnEQiR43/NlgsEKazYPJtzFEdA2mc/HpPnad/lvGaFu12f7GlrH4T4kCDI82eNZVEyv0HZyN5clprGiiCqVS048T9MPdCFKj4VydAEhVGK+lOoVLHVPtkGAaLCUdsllFmj7XTCYyUMvGqe57GcU/dEzN4UvisbvFbBNKb8lF2dRu2jb4HHd36nDyOKp3zO+Yw7exuw1RioNIlcPI9QVD+2wd+Gfj1pqE4rgpHW/QQ3U4WUw+bdpvlpXnabaF8XbycOhzogHysYdaBQhXrzS4CMcYbeHFZvJQGC81aqx4oplTND6F3R37A5AsYQeTuxqZS3UCJL0N433ysAdkC48hHYivSRRGShqES0GavnA35Go4RsulDS7XqFtMcqlDZfvNPRbHDl6i7TZ//PWaSNWZA1S2u7f3MEHX5jnktz9GR1lWioMxPz1Dl2PGAC+1yXk7Tcc/wbFTyv9u7UgfxUOb6Kg2ucBqs6wsQmKy7x+d0KdyaIsDGoAFGi9/BeD9jLZxcfeOoXipUWPFEk1BO75+pEY7Dn3cT1IFfDK5i3EPb0uP2er4SHqQC1RxPTVbuKmjj3QgWZViSklHgDR94W7AX7b/iQKlDeagG0k/jQQv9Ojvi/bUYsQkqpjPKCwQom8fgoaIIQAr1gWVBlG/iNJH0nfwdC+CYpWP/E9Aoe//rMHNcGWaT3fQmB2U5bGOmBHzieGMzkf/oko3uYB8RGDU6bqgvCLX5QHSO0X7zgHSN+rsqkAKVi4NjvwfrPb/Fk/6kUjldrJY6Kx0IiZWRX1GYEN09VX6Aymgq4q3Pui2IOzoaYA8lb0H+14B0jfq7KpAClZ+kwjf0+mG8VwscnLodKNyaH/MxLqozyhwiD60bQh5VfHWB9WGAo9SBMj8zI1XiYzDgnBxUKC1v++WhqsRADdwev0Gxv5/Dz0503j3S2AO7+9+Htcy2JgKjF3qbEzga1FtBf4NbgI9dnSphTMAAAAASUVORK5CYII=\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.02, \\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.02, θ₀: 0.2, κ: 1.8, τ: 0.0003}" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "discretize = ps.fd.Discretization2ndOrder(dx=0.03, dt=1e-5)\n", "parameters = {\n", " τ: 0.0003,\n", " κ: 1.8,\n", " εb: 0.01,\n", " δ: 0.02,\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", "φ_eqs = ps.simp.sympy_cse_on_assignment_list([ps.Assignment(φ_delta_field.center, \n", " discretize(dφ_dt.subs(parameters)))])\n", "φ_eqs.append(ps.Assignment(φ, discretize(ps.fd.transient(φ) - φ_delta_field.center)))\n", "\n", "temperature_evolution = -ps.fd.transient(T) + ps.fd.diffusion(T, 1) + κ * φ_delta_field.center\n", "temperature_eqs = [\n", " ps.Assignment(T, discretize(temperature_evolution.subs(parameters)))\n", "]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When creating the kernels we pass as target (which may be 'cpu' or 'gpu') the default target of the target handling. This enables to switch to a GPU simulation just by changing the parameter of the data handling.\n", "\n", "The rest is similar to the previous tutorial." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "φ_kernel = ps.create_kernel(φ_eqs, cpu_openmp=4, target=dh.default_target).compile()\n", "temperature_kernel = ps.create_kernel(temperature_eqs, target=dh.default_target).compile()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def timeloop(steps=200):\n", " φ_sync = dh.synchronization_function(['phi'])\n", " temperature_sync = dh.synchronization_function(['T'])\n", " dh.all_to_gpu() # this does nothing when running on CPU\n", " for t in range(steps):\n", " φ_sync()\n", " dh.run_kernel(φ_kernel)\n", " temperature_sync()\n", " dh.run_kernel(temperature_kernel)\n", " dh.all_to_cpu()\n", " return dh.gather_array('phi')\n", " \n", "def init(nucleus_size=np.sqrt(5)):\n", " for b in dh.iterate():\n", " x, y = b.cell_index_arrays\n", " x, y = x-b.shape[0]//2, y-b.shape[0]//2\n", " bArr = (x**2 + y**2) < nucleus_size**2\n", " b['phi'].fill(0)\n", " b['phi'][(x**2 + y**2) < nucleus_size**2] = 1.0\n", " b['T'].fill(0.0)\n", " \n", "def plot():\n", " plt.subplot(1,3,1)\n", " plt.scalar_field(dh.gather_array('phi'))\n", " plt.title(\"φ\")\n", " plt.colorbar()\n", " plt.subplot(1,3,2)\n", " plt.title(\"T\")\n", " plt.scalar_field(dh.gather_array('T'))\n", " plt.colorbar()\n", " plt.subplot(1,3,3)\n", " plt.title(\"∂φ\")\n", " plt.scalar_field(dh.gather_array('phidelta'))\n", " plt.colorbar()" ] }, { "cell_type": "code", "execution_count": 12, "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| (nan,nan)| (nan,nan)\n", "\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7kAAAF1CAYAAAAtEi0mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+0ZWV95/n3x+KHRg2/NTRgQ1rsjjHdqDXITKZtWwTRlVh2BzsYWys9MCzT0knGZFZwnGhCzFqS7rQdZ5ikMRLRMQLBpK2kSdOIMklnRawiIgo0oSRpucASsJCoaSRV9Z0/9r73Hm6dW/cW55577t3P+7XWXvfsH+ec55zCj+e7n2c/O1WFJEmSJElD8IxZN0CSJEmSpLVikStJkiRJGgyLXEmSJEnSYFjkSpIkSZIGwyJXkiRJkjQYFrmSJEmSpMGwyJUkSZIkDYZFriSpaUm+NbLsT/LfR9bfMuv2SZKkQ2ORK0lqWlU9Z34Bvgr88Mi2j8+6fZI0iSQvS3Jfkq8meees2yOtB4tcHVSS/ynJF5O8MskjSW5L8uJZt0uSJEmr8ufAGcA24N1J/mGSE5P8pyRvSHJXkvuT/NMZt1NaMxa5WlaS7wJ+B7gCeAvwceAa4ONJMsu2SZIkaWVV9a2q+quq+gLdb7nXAv8P8AhwFPAkcD5wdZKTZtdSae1Y5OpgzgQCfAh4NvAN4P8CfgB4wQzbJUmSpFVIcmGSW5NcDewDngn8EPDLwJHAN6rqVuCLwGtm11Jp7Vjk6mCeDzxQVTW/oaqeAB4DvmdmrZIkSdKKkvwA8D7gjcC/By6i68E9DLh/yeEP4e87DYRFrg7mAeDk0aHJSZ4FHAPMzaxVkiRJWo3XAL9fVQ8BO4Fv01169iQHjso7CX/faSAscnUwt9KF4Tvohi1vAX4R+JOqemCWDZMkSdKK/oru9xvAe4A/qqr/Rlfo/jxwBECSbcDfA26YRSOltXbYrBugjauq/qYPvauAl9Fdx/FfgLfOtGGSJElajY8D/yTJl4F7gf+l3/5TdHOu/ArwLOAE4J9W1WMzaaW0xjJyuaW0rCT/L7C7qn5h1m2RJEnS5JJcBPzzqnrVrNsirSWHK0uSJEmSBsMiV9JUJbkqycP9UKlx+5Pkg0l2J7kjycvWu42SNKkk5yW5p8+yS8fsf2WSP0uyN8n5S/ZtT3Jvv2wf2f7yJF/qX/OD3qNe0qxtlqxzuLKkqUrySuBbwEer6iVj9r8e+FfA64FXAL9WVa9Y31ZK0tOXZAvw58A5dLPT7gTeXFV3jRxzKvDdwM8CO6rq+n77scAuYCtQwG3Ay6vqsSSfp7t28nN0EwJ9sKr+cJ0+liQ9xWbKOntyJU1VVf0RsOcgh2yjK4Crqj4HHJ3kxPVpnSStiTPp5q24r6qepJu5dtvoAVX1l1V1B7B/yXNfC9xUVXv6SX9uAs7rc/C7q+pP+/vVf5TuXqeSNCubJussciXN2kk89Yb0c/02SdosJsmx5Z679J6lZqOkWds0WbchbiF0/PHH16mnnjrrZkib2m233fZoVZ1wKM957T9+dn19z77J3veO79wJPDGy6cqquvIQXmLcdReDvI7CrJMmt0GzbpIcW+65mzYbzTppcmbdZDZEkXvqqaeya9euWTdD2tSS/LdDfc6je/Zx640nT/S+h5/4lSeqausELzEHnDKyfjLw4ESN2qDMOmlyGzTrJsmxOeBVS557S7/95CXbN0U2mnXS5My6yThcWdKs7QDe1s+yfBbweFU9NOtGSdIh2AmcnuS0JEcAF9Bl22rcCJyb5JgkxwDnAjf2OfjNJGf1M42+DfjUNBovSau0abJuQ/TkSpqVYl8tnRdgbSX5BN2Zu+OTzAHvBQ4HqKrfoJtF7/XAbuCvgX8x1QZJatB0s66q9ia5hO5H3Bbgqqq6M8llwK6q2pHkfwB+DzgG+OEkv1hV319Ve5L8Et2PR4DLqmp+sr6fAD4CPAv4w36RpGWYdfMscqWGFbB/ypd4VdWbV9hfwDum2ghJTVunrLuB7qTd6Lb3jDzeyVOH5I0edxVw1Zjtu4ADbr0mSeOYdYsscqXG7T9ghndJGh6zTlILzLqO1+RKkiRJkgbDnlypYUWxrzbFHSkk6Wkz6yS1wKxbZJErNW7a125I0kZg1klqgVnXsciVGlbAPsNQ0sCZdZJaYNYtssiVGucZP0ktMOsktcCs6zjxlCRJkiRpMOzJlRpW4AQFkgbPrJPUArNukUWu1DjvpiapBWadpBaYdR2LXKlhRTlBgaTBM+sktcCsW2SRK7WsYJ9ZKGnozDpJLTDrFjjxlCRJkiRpMOzJlRpWeO2GpOEz6yS1wKxbZJErNS3sI7NuhCRNmVknqQVm3TyLXKlhBez32g1JA2fWSWqBWbfIa3IlSZIkSYNhT67UOIe1SGqBWSepBWZdxyJXalhhGEoaPrNOUgvMukUWuVLj9pdhKGn4zDpJLTDrOha5UsM84yepBWadpBaYdYuceEqSJEmSNBj25EoNK8I+z3VJGjizTlILzLpFFrlS47x2Q1ILzDpJLTDrOha5UsO8dkNSC8w6SS0w6xZZ5EpNC/vKYS2Shs6sk9QCs26e34IkSZIkaTDsyZUaVsB+z3VJGjizTlILzLpFFrlS47x2Q1ILzDpJLTDrOha5UsOqvHZD0vCZdZJaYNYt8luQJEmSJA2GPblS4/Y7rEVSA8w6SS0w6zoWuVLDuvupOaBD0rCZdZJaYNYt8luQmtZduzHJIkkb3/SzLsl5Se5JsjvJpWP2H5nk2n7/rUlO7be/JcntI8v+JGf0+27pX3N+3/PW+IuRNChm3Tx7cqWGOdW8pBZMO+uSbAGuAM4B5oCdSXZU1V0jh10IPFZVL0xyAXA58KNV9XHg4/3r/ADwqaq6feR5b6mqXVNrvKTBMOsW+etWkiRpMmcCu6vqvqp6ErgG2LbkmG3A1f3j64Gzkyy9eO7NwCem2lJJevo2TdbZkys1bl85QYGk4VuDrDs+yWgvw5VVdWX/+CTg/pF9c8Arljx/4Ziq2pvkceA44NGRY36UA38w/laSfcAngfdVVU32MSQNmVnXsciVGlbECQokDd4aZd2jVbV1mX3jflUu/YF20GOSvAL466r68sj+t1TVA0meS/fD763ARw+hzZIaYtYt8tet1Lj99YyJFknaDKacdXPAKSPrJwMPLndMksOAo4A9I/svYMnwvap6oP/7TeC36YYKStKyzLqOPblSw5xqXlIL1iHrdgKnJzkNeIDuR9yPLTlmB7Ad+FPgfOAz88PxkjwDeBPwyvmD+x+HR1fVo0kOB34I+PQ0P4Skzc2sW2SRK0mSNIH+urNLgBuBLcBVVXVnksuAXVW1A/gw8LEku+l6NS4YeYlXAnNVdd/ItiOBG/sffVvofvR9aB0+jiSNtZmyziJXalgRJ56SNHjrkXVVdQNww5Jt7xl5/ARdD8a4594CnLVk27eBl695QyUNllm3yCJXapz3yZXUArNOUgvMuo5FrtSwKtjn5FGSBs6sk9QCs26R34IkSZIkaTDsyZWaFvaPvZ2ZJA2JWSepBWbdPItcqWGFw1okDZ9ZJ6kFZt0ii1ypcd4nV1ILzDpJLTDrOha5UsOKsN9bCEkaOLNOUgvMukWW+pIkSZKkwVixyE3yzCSfT/LFJHcm+cV++2lJbk1yb5JrkxzRbz+yX9/d7z91uh9B0iT28YyJlqEw66RhM+s6Zp00bGZdZzWf5DvAq6vqHwBnAOclOQu4HPhAVZ0OPAZc2B9/IfBYVb0Q+EB/nKQNqID99YyJlgEx66SBMuuewqyTBsqsW7TiJ6nOt/rVw/ulgFcD1/fbrwbe2D/e1q/T7z87iYPDpQ0p7JtwGQqzThoys26eWScNmVk3b1UTTyXZAtwGvBC4AvgK8I2q2tsfMgec1D8+CbgfoKr2JnkcOA54dMlrXgxcDPBMvotznvGmyT6J1LjncszLD/U582f81Jl21m055hhO+7VfnfbHkAbtiFNONusmNO2se8ELXjDtjyBpDLNu0aq+haraV1VnACcDZwLfN+6w/u+4UwB1wIaqK6tqa1VtPZwjV9teSZqaaWfdluc8e+0aK0lP07Sz7oQTTli7xkrS03BItxCqqm8kuQU4Czg6yWH9Wb+TgQf7w+aAU4C5JIcBRwF71q7JktbSkIamrBWzThoes+5AZp00PGZdZzWzK5+Q5Oj+8bOA1wB3A58Fzu8P2w58qn+8o1+n3/+ZqjrgjJ+k2auKExT0zDppuMy6RWadNFxm3aLV9OSeCFzdX7/xDOC6qvqDJHcB1yR5H/AF4MP98R8GPpZkN92Zvgum0G5Ja2TflAMtyXnArwFbgN+sqvcv2f8CuklNju6PubSqbphqo8Yz66QBm3bWbSJmnTRgZl1nxSK3qu4AXjpm+31013Es3f4E4CxSkuYnN7kCOIduyNvOJDuq6q6Rw/5Puh9Zv57kxcANwKnr3VazTlILzDpJLTika3IlDUsB+6d77caZwO7+xxNJrqG7HcVokVvAd/ePj2LxOjBJWhPrkHWSNHNm3SKLXKlpmfawloVbT/TmgFcsOeYXgP+c5F8Bz6a7PkyS1tDUs06SNgCzbp5FrtSw7n5qE5/xOz7JrpH1K6vqyv7xam498WbgI1X1q0n+R7prv15SVfsnbZgkwZplnSRtaGbdIotcqXH7Vne77IN5tKq2LrNv/tYT80ZvSzHvQuA8gKr60yTPBI4HHp60YZI0bw2yTpI2PLOu47cgaZp2AqcnOS3JEXSzcu5YcsxXgbMBknwf8EzgkXVtpSRJkgbDnlypYUWmOqylqvYmuQS4ke72QFdV1Z1JLgN2VdUO4GeADyX53+hG2vy492CUtJamnXWStBGYdYsscqXG7Z/ygI7+nrc3LNn2npHHdwE/ONVGSGretLNOkjYCs65jkSs1rAr2ecZP0sCZdZJaYNYtssiVGuewFkktMOsktcCs69ifLUmSJEkaDHtypYZ1ExR4rkvSsJl1klpg1i2yyJUatw+HtUgaPrNOUgvMuo5FrtSwwms3JA2fWSepBWbdIvuzJUmSJEmDYU+u1DSv3ZDUArNOUgvMunl+C1Lj9pOJFknaDKaddUnOS3JPkt1JLh2z/8gk1/b7b01yar/91CT/Pcnt/fIbI895eZIv9c/5YBJDV9JBmXUde3KlhnnTcEktmHbWJdkCXAGcA8wBO5PsqKq7Rg67EHisql6Y5ALgcuBH+31fqaozxrz0rwMXA58DbgDOA/5wSh9D0iZn1i2yJ1dq3P56xkSLJG0GU866M4HdVXVfVT0JXANsW3LMNuDq/vH1wNkH661IciLw3VX1p1VVwEeBNz6dzy6pHWZdx1+okiRJKzs+ya6R5eKRfScB94+sz/XbGHdMVe0FHgeO6/edluQLSf6/JP9w5Pi5FV5TktbaILLO4cpSw7qbhjtcWdKwrVHWPVpVW5fZN+7Fa5XHPAS8oKq+nuTlwH9I8v2rfE1JWmDWLbLIlRrn5FGSWjDlrJsDThlZPxl4cJlj5pIcBhwF7OmH530HoKpuS/IV4EX98Sev8JqS9BRmXcfhylLD5m8aPskiSRvdOmTdTuD0JKclOQK4ANix5JgdwPb+8fnAZ6qqkpzQT+ZCku8FTgfuq6qHgG8mOau/nu1twKfW5AuRNEhm3SJ7ciVJkiZQVXuTXALcCGwBrqqqO5NcBuyqqh3Ah4GPJdkN7KH7cQjwSuCyJHuBfcDbq2pPv+8ngI8Az6KbadSZlSXNzGbKOotcqXHOkCypBdPOuqq6ge7WF6Pb3jPy+AngTWOe90ngk8u85i7gJWvbUklDZtZ1LHKlljnkWFILzDpJLTDrFljkSg0rnHhK0vCZdZJaYNYtssiVGucZP0ktMOsktcCs63gxniRJkiRpMOzJlRo2P9W8JA2ZWSepBWbdIotcqXGGoaQWmHWSWmDWdSxypYYVzsInafjMOkktMOsWWeRKjXMWPkktMOsktcCs6zjxlCRJkiRpMOzJlVpWXrshqQFmnaQWmHULLHKlhjkLn6QWmHWSWmDWLbLIlRpnGEpqgVknqQVmXcdrciVJkiRJg2FPrtQwp5qX1AKzTlILzLpFFrlS48owlNQAs05SC8y6jkWu1DjvpyapBWadpBaYdR2LXKlh5VTzkhpg1klqgVm3yImnJEmSJEmDYU+u1Div3ZDUArNOUgvMuo5FrtQ0Z+GT1AKzTlILzLp5FrlS4zzjJ6kFZp2kFph1HYtcqWGFExRIGj6zTlILzLpFTjwlSZIkSRoMe3KlllU33bwkDZpZJ6kFZt0Ci1ypcd40XFILzDpJLTDrOha5UsMKJyiQNHxmnaQWmHWLvCZXkiRJkjQY9uRKTfN+apJaYNZJaoFZN88iV2qcExRIaoFZJ6kFZl3HIldqnNduSGqBWSepBWZdxyJXaliVYShp+Mw6SS0w6xY58ZQkSZIkaTAscqXG7a9MtEjSZjDtrEtyXpJ7kuxOcumY/Ucmubbff2uSU/vt5yS5LcmX+r+vHnnOLf1r3t4vz1vDr0TSAJl1HYcrS41zggJJLZhm1iXZAlwBnAPMATuT7Kiqu0YOuxB4rKpemOQC4HLgR4FHgR+uqgeTvAS4EThp5Hlvqapd02u9pCEx6zr25EqNq8pEiyRtBlPOujOB3VV1X1U9CVwDbFtyzDbg6v7x9cDZSVJVX6iqB/vtdwLPTHLkGn1sSY0x6zoWuVLDismC0CJX0mawRll3fJJdI8vFI29xEnD/yPocT+2heMoxVbUXeBw4bskxPwJ8oaq+M7Ltt/rhez+fxNCVtCyzbpHDlSVJklb2aFVtXWbfuB9kSwcNHvSYJN9PN6zv3JH9b6mqB5I8F/gk8Fbgo6tvsiQdskFk3Yo9uUlOSfLZJHcnuTPJT/Xbj01yU5J7+7/H9NuT5IP9xcZ3JHnZJA2UNF014TIUZp00bFPOujnglJH1k4EHlzsmyWHAUcCefv1k4PeAt1XVVxbaXPVA//ebwG/TDRWciFknDZtZ11nNcOW9wM9U1fcBZwHvSPJi4FLg5qo6Hbi5Xwd4HXB6v1wM/PqkjZQ0JeU1uSPMOmmopp91O4HTk5yW5AjgAmDHkmN2ANv7x+cDn6mqSnI08B+Bd1XVn8wfnOSwJMf3jw8Hfgj48sTfhVknDZdZt2DFIreqHqqqP+sffxO4m26s9ehFxVcDb+wfbwM+Wp3PAUcnOXHShkqaErtyAbNOGrwpZl1/3dkldLOF3g1cV1V3JrksyRv6wz4MHJdkN/BOFovIS4AXAj+/5PYZRwI3JrkDuB14APjQhN+CWScNnVkHHOI1uenuc/RS4Fbg+VX1EHSBOXI/o+UuSH5oyWtdTHdGkGfyXU+j6ZI2gyTnAb8GbAF+s6reP+aYfwb8Al28frGqfmxdG3lge05lClm35ZhjptpuSbNTVTcANyzZ9p6Rx08AbxrzvPcB71vmZV++lm1calpZ94IXvGCazZY0Q5sl61Zd5CZ5Dt2FwD9dVX91kEmvVnNBMlV1JXAlwHfn2AH1B0mbyzSHHGcV91NLcjrwLuAHq+qxrMENwCcxzaw78gWnmHXSjAzs8oqJTTPrtm7datZJM2LWdVZ1C6F+fPQngY9X1e/2m782P1yl//twv301FyRL2iCqJltWsJr7qf2vwBVV9VjXnnqYGTHrpOGactZtKmadNFxmXWc1syuHbmz13VX1b0d2jV5UvB341Mj2t/Wz8Z0FPD4//EXSxlKsyQQFk95P7UXAi5L8SZLP9cOb151ZJw3XGmXdIJh10nCZdYtWM1z5B+nuVfSlJLf32/4P4P3AdUkuBL7K4tjrG4DXA7uBvwb+xZq2WNLaKWDyQJv0fmqH0c3a+Sq6HoI/TvKSqvrGpA07RGadNFRrk3VDYdZJQ2XWLVixyK2q/8L4H6oAZ485voB3TNguScOw2vupfa6q/gb4iyT30BW9O9eniR2zTlILzDpJLVjVNbmShmvK126s5n5q/wH4xwD9fdJeBNy3tp9SUuu8Tk1SC8y6ziHdQkjSAE0x0Kpqb5L5+6ltAa6av58asKuqdvT7zk1yF7AP+N+r6uvTa5WkJg3ox5skLcusAyxypcZNf5KBVdxPrehuFv7OqTZEUsOGNaGKJI1n1s2zyJVa5xk/SS0w6yS1wKwDvCZXkiRJkjQg9uRKLSsc1iJp+Mw6SS0w6xZY5Eqtc1iLpBaYdZJaYNYBFrmSlr1doiQNiVknqQVmHXhNriRJkiRpQOzJlVrnsBZJLTDrJLXArAMsciUZhpJaYNZJaoFZB1jkSm0rwFn4JA2dWSepBWbdAotcqXHlGT9JDTDrJLXArOs48ZQkSZIkaTDsyZVa5xk/SS0w6yS1wKwDLHIlee2GpBaYdZJaYNYBFrlS8+IZP0kNMOsktcCs61jkSi0rHNYiafjMOkktMOsWOPGUJEmSJGkw7MmVmhav3ZDUALNOUgvMunkWuVLrHNYiqQVmnaQWmHWARa4kw1BSC8w6SS0w6wCvyZUkSZIkDYg9uVLrPOMnqQVmnaQWmHWARa7UtsIJCiQNn1knqQVm3QKLXKlx3jRcUgvMOkktMOs6XpMrta4mXCRpM5hy1iU5L8k9SXYnuXTM/iOTXNvvvzXJqSP73tVvvyfJa1f7mpJ0ALMOsCdXq3Tjg7cvPH7t3zpjhi2RJGljSbIFuAI4B5gDdibZUVV3jRx2IfBYVb0wyQXA5cCPJnkxcAHw/cDfAj6d5EX9c1Z6TUlaN5sp6+zJlSRJmsyZwO6quq+qngSuAbYtOWYbcHX/+Hrg7CTpt19TVd+pqr8Advevt5rXlKT1tGmyziJXalxqskWSNoMpZ91JwP0j63P9trHHVNVe4HHguIM8dzWvKUlPYdZ1HK6sVXGI8oA5C5+kFkyedccn2TWyfmVVXdk/HvfiS38uLnfMctvHdUR4alHSwZl1gEWu1DYnj5LUgrXJukerausy++aAU0bWTwYeXOaYuSSHAUcBe1Z47kqvKUmLzLoFDleWJEmazE7g9CSnJTmCbnKVHUuO2QFs7x+fD3ymqqrffkE/I+lpwOnA51f5mpK0njZN1tmTK7XOnlxJLZhi1lXV3iSXADcCW4CrqurOJJcBu6pqB/Bh4GNJdtP1alzQP/fOJNcBdwF7gXdU1T6Aca85vU8haRDMOsAiV2qek0dJasG0s66qbgBuWLLtPSOPnwDetMxzfxn45dW8piQdjFnXsciVWmeRK6kFZp2kFph1gNfkSpIkSZIGxJ5cqXWe8ZPUArNOUgvMOsAiV2raKm/8LUmbmlknqQVm3SKLXKl1k980XJI2PrNOUgvMOsAiV5Jn/CS1wKyT1AKzDnDiKUmSJEnSgNiTKzXOazcktcCsk9QCs65jkSu1zjCU1AKzTlILzDrAIldqm7PwSWqBWSepBWbdAq/JlSRJkiQNhj25Uus84yepBWadpBaYdYBFriTDUFILzDpJLTDrAItcqXleuyGpBWadpBaYdR2vyZUkSZIkDYZFriRJkiRpMByuLLXOYS2SWmDWSWqBWQdY5Ept835qklpg1klqgVm3wCJXap1hKKkFZp2kFph1gEWuJMNQUgvMOkktMOsAJ56SJEmSJA2IPblSw4LXbkgaPrNOUgvMukUWuVLrDENJLTDrJLXArAMscqW2OQufpBaYdZJaYNYt8JpcSVOV5Lwk9yTZneTSgxx3fpJKsnU92ydJkqRhsSdXat0Uz/gl2QJcAZwDzAE7k+yoqruWHPdc4CeBW6fXGklNs3dDUgvMOmAVPblJrkrycJIvj2w7NslNSe7t/x7Tb0+SD/Y9Nnckedk0Gy9pDdSEy8GdCeyuqvuq6kngGmDbmON+CfgV4InJPszTZ9ZJAzfdrNtUzDtpwMw6YHXDlT8CnLdk26XAzVV1OnBzvw7wOuD0frkY+PW1aaakaUlNtgDHJ9k1slw88vInAfePrM/12xbfP3kpcEpV/cG0P+sKPoJZJw3WGmTdkHwE804aJLOus2KRW1V/BOxZsnkbcHX/+GrgjSPbP1qdzwFHJzlxrRoraQomP+P3aFVtHVmuHHn1LPOO3c7kGcAHgJ9Z2w916Mw6aeDs3Vhg3kkDZtYBT3/iqedX1UMA/d/n9dtX7LWZl+Ti+Z6fv+E7T7MZkja4OeCUkfWTgQdH1p8LvAS4JclfAmcBOzbQ5FNrmnX7vvXtqTZWkiYwUd6NZt0jjzwy9cZK0sGs9ezKB+21ecrGqivne34O58g1boakVZn0bN/KZ/x2AqcnOS3JEcAFwI6Ft696vKqOr6pTq+pU4HPAG6pq1xp9wml5Wlm35TnPnnKzJI01/awbslXl3WjWnXDCCevQLEkHMOsWPN0i92vzQ1X6vw/321fqtZG0wUzz2o2q2gtcAtwI3A1cV1V3JrksyRum/+kmZtZJA+F1aisy76QBMOs6T7fI3QFs7x9vBz41sv1t/Ux8ZwGPzw99kbRBTfmMX1XdUFUvqqq/U1W/3G97T1XtGHPsqzZYL65ZJw3FDHs3lpu5eMxx2/tj7k2yvd/2XUn+Y5L/muTOJO8fOf7HkzyS5PZ+uWiCZpp30hCYdcDqbiH0CeBPgb+bZC7JhcD7gXOS3Et3/8v5RtwA3AfsBj4E/MuVXl/SbHnGr2PWScM246xbbubixfYlxwLvBV5Bd/u19478QPw3VfX3gJcCP5jkdSNPvbaqzuiX31xNY8w7abjMus5hKx1QVW9eZtfZY44t4B0rvaYkbTRmnaQp2ga8qn98NXAL8HNLjnktcFNV7QFIchNwXlV9AvgsQFU9meTP6IYMP23mnaQp2TBZt9YTT0nabGY4rEWS1s3kWXewe4KvZLmZi0et5r7iRwM/TNdDMu9HktyR5Poko9fOSmqRWQesoidX0oBZqEpqwdpk3aNVteztzZJ8GvieMbvevcrXX+m+4ocBnwA+WFX39Zt/H/hEVX0nydvpek5evcr3kzQ0Zt0Ci1ypYWF80kjSkKxH1lXVa5Z9/+RrSU6sqoeWzFw8ao7FYX7QDdO7ZWT9SuDeqvp3I+/59ZH9HwIufxpNlzQQZt0ihytLkiRN13IzF4+6ETg3yTH9JCwKWF4mAAAUTUlEQVTn9ttI8j7gKOCnR58wf8uf3hvobtUmSbOyYbLOnlypdQ5XltSC2Wbd+4Hr+lmMvwq8CSDJVuDtVXVRVe1J8kvAzv45l/XbTqYbBvhfgT9LAvB/97OL/mR/z/G9wB7gx9fzQ0nagMw6wCJXat6QbgMkScuZZdb1Q+3GzVy8C7hoZP0q4Kolx8yxzAjEqnoX8K41baykTc2s61jkSq2zyJXUArNOUgvMOsAiV5JhKKkFZp2kFph1gBNPSZIkSZIGxJ5cqWXlNbmSGmDWSWqBWbfAIldqnWEoqQVmnaQWmHWARa7UPM/4SWqBWSepBWZdxyJXap1hKKkFZp2kFph1gBNPSZIkSZIGxJ5cqXEOa5HUArNOUgvMuo5FrtSywmEtkobPrJPUArNugUWu1DrDUFILzDpJLTDrAK/JlSRJkiQNiD25UsOC125IGj6zTlILzLpFFrlS6wxDSS0w6yS1wKwDLHKl5qVMQ0nDZ9ZJaoFZ17HIlVrmLHySWmDWSWqBWbfAiackSZIkSYNhT67UOCcokNQCs05SC8y6jkWu1DrDUFILzDpJLTDrAItcqXme8ZPUArNOUgvMuo5FrtQ6w1BSC8w6SS0w6wAnnpIkSZIkDYg9uVLLymEtkhpg1klqgVm3wCJXap1hKKkFZp2kFph1gEWu1LTgGT9Jw2fWSWqBWbfIa3IlSZIkSYNhT67UuvKUn6QGmHWSWmDWARa5UvMc1iKpBWadpBaYdR2LXKllhRMUSBo+s05SC8y6BRa5UuOyf9YtkKTpM+sktcCs6zjxlCRJkiRpMOzJlVrnsBZJLTDrJLXArAPsyZWal5pskaTNYJZZl+TYJDclubf/e8wyx23vj7k3yfaR7bckuSfJ7f3yvH77kUmuTbI7ya1JTp2spZI2O7OuY5ErtazoppqfZJGkjW72WXcpcHNVnQ7c3K8/RZJjgfcCrwDOBN675AfiW6rqjH55uN92IfBYVb0Q+ABw+aQNlbSJmXULLHKlxtmTK6kFM866bcDV/eOrgTeOOea1wE1VtaeqHgNuAs47hNe9Hjg7SSZuraRNy6zrWORKkiSt7Pgku0aWiw/huc+vqocA+r/PG3PMScD9I+tz/bZ5v9UP3/v5kR93C8+pqr3A48Bxh9AuSVpqEFnnxFNS6+yNldSCybPu0arautzOJJ8GvmfMrnev8vXH9UrMt/otVfVAkucCnwTeCnx0hedIapFZB1jkSk0LDjmWNHzrkXVV9Zpl3z/5WpITq+qhJCcCD485bA541cj6ycAt/Ws/0P/9ZpLfpruO7aP9c04B5pIcBhwF7Jn800jajMy6RQ5Xllo26eQETjwlaTOYfdbtALb3j7cDnxpzzI3AuUmO6SdhORe4MclhSY4HSHI48EPAl8e87vnAZ6oMZqlZZt0Ce3IlSZKm6/3AdUkuBL4KvAkgyVbg7VV1UVXtSfJLwM7+OZf1255N9wPwcGAL8GngQ/0xHwY+lmQ3Xa/GBev3kSTpABsm6yxypcY5XFlSC2aZdVX1deDsMdt3AReNrF8FXLXkmG8DL1/mdZ+g/xEpSWDWzbPIlVpnkSupBWadpBaYdYBFrtQ8e3IltcCsk9QCs65jkSu1rID9pqGkgTPrJLXArFvg7MqSJEmSpMGwJ1dqnSf8JLXArJPUArMOsMiVmue1G5JaYNZJaoFZ17HIlVo3+Y2/JWnjM+sktcCsA7wmV2pearJlxddPzktyT5LdSS4ds/+dSe5KckeSm5P87Wl8Tkltm3bWSdJGYNZ1LHIlTU2SLcAVwOuAFwNvTvLiJYd9AdhaVX8fuB74lfVtpSRJkobEIldqWa3BcnBnArur6r6qehK4Btj2lCZUfbaq/rpf/Rxw8qQfS5KeYvpZJ0mzZ9Yt8JpcqWEBMt1rN04C7h9ZnwNecZDjLwT+cJoNktSedcg6SZo5s26RRa7Uuv0Tv8LxSXaNrF9ZVVf2jzPm+LHpm+SfA1uBfzRxiyRpqcmzTpI2PrMOsMiVNLlHq2rrMvvmgFNG1k8GHlx6UJLXAO8G/lFVfWftmyhJkqRWWORKjZvysJadwOlJTgMeAC4Afuwp75+8FPj3wHlV9fA0GyOpXQ7hk9QCs64zlYmnVrpliKQNYsoTFFTVXuAS4EbgbuC6qrozyWVJ3tAf9q+B5wC/k+T2JDvW8BNOnXknbQJOxjIxs07aBMy6BWvekztyy5Bz6IYq7kyyo6ruWuv3kjSpmvpNw6vqBuCGJdveM/L4NVNtwBSZd9JmMf2sGzKzTtoszLp50+jJXfGWIZI2Dm8aPhHzTtokzLqJmHXSJmHWdaZR5I67ZchJSw9KcnGSXUl2/Q3OMyNpU1ox70azbt+3vr2ujZOkNXJIWffII4+sa+MkaalpTDy1qluG9LcYuRJg69atddOu35lCU6R2JLntaT3RYS2TWDHvlmbdrp/6mfVolzRY+emfNevW3yFn3Xo0StIYZh0wnSJ3VbcMkbQBFMT7qU3CvJM2A7NuUmadtBmYdQumMVx54ZYhSY6gu2XIppotVWpK1WRL28w7abMw6yZh1kmbhVkHTKEnt6r2Jpm/ZcgW4KqqunOt30eSZs28k9QCs07SZjON4cpjbxkiaYMazkm7mTDvpE3CrJuIWSdtEmYdMKUiV9LmkQENTZGk5Zh1klpg1nUscqXWGYaSWmDWSWqBWQdY5EptK8BZ+CQNnVknqQVm3YJpzK4sSZIkSdJM2JMrNSyU125IGjyzTlILzLpFFrlS6wxDSS0w6yS1wKwDLHIlGYaSWmDWSWqBWQdY5Eptc4ICSS0w6yS1wKxb4MRTkiRJkqTBsMiVGpeqiRZJ2gxmmXVJjk1yU5J7+7/HLHPc9v6Ye5Ns77c9N8ntI8ujSf5dv+/Hkzwysu+iiRoqadMz6zoOV5ZaZ6EqqQWzzbpLgZur6v1JLu3Xf270gCTHAu8FttINOrwtyY6qegw4Y+S424DfHXnqtVV1ybQ/gKRNwqwD7MmVGlddGE6ySNKGN/Os2wZc3T++GnjjmGNeC9xUVXv6H3s3AeeNHpDkdOB5wB9P2iBJQ2TWzbPIlSRJWtnxSXaNLBcfwnOfX1UPAfR/nzfmmJOA+0fW5/pto95M15sx+kv0R5LckeT6JKccQpskaZxBZJ3DlaWWFfbGShq+tcm6R6tq63I7k3wa+J4xu969ytfPmG1LG30B8NaR9d8HPlFV30nydrqek1ev8v0kDY1Zt8AiV2qdU81LasGUs66qXrPcviRfS3JiVT2U5ETg4TGHzQGvGlk/Gbhl5DX+AXBYVd028p5fHzn+Q8DlT6/1kgbDrAMcriw1z9mVJbVgxlm3A9jeP94OfGrMMTcC5yY5pp+R9Nx+27w3A594ymfqfkTOewNw96QNlbS5mXUde3Kl1lmoSmrBbLPu/cB1SS4Evgq8CSDJVuDtVXVRVe1J8kvAzv45l1XVnpHX+GfA65e87k8meQOwF9gD/PgUP4OkzcCsAyxyJUmSpqofanf2mO27gItG1q8CrlrmNb53zLZ3Ae9au5ZK0tO3kbLOIldqWQH77cmVNHBmnaQWmHULLHKlpnmvW0ktMOsktcCsm2eRK7XOMJTUArNOUgvMOsAiV5JhKKkFZp2kFph1gLcQkiRJkiQNiD25UsucoEBSC8w6SS0w6xZY5EpNK6j9s26EJE2ZWSepBWbdPItcqXVeuyGpBWadpBaYdYDX5EqSJEmSBsSeXKllXrshqQVmnaQWmHULLHKl1jmsRVILzDpJLTDrAItcSYahpBaYdZJaYNYBFrlS48owlNQAs05SC8y6eU48JUmSJEkaDHtypZYVsN/7qUkaOLNOUgvMugUWuVLrHNYiqQVmnaQWmHWARa4kw1BSC8w6SS0w6wCLXKlx5f3UJDXArJPUArNunhNPSZIkSZIGw55cqWUFVU5QIGngzDpJLTDrFljkSq1zWIukFph1klpg1gEWuZKcoEBSC8w6SS0w6wCvyZUkSZIkDYg9uVLLqrxpuKThM+sktcCsW2CRK7XOYS2SWmDWSWqBWQdY5ErNK8/4SWqAWSepBWZdxyJXalp5xk9SA8w6SS0w6+Y58ZQkSZIkaTDsyZVaVng/NUnDZ9ZJaoFZt8AiV2pdee2GpAaYdZJaYNYBFrlS0wooz/hJGjizTlILzLpFXpMrtayqO+M3ybKCJOcluSfJ7iSXjtl/ZJJr+/23Jjl1Cp9UUsvWIesOJsmxSW5Kcm//95hljvtPSb6R5A+WbD+tz8d7+7w8ot9ufkpaZNYtsMiVNDVJtgBXAK8DXgy8OcmLlxx2IfBYVb0Q+ABw+fq2UpKm7lLg5qo6Hbi5Xx/nXwNvHbP9cuAD/fMfo8tNMD8lbSwbJusscqXG1f6aaFnBmcDuqrqvqp4ErgG2LTlmG3B1//h64OwkWdMPKal5U866lYzm3NXAG8e2sepm4Juj2/o8fDVdPi59vvkp6SnMuo7X5Eqtm+4EBScB94+szwGvWO6Yqtqb5HHgOODRaTZMUmNmOxnL86vqIYCqeijJ8w7huccB36iqvf36HF1ugvkpaSmzDtggRe5tt932rST3zLodvePZOP/nYFvGsy3j/d1DfcI3eezGT9f1x0/4vs9Msmtk/cqqurJ/PO4s29LThKs5ZhBuu+22R5N8m43z38xG+u/XtoxnWw70t5NcPJIzK1qHrCPJp4HvGfO8d0/4vgfLyA2Zn2bdQdmW8WzLgcy6CbJuQxS5wD1VtXXWjQBIssu2HMi2jLfR2nKoz6mq86bRlhFzwCkj6ycDDy5zzFySw4CjgD1TbtdMVNUJG+2/GdtyINsy3kZrC7DqH37rkHVU1WuW25fka0lO7Hs2TgQePoSXfhQ4OslhfQ/HaI5uyPw065ZnW8azLeOZdU8/67wmV9I07QRO72fLOwK4ANix5JgdwPb+8fnAZ6pq5j0RkrSGRnNuO/Cp1T6xz8PP0uXj0uebn5I2kg2TdRa5kqamPxN3CXAjcDdwXVXdmeSyJG/oD/swcFyS3cA7WX4mPknarN4PnJPkXuCcfp0kW5P85vxBSf4Y+B26SVXmkry23/VzwDv7nDyOLjfB/JS0sWyYrNsow5VX3Q2/DmzLeLZlPNuygqq6Abhhybb3jDx+AnjTerdrhjbSv5NtGc+2jGdbnqaq+jpw9pjtu4CLRtb/4TLPv49utvql2zdyfm6kfyPbMp5tGc+2PE0bKeviqBZJkiRJ0lA4XFmSJEmSNBgzL3KTnJfkniS7k6z7tSRJ/jLJl5LcPj87bZJjk9yU5N7+7zFTeu+rkjyc5Msj28a+dzof7L+nO5K8bB3a8gtJHui/m9uTvH5k37v6ttwzMo5+rdpySpLPJrk7yZ1Jfqrfvu7fzUHasu7fTZJnJvl8ki/2bfnFfvtpSW7tv5dr+wmeSHJkv76733/qWrVFh86sM+vGtMWsG98Ws24TM+vMujFtMevGt8Wsm6aqmtkCbAG+AnwvcATwReDF69yGvwSOX7LtV4BL+8eXApdP6b1fCbwM+PJK7w28HvhDuvtEnQXcug5t+QXgZ8cc++L+3+pI4LT+33DLGrblROBl/ePnAn/ev+e6fzcHacu6fzf953tO//hw4Nb+814HXNBv/w3gJ/rH/xL4jf7xBcC10/jv2GVV/3ZmnVk3ri1m3fi2mHWbdDHrzLpl2mLWjW+LWTfFZdY9uWcCu6vqvqp6ErgG2DbjNkHXhqv7x1cDb5zGm1TVH3HgPZ6We+9twEer8zm6+0idOOW2LGcbcE1Vfaeq/gLYzZiLxCdoy0NV9Wf942/Szcp7EjP4bg7SluVM7bvpP9+3+tXD+6WAVwPX99uXfi/z39f1dDPYjbuZtqbPrDPrxrXFrBvfFrNu8zLrzLpxbTHrxrfFrJuiWRe5JwH3j6zPcfD/0KahgP+c5LYkF/fbnl9VD0H3PwbgeevYnuXee1bf1SX9UJGrRob3rFtb+qEYL6U7uzXT72ZJW2AG302SLUlup7u59k10ZxS/Ud2tepa+30Jb+v2P003HrvVn1h3IrBth1h3QBrNuczLrDmTWjTDrDmiDWTclsy5yx519WO/pnn+wql4GvA54R5JXrvP7r9YsvqtfB/4OcAbwEPCr69mWJM8BPgn8dFX91cEOnXZ7xrRlJt9NVe2rqjOAk+nOJH7fQd5vI/zvS52N8G9h1i3PrFu+LWadDsVG+Lcw65Zn1i3fFrNuYGZd5M4Bp4ysnww8uJ4NqKoH+78PA79H9x/Y1+aHRfR/H17HJi333uv+XVXV1/r/8e0HPsTi8IyptyXJ4XTh8/Gq+t1+80y+m3FtmeV307//N4Bb6K7dODrJ/D2vR99voS39/qNY/dAlrS2z7kBmHWbdSsy6TcesO5BZh1m3ErNu7c26yN0JnN7PInYE3UXUO9brzZM8O8lz5x8D5wJf7tuwvT9sO/Cp9WrTQd57B/C2dM4CHp8f4jEtS65/+Cd03818Wy7oZ3k7DTgd+Pwavm+ADwN3V9W/Hdm17t/Ncm2ZxXeT5IQkR/ePnwW8hu5aks8C5/eHLf1e5r+v84HPVJVn/GbDrDuQWWfWLdcWs27zMusOZNaZdcu1xaybpprxzFd0M6j9Od0Y9Hev83t/L92MaV8E7px/f7rx7TcD9/Z/j53S+3+CbkjE39CdnblwufemG6JwRf89fQnYug5t+Vj/XnfQ/Q/rxJHj39235R7gdWvclv+ZbvjFHcDt/fL6WXw3B2nLun83wN8HvtC/55eB94z8d/x5uskQfgc4st/+zH59d7//e6f9vymXg/77mXVm3dK2mHXj22LWbeLFrDPrxrTFrBvfFrNuikv6L02SJEmSpE1v1sOVJUmSJElaMxa5kiRJkqTBsMiVJEmSJA2GRa4kSZIkaTAsciVJkiRJg2GRK0mSJEkaDItcSZIkSdJgWORKkiRJkgbj/wfRXnNz8z/exgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "timeloop(10)\n", "init()\n", "plot()\n", "print(dh)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result = None\n", "if 'is_test_run' not in globals():\n", " ps_notebook.set_display_mode('video')\n", " ani = ps.plot2d.scalar_field_animation(timeloop, rescale=True, frames=300)\n", " result = ps_notebook.display_animation(ani) \n", "result" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "assert np.isfinite(dh.max('phi'))\n", "assert np.isfinite(dh.max('T'))" ] } ], "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 }