{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from pystencils.session import *\n", "import pystencils as ps\n", "sp.init_printing()\n", "frac = sp.Rational" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial 06: 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=ps.Target.CPU)\n", "φ_field = dh.add_array('phi', latex_name='φ')\n", "φ_field_tmp = dh.add_array('phi_temp', latex_name='φ_temp')\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": [ "$\\displaystyle \\frac{{φ}_{(0,0)}^{4}}{4} - {φ}_{(0,0)}^{3} \\cdot \\left(\\frac{1}{2} - \\frac{m}{3}\\right) + {φ}_{(0,0)}^{2} \\cdot \\left(\\frac{1}{4} - \\frac{m}{2}\\right) + \\frac{ε^{2} \\left({\\partial_{0} {φ}_{(0,0)}}^{2} + {\\partial_{1} {φ}_{(0,0)}}^{2}\\right)}{2}$" ], "text/plain": [ " 4 2 ⎛ 2 2⎞\n", "φ_C 3 ⎛1 m⎞ 2 ⎛1 m⎞ ε ⋅⎝D(φ[0,0]) + D(φ[0,0]) ⎠\n", "──── - φ_C ⋅⎜─ - ─⎟ + φ_C ⋅⎜─ - ─⎟ + ────────────────────────────\n", " 4 ⎝2 3⎠ ⎝4 2⎠ 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", "φ_tmp = φ_field_tmp.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": "iVBORw0KGgoAAAANSUhEUgAAAboAAAEWCAYAAAAQKVIQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0DElEQVR4nO3deXyU5b338c8vewJZyAqENewICBhZtIp1t1pQqy1uRasFPbU95/Q5p9s5rT22Ptb2sZ4utorWirWtexUVtZQqboBE9n0JERK2LCSB7Mv1/DGDTWOQhEzmnpl836/XvDJzLzO/3IR8c9/XdV+XOecQERGJVFFeFyAiItKTFHQiIhLRFHQiIhLRFHQiIhLRFHQiIhLRFHQiIhLRFHQip8jMnJmN9D9/3Mx+3Mn9cszsbTM7amb392yVIhLjdQEiXjGzIiAHaAGagPeB251z+3r4o+cDZUCK042sIj1OZ3TS233eOdcXGAAcAn4VhM8cCmw5UciZWdj9AWo++n0iIUk/mCKAc64eeA4Yf3yZmb1lZre1eX2zmb17svcys2Qze9PMfmlm1m7d48A84FtmdszMLjSzH5rZc2b2pJlVAzebWaqZ/c7MDphZiZn92Myi27zPV8xsq5kdMbM3zGzop9Qzw8zeN7NKM1tvZue1+x5/ZGbv+S+l/tXMMruw7z1m9h5QC+SZ2cVmtt3MqszsN2a23MxuM7M4M6sws4lt9s82s1ozyzrZMRXpDgWdCGBmScCXgJXdfJ8MYBnwnnPuG+3P2pxzNwN/BH7qnOvrnPubf9UcfEGb5l//ONAMjASmABcDt/k/Yw7wPeBqIAt4B/jzCerJBV4FfgykA/8BPN8uXK4HbgGygTj/Np3d9yZ8l2KTgSr/9/BdIAPYDpzl/74bgaeAG9vsex2wzDlX2lHtIoGioJPe7kUzq8T3S/oi4GfdeK+BwHLgWefcf3dx3xXOuRedc61ACvA54N+cczXOucPAA8Bc/7a3A/c657Y655qB/wtMPsFZ3Y3AEufcEudcq3NuKVDgf//jfu+c2+GcqwOeASZ3Yd/HnXOb/XVcBmx2zr3gf/1L4GCbbRcB17U5y70J+EMXj5NIlynopLe70jmXBiQAdwLLzaz/Kb7X5UAi8NAp7Nu2A8xQIBY44L9kWAk8jO+M6/j6X7RZVwEYkNvB+w4Frj2+rX/7z+BrkzyubRjVAn27sG/buge2fe0/my1u83qV//3PM7Ox+M5WF3d4NEQCKOwavUV6gnOuBXjBzB7G98v8OaAGSGqz2ckC8BGgH7DEzC51ztV0pYQ2z/cBDUCm/8yovX3APc65P3biffcBf3DOfbULtXRl37Z1HwAGHX/hP3Mb1G77RfjOFA8Cz/nbRkV6lM7oRPi41+AcfEG11b94HXC1mSX575e7tRNvdSe+tqmXzSzxVGpxzh0A/grcb2YpZhZlZiPMbJZ/k4eA75rZaf7aU83s2hO83ZPA583sEjOLNrMEMzvPzNoHUCD2fRWYaGZX+nuOfo1P/nHwJHAVvrB7ohM1iHSbgk56u5fN7BhQDdwDzHPObfavewBoxHfbwSJ8nUQ+lf9y3Xx8l+xeMrOEU6zry/g6hmwBjuA7wxzg/4y/APcBT/l7aW7C1z7WUT378HV0+R5Qiu8s7T/pxP/9ru7rnCsDrgV+CpTj68FagO/stO17rsF3JvjOyWoQCQTT/aoi0hP899UVAzc4595ss/wxYP8pdNgROSVqoxORgDGzS4BVQB2+sz+jzS0bZjYM320RU7yoT3onXboUkUCaCezGN8TZ5/H1aq0DMLMf4bvM+jPn3B7vSpTeRpcuRUQkoumMTkREIlpYttFlZma6YcOGeV2GiIiEiA8//LDMOdfhuKlhGXTDhg2joKDA6zJERCREmNlHJ1qnS5ciIhLRFHQiIhLRFHQiIhLRFHQiIhLRFHQiIhLRFHQiIhLRFHQiIhLRemXQHayq5/6/buej8q7MiykiIuGoVwZdY3Mrv/r7LpbvKPW6FBER6WG9MugGpycyMDWBlYXlXpciIiI9rFcGnZkxIy+DlYUVaPYGEZHI1iuDDmDGiAwqahrZceiY16WIiEgP6rVBNzMvA0CXL0VEIlyvDbpB/RLJTUtU0ImIRLheG3T/aKcrp7VV7XQiIpGq1wYdwIy8dI7UNrHj8FGvSxERkR4SkKAzs0vNbLuZ7TKz73Sw/ptmtsXMNpjZMjMb2mbdPDPb6X/MC0Q9nTXjeDvdbl2+FBGJVN0OOjOLBh4ELgPGA9eZ2fh2m60F8p1zk4DngJ/6900H7gKmA9OAu8ysX3dr6qzB6UkM6pfICrXTiYhErECc0U0DdjnnCp1zjcBTwJy2Gzjn3nTO1fpfrgQG+Z9fAix1zlU4544AS4FLA1BTp83Iy2DVngq104mIRKhABF0usK/N62L/shO5FXitq/ua2XwzKzCzgtLSwA3dNSMvg8raJrYfUjudiEgkCmpnFDO7EcgHftbVfZ1zC51z+c65/KysrIDVNCMvHYAVaqcTEYlIgQi6EmBwm9eD/Mv+iZldCPwXMNs519CVfXvSoH5JDE7X/XQiIpEqEEG3GhhlZsPNLA6YCyxuu4GZTQEexhdyh9usegO42Mz6+TuhXOxfFlQzhqudTkQkUnU76JxzzcCd+AJqK/CMc26zmd1tZrP9m/0M6As8a2brzGyxf98K4Ef4wnI1cLd/WVDNHJFBVV0T2w6qnU5EJNLEBOJNnHNLgCXtlv2gzfMLP2Xfx4DHAlHHqZruv59uRWE54wemeFmKiIgEWK8eGeW43LREhqQnqZ1ORCQCKej8ZuSl84Ha6UREIo6Czu94O92WA9VelyIiIgGkoPObPlzz04mIRCIFnd/AtESGZiSxsjDonT5FRKQHKejamJmXwao95bSonU5EJGIo6NqYkZfB0fpmtqqdTkQkYijo2vh4fjq104mIRAwFXRv9UxMYntlHAzyLiEQQBV07x++nUzudiEhkUNC1MyMvg6MNzWzZr3Y6EZFIoKBrZ8bH416WeVyJiIgEgoKunZyUBEZm92X5jsDNYi4iIt5R0HXggnHZrCqsoLq+yetSRESkmxR0HbhwXA7NrY63dVYnIhL2FHQdmDqkH/2SYvnblkNelyIiIt0UkKAzs0vNbLuZ7TKz73Sw/lwzW2NmzWZ2Tbt1Lf5Zxz+eedxr0VHGZ8dk8+b2UppbWr0uR0REuqHbQWdm0cCDwGXAeOA6MxvfbrO9wM3Anzp4izrn3GT/Y3Z36wmUC8fnUFXXxIcfHfG6FBER6YZAnNFNA3Y55wqdc43AU8Ccths454qccxuAsDk9OmdUJrHRxrJth70uRUREuiEQQZcL7Gvzuti/rLMSzKzAzFaa2ZUBqCcgkhNimZGXwd+2qp1ORCSchUJnlKHOuXzgeuB/zWxERxuZ2Xx/IBaUlganN+QFY7MpLK2hsPRYUD5PREQCLxBBVwIMbvN6kH9ZpzjnSvxfC4G3gCkn2G6hcy7fOZeflZV16tV2wQXjcgBYtlWXL0VEwlUggm41MMrMhptZHDAX6FTvSTPrZ2bx/ueZwNnAlgDUFBCD05MY2z9Zly9FRMJYt4POOdcM3Am8AWwFnnHObTazu81sNoCZnWlmxcC1wMNmttm/+zigwMzWA28CP3HOhUzQgW+UlIKPjlBVq1FSREQCyTnHI28XsuPQ0R79nJhAvIlzbgmwpN2yH7R5vhrfJc32+70PTAxEDT3lgnE5PPjmbt7acZg5k7vSx0ZERD5N8ZE67lmylaT4aEbnJPfY54RCZ5SQNnlQGpl94/ib2ulERAJqY0kVABNzU3v0cxR0JxEVZZw/Npu3th+mSaOkiIgEzKaSKmKijDH9e+5sDhR0nXLBuByO1jezek+F16WIiESMjSVVjM5JJj4mukc/R0HXCeeMyiQuJkqXL0VEAsQ5x+b91T1+2RIUdJ2SFBfDWSMyWLbtEM45r8sREQl7+6vqqahpZMIgBV3IuHBcDh+V17Jbo6SIiHTbxmJfR5QJA1N6/LMUdJ10wbhsAJZu0eVLEZHu2ry/iugoY9wABV3IGJCayGkDU1imUVJERLptY0kVo7L7khDbsx1RQEHXJReMy2HN3iNU1DR6XYqISNhyzrGppIoJQeiIAgq6LrloXA6tDt7UHHUiIqfsUHUDZccag9LjEhR0XTIhN4WclHgN8iwi0g3HR0SZkNvz7XOgoOsSM+Pi8f15c/thjtZrkGcRkVOxsaSKKIPxA3RGF5KumppLfVMrr2066HUpIiJhaXNJFSOz+5IY1/MdUUBB12VTBqcxPLMPL6wp9roUEZGwtLGkigkDg3M2Bwq6LjMzrp6Sy8rCCoqP1HpdjohIWDlcXc/how1B63EJCrpTcuUU37x0L64t8bgSEZHwsmm/f2qeIAz9dZyC7hQMTk9i+vB0nl9TorEvRUS6YGNxNWYwPggjohwXkKAzs0vNbLuZ7TKz73Sw/lwzW2NmzWZ2Tbt188xsp/8xLxD1BMMXpg5iT1kNa/dVel2KiEjY2FhSRV5mH/rExwTtM7sddGYWDTwIXAaMB64zs/HtNtsL3Az8qd2+6cBdwHRgGnCXmfXrbk3BcNnE/iTERqlTiohIF2zeXxW0G8WPC8QZ3TRgl3Ou0DnXCDwFzGm7gXOuyDm3AWg/RfclwFLnXIVz7giwFLg0ADX1uOSEWC45rT8vrz9AQ3OL1+WIiIS8smMNHKiqD2pHFAhM0OUC+9q8LvYvC+i+ZjbfzArMrKC0tPSUCg20q6cOoqquSUOCiYh0wj9GRAm/oAsK59xC51y+cy4/KyvL63IAOHtEBtnJ8Ty/Rr0vRUROZrM/6E4Lwhx0bQUi6EqAwW1eD/Iv6+l9PRcTHcWVU3J5c9thyo81eF2OiEhI21hSxfDMPiQnxAb1cwMRdKuBUWY23MzigLnA4k7u+wZwsZn183dCudi/LGxcPTWX5lbHy+v3e12KiEhI21RSHfTLlhCAoHPONQN34guorcAzzrnNZna3mc0GMLMzzawYuBZ42Mw2+/etAH6ELyxXA3f7l4WNsf1TGD8ghRd087iIyAlV1DRSUlnHhCBftgQIyI0MzrklwJJ2y37Q5vlqfJclO9r3MeCxQNThlaun5vLjV7ey89BRRuUke12OiEjI2eRvnwv2rQUQRp1RQtmcyblER5nO6kRETuB4j8vTFHThKSs5nlmjs3hxbQktrRoSTESkvc37qxiSnkRqYnA7ooCCLmCunprLgap6Vuwu97oUEZGQs7Ek+COiHKegC5ALx+WQnBCjIcFERNqpqm1iX0WdJz0uQUEXMAmx0VwxaQCvbTpIdX2T1+WIiISM41PzTMgNfo9LUNAF1PXThlLX1MJTH+z1uhQRkZDx8dBfQZxVvC0FXQBNHJTKzLwMHnu3iMbm9uNXi4j0TptKqhjUL5F+feI8+XwFXYDNn5XHwep6FmukFBERwBd0Xp3NgYIu4M4bncWYnGQeebtQs4+LSK9XXd9EUXktEwcp6CKGmfHVc/PYfugob+0IjemERES8ssmjqXnaUtD1gNmnD6R/SgILlxd6XYqIiKfW7q0EYJKCLrLExURxy9nDWFFYzsbiKq/LERHxzOqiCkZl9/WsIwoo6HrMddOH0Dc+hoff3u11KSIinmhpdXxYdIQzh6d7WoeCroekJMRy/fQhLNl4gH0VtV6XIyISdNsOVnO0oZlpwxR0EeuWs4cRZcbv3t3jdSkiIkFXUHQEgPxh/TytQ0HXgwakJjJ78kCeXr2PIzWNXpcjIhJUHxRVMDA1gUH9kjytIyBBZ2aXmtl2M9tlZt/pYH28mT3tX7/KzIb5lw8zszozW+d/PBSIekLJ/HPzqGtq4cmVH3ldiohI0DjnKCiqIN/jy5YQgKAzs2jgQeAyYDxwnZmNb7fZrcAR59xI4AHgvjbrdjvnJvsft3e3nlAztn8Ks0ZnsWhFEfVNLV6XIyISFPsq6jhU3eB5RxQIzBndNGCXc67QOdcIPAXMabfNHGCR//lzwAVmZgH47LCw4Nw8yo418sIazUAuIr3D6qIKAM70uH0OAhN0ucC+Nq+L/cs63MY51wxUARn+dcPNbK2ZLTezcwJQT8iZOSKDCbkpPPpOoWYgF5FeYXVRBamJsYzOTva6FM87oxwAhjjnpgDfBP5kZh1OWGRm882swMwKSkvDa2gtM2P+uSMoLKvh1Y0HvC5HRKTHrS6qIH9oP6KivL94F4igKwEGt3k9yL+sw23MLAZIBcqdcw3OuXIA59yHwG5gdEcf4pxb6JzLd87lZ2VlBaDs4Lp84gDGD0jhJ0u2UteotjoRiVzlxxrYXVoTEh1RIDBBtxoYZWbDzSwOmAssbrfNYmCe//k1wN+dc87MsvydWTCzPGAUEJEDREZHGT+cfRr7q+p5aLlGSxGRyLXaf//ctOHet89BAILO3+Z2J/AGsBV4xjm32czuNrPZ/s1+B2SY2S58lyiP34JwLrDBzNbh66Ryu3Ouors1happw9O5YtIAHlq+m+IjGi1FRCJTQVEFcTFRns5Y0JaF45xp+fn5rqCgwOsyTsn+yjrOv/8tzh+bzW9uOMPrckREAm7Or98lPjaaZxbMDNpnmtmHzrn8jtZ53Rml1xmYlsgds0ayZONB3t9d5nU5IiIBVdvYzKb91SFxW8FxCjoPLJiVR25aIne/vIXmllavyxERCZi1eytpaXWcGSIdUUBB54mE2Gj++/JxbDt4lD9/sNfrckREAmZ1UQVRBmcM1Rldr3fphP7MzMvg/qU7NOCziESM1UUVjO2fQnJCrNelfExB5xEz467Z46mua+LnS3d4XY6ISLc1tbSydm8l00JgfMu2FHQeGts/hRtnDOWPqz5i64Fqr8sREemWLfurqW1s8Xz+ufYUdB775kWjSUmM5X9e3kw43uohInLcPwZy1hmdtJGWFMf/uWg0KwsreGWDxsEUkfC1uqiCIelJ5KQkeF3KP1HQhYDrpg1hYm4q//WXjeyr0IgpIhJ+fBOtHgm5szlQ0IWEmOgoHrx+KgB3/PFDTdAqImGnsKyG8prGkBnfsi0FXYgYkpHE/V+czKaSau5+ZYvX5YiIdMnqPb72uVCZsaAtBV0IuWh8DrfPGsGfVu3lhTXFXpcjItJpq4uOkNEnjrzMPl6X8gkKuhDzHxePZvrwdP7rL5vYfvCo1+WIiHTK6qIK8of1w8z7iVbbU9CFmJjoKH51/RT6JsRwxx8/5FhDs9cliYh8qkPV9eytqA3JjiigoAtJ2ckJ/Oq6KRSV1fDt5zfo/joRCWmhev/ccQq6EDUjL4P/vGQsr244wKL3i7wuR0TkhFYVVpAUF834gSlel9IhBV0IW3BuHheOy+GeJVtZs/eI1+WIiHyCc47lO0o5a0QGsdGhGSkBqcrMLjWz7Wa2y8y+08H6eDN72r9+lZkNa7Puu/7l283skkDUEymiooz7rz2d/qkJzH+igC37NR6miISWwrIa9lbUct6YbK9LOaFuB52ZRQMPApcB44HrzGx8u81uBY4450YCDwD3+fcdD8wFTgMuBX7jfz/xS02K5fFbphEbHcXchSt0ZiciIeXNbYcBOG9MlseVnFhMAN5jGrDLOVcIYGZPAXOAtnc9zwF+6H/+HPBr8/VBnQM85ZxrAPaY2S7/+60IQF0RY0RWX55ZMJMbf7eKGx9dxaPz8jlrRKbXZUW0llbH/so6SirrqK5rorq+2f+1iaP+5zWNzcRERREfE0VCbDTxMVHEx0YRH+N7np0Sz5D0PgxJTyKzb1xIdrsW6a7lO0oZld2XQf2SvC7lhAIRdLnAvjavi4HpJ9rGOddsZlVAhn/5ynb75nb0IWY2H5gPMGTIkACUHV4GpyfxrD/sbv79ah66cSrnj83xuqywd7Cqnh2HjlJUXkNRWa3va3kN+ypqaWrpuLdrn7hoUhJj6RMfQ0uro76phYbmVhr8X5tbP7lfUlw0Q9KTGJyexJD0JMbkJDNteDpDM5IUgBK2ahqaWVVYwc1nD/O6lE8ViKALCufcQmAhQH5+fq/sb5+dksBT82cy77EPmP/Eh/zv3MlcMWmg12WFjaraJjaUVLJ+XyXri6tYv6+Sw0cbPl6fGBvNsMw+jMlJ5pLT+jMsI4lB/ZJITYwlJSGWlMQY+sbHEHOSBvfmllbqm1s5WFXPvopaPiqvYW9FHXv9z9/ZWUp9UysA2cnxTBuezvTh6UwbnsGo7L5ERSn4JDy8v7ucxpZWzhsdupctITBBVwIMbvN6kH9ZR9sUm1kMkAqUd3JfaSO9Txx/+up0bn28gG/8eS21DS188czBJ9+xF9pXUcuK3eWsKCxn7d4jFJX/Y2aIvKw+nD0yk0mDUhk3IIXhmX3ITo4PyNlVTHQUfaOjGJndl5HZfT+x3jnH7tJjrNpTwQd7KljVZoqmtKRYZuZlcMWkgVwwLpuEWDVZS+h6c/th+sRFh+T4lm0FIuhWA6PMbDi+kJoLXN9um8XAPHxtb9cAf3fOOTNbDPzJzH4ODARGAR8EoKaIlpwQy6KvTGPBkx/yrec3UF3fxK2fGd7rL4Edqq5nxe5y3t9dxvu7yyk+UgdAZt84zhjaj2vzBzN5cBoTclNJTYz1rE4zY2R2MiOzk7lh+lCccxQfqfMHXzlvbS/ltU0H6RsfwyWn9WfO5IGcNSLjpGeSIsHknGP59lI+MyqTuJjQ/tnsdtD529zuBN4AooHHnHObzexuoMA5txj4HfAHf2eTCnxhiH+7Z/B1XGkGvuac0xw1nZAYF80jXz6Df/3zOn786lbe21XGj6+aSG5aotelBU19UwsrC8tZvqOU5TtKKSytASA1MZYZeel89Zw8Zo7wXQ4M5T8CzIzB/va7a84YREurY1VhOS+uK+G1TQd5fk0xmX3juWLSAOZMHsjkwWkh/f1I77Dz8DFKKuv4+vkjvS7lpCwch5fKz893BQUFXpcRElpaHYveL+Jnb2wnyuBbl47lphlDI7KdxzlHUXktb20/zFvbS1lZWE5DcyvxMVFMz8vgnJGZzByRwbgBKURHyPdf39TCW9tLeWldCcu2HaaxuZXTB6dxx6w8LhrfP2K+Twk/Dy/fzb2vbWPFd89nQKr3f2Cb2YfOufwO1ynoIsO+ilq+95eNvLOzjKlD0rjvC5MYlZPsdVndVtfYworCMt7a7jtr+8jfzpaX2YdZY7KYNTqLGXkZvaItq7q+iZfW7efRdwr5qLyWvMw+fPXcPK6aktsrvn8JLXMXrqCytonX/+1cr0sBFHS9hnOOv6wt4e5XtlDT0MzXPjuSfzlvZMhfP2/LOUdhWQ1vbS/lre2HWbWngsbmVhJjozlrRAbnjcli1uhshmSE7j07Pa2l1fH6poM8tHw3G0uqyEqO55azh3HD9KGetj1K73G0vokpdy/lq+fm8e1Lx3pdDqCg63XKjjVw98tbWLx+P6Oy+3L7rBF8buIAEuNC86/+IzWNrCj0dSJ5e0cZeyt8Z20jsvpw3phszhuTxZnD0nXW0o5zjhW7y3no7ULe3lFK3/gYbv3McBbMyiMpLmzuHJIw9PqmA9z+5Bqenj+D6XkZXpcDKOh6rb9vO8SPXtnKnrIakuNjmD15IF86czATc1M97cxQ09DMB0UVvL+rjPd2lbP1YDXO+W7EnpHnO2s7b0w2g9N771lbV23ZX82Db+7i1Y0H6J+SwLcvG8Oc03Mjsq1WvPft5zawZNMB1nz/opAZyFlB14s551i1p4JnVu9jyaYD1De1MrZ/Ml86czBXTcklLSmuRz+/tdV3KXJDcSUbiqtYX1zJxuIqmlsdcdFRTBmSxtkjMzl7ZAaTBqWFzH+acFVQVMGPXtnC+uIqTh+Uyg8+P54zhob2PU4SXpxzzLh3GflD03nwhqlel/MxBZ0Avs4Mi9ft55mCfWworiIuOorJg9MY0z+ZsQOSGds/mdE5ySQndL2dp6XVUXasgZLKOvZV1LJ5fzXr91WyeX/1x7OkJ8ZGc9rAFPKHpXP2yAzyh6aH7OXUcNba6nhxXQn3vb6NQ9UNfP70gXz70jEhPRahhI8t+6v53C/f4afXTOKL+aEzWIWCTj5h64Fqnv+wmHX7Ktl+8ChH/WEEMKhfImP7JzMgNZHoKCPKjOgoiI6K8n01o7HFcbCqjv2V9eyvquNgVf0/jfEYFxPFuAEpnD4olYm5qUwalMaIrD666TmIahubeWh5IQ8v3w3AHeeNCLvOSRJ6HnxzFz97YzsffO8CslMSvC7nYwo6+VTOOUoq69h+8Cjbjj8OVFN2rIGWVker852xtbQ6Wpzva2y00T81gYGpiQxMS2RAagID0xIZmOb7mpfZV79QQ0RJZR33LtnKKxsOMG5ACvdfe3rIzgQtoe/ah96nrqmFV75+jtel/JNPCzp1zRLMjEH9fAMYXzCuczMiOOc0OkeYyE1L5NfXT2X26Qf53l82MfvX7/L180fxL58doTZR6ZKq2ibW7K3kjlkjvC6lS/RTLqdEIRd+Lj6tP0v//VwunzSAB/62gysffI9tBzVrvXTeO7tKaWl1fHZsaM9W0J6CTqQX6dcnjl/MncJDN57Boep6Pv+rd/nVsp00tbR6XZqEgTe3lZKWFMvkwf28LqVLFHQivdClE/rz13+fxSWn9ef+pTv4wm/fZ2+baYxE2mttdSzfcZhzRmWF3RirCjqRXiq9Txy/vn4qD14/laKyGi7/1Tu8sfmg12VJiNq8v5qyY418dkx4XbYEBZ1Ir3f5pAG8+o1zGJbRhwV/+JB7Xt2iS5nyCW9uP4wZnBvis4l3REEnIgxOT+K5O2Zy04yhPPLOHq5buJIDVXVelyUhZNm2w0zKTSWzb7zXpXSZgk5EAIiPieZHV07gF3Mns+VANZf/8l3e2VnqdVkSAj4qr2H9vkounTDA61JOSbeCzszSzWypme30f+2wK46ZzfNvs9PM5rVZ/paZbTezdf5HdnfqEZHumzM5l8V3fobMvnF8+bEPeGDpDlpaw29gCQmcl9btB2DO5IEeV3JquntG9x1gmXNuFLDM//qfmFk6cBcwHZgG3NUuEG9wzk32Pw53sx4RCYCR2X158Wtnc9WUXH6xbCcL/vAhNW2GiZPewznHi2tLmJGXzsA072cSPxXdDbo5wCL/80XAlR1scwmw1DlX4Zw7AiwFLu3m54pID0uKi+H+a0/nf2afxt+3HeKah1ZQUql2u95mQ3EVhWU1XDk51+tSTll3gy7HOXfA//wg0NH4UbnAvjavi/3Ljvu9/7Ll903DbYiEFDNj3lnD+P0t0yiuqGXOr99jzd4jXpclQfTiuhLioqO4bGJ4ts9BJ4LOzP5mZps6eMxpu53zjQ7d1Qv5NzjnJgLn+B83fUod882swMwKSkvVQC4STLNGZ/HCv5xFUlw0cxeu5KV1JV6XJEHQ3NLKy+v3c8G4bFITuz59V6g4adA55y50zk3o4PEScMjMBgD4v3bUxlYCtJ20aJB/Gc6541+PAn/C14Z3ojoWOufynXP5WVnhdx+HSLgblZPMi187m8mD0vjXp9bx86U7aFUnlYj27q4yyo41MieML1tC9y9dLgaO96KcB7zUwTZvABebWT9/J5SLgTfMLMbMMgHMLBa4AtjUzXpEpAel94njD7dN45ozBvHLZTv5+lNrqWts8bos6SEvrdtPSkJM2A3i3F53g+4nwEVmthO40P8aM8s3s0cBnHMVwI+A1f7H3f5l8fgCbwOwDt9Z3iPdrEdEelh8TDQ/u2YS371sLEs2HuCGR1dypKbR67IkwGoamnl900EunzSQ+Jhor8vplm7NR+ecKwcu6GB5AXBbm9ePAY+126YGOKM7ny8i3jAzFswawZD0JP716XVc89D7PHHrdHLDtPu5fNLSLYeoa2rhyjC9d64tjYwiIqfssokDeOIr0zh8tIGrf6P57SLJX9aWkJuWyJnD0r0updsUdCLSLTPyMnj29pkAXPvQClYVlntckXRX6dEG3t1VxpzJA4kKsyl5OqKgE5FuG9s/hefvOIus5HhueuwDXt904OQ7Sch6ZcN+WlodV00J796WxynoRCQgBvVL4vnbz+K0gSnc8cc1/GHlR16XJKfoxbUljB+QwqicZK9LCQgFnYgETL8+cfzpthmcPyab77+4iZ8v3YFvLAkJF4Wlx1hfXBUxZ3OgoBORAEuMi+bhm87gWv+9dv/z8hbdWB5GXly3HzOYHQG9LY/r1u0FIiIdiYmO4r4vTKJvQgy/f6+ImoZmfvKFSURHQMeGSHZ8poKzRmSQk5LgdTkBo6ATkR4RFWX84IrxJCfE8stlO6ltbOGBL00mLkYXkkLV2n2V7K2o5evnj/S6lIBS0IlIjzEzvnnRaJLjY7hnyVZqGpv57Q1nkBgX3iNtRKoX15YQHxPFpRP6e11KQOlPKxHpcV89N4//e9VElu8oZd7vP+BofZPXJUk7jc2tvLLhABeOzyE5IXxnKuiIgk5EguL66UP43y9NZs1HR7jh0VUaHzPEvLSuhIqaRr6YP/jkG4cZBZ2IBM2cybk8fNMZbDt4lC8tXMHho/VelyRAa6vjkXcKGds/mXNHZXpdTsAp6EQkqC4Yl8PjN59J8ZE65j68kgNVdV6X1Ou9teMwOw4dY8GsPMwir2esgk5Egu6skZkfDwb9xYdXsK+i1uuSerWHlxcyMDWBKyZFzr1zbSnoRMQT+cPS+eNt06mua+aLD6+gsPSY1yX1Suv2VbJqTwVf+cxwYqMjMxIi87sSkbBw+uA0/vzVGTQ2t/KlhSvZceio1yX1Ogvf3k1yQgxzpw3xupQeo6ATEU+NH5jCU/NnYMDchSvZvL/K65J6jaKyGl7bdJCbZgylb3zk3lbdraAzs3QzW2pmO/1f+51gu9fNrNLMXmm3fLiZrTKzXWb2tJnFdaceEQlPo3KSeWbBTBJiorhu4UrW7av0uqRe4dF3C4mNiuLms4Z5XUqP6u4Z3XeAZc65UcAy/+uO/Ay4qYPl9wEPOOdGAkeAW7tZj4iEqWGZfXh6wUzSkuK48dFVfLCnwuuSIlrZsQaeLSjm6qm5ZEfQuJYd6W7QzQEW+Z8vAq7saCPn3DLgny6+m68P6/nAcyfbX0R6h8HpSTyzYCbZKfHMe+wD3t1Z5nVJEeuJFR/R0NzKbefkeV1Kj+tu0OU4545PJXwQyOnCvhlApXOu2f+6GDjhBEhmNt/MCsysoLS09NSqFZGQ1z81gafnz2RoRhJfWbSav2875HVJEae2sZknVhRx0fgcRmb39bqcHnfSoDOzv5nZpg4ec9pu53yzK/bYpFPOuYXOuXznXH5WVlZPfYyIhICs5Hj+/NUZjO2fzPwnPmTJxgMn30k67dmCYiprm1hwbuSfzUEngs45d6FzbkIHj5eAQ2Y2AMD/9XAXPrscSDOz4119BgElXf0GRCQy9esTx5O3Tef0wWnc+ac1/GVtsdclRYTmllYeeaeQM4b2I39YutflBEV3L10uBub5n88DXursjv4zwDeBa05lfxGJfCkJsTzxlWnMyMvgm8+s588f7PW6pLD32qaDFB+pY34vOZuD7gfdT4CLzGwncKH/NWaWb2aPHt/IzN4BngUuMLNiM7vEv+rbwDfNbBe+NrvfdbMeEYkwfeJjeOzmM5k1OovvvrCRx97d43VJYcs5x8Nv7yYvsw8XjetKl4rw1q07BJ1z5cAFHSwvAG5r8/qcE+xfCEzrTg0iEvkSYqN5+KYz+Maf13L3K1uobWzma58dGZEDEPek5TtK2VRSzb1XTyQqqvccO42MIiJhIT4mmgevn8pVU3L5f3/dwT2vbsXXAiKdUd/Uwg8Xb2Z4Zh+unnrCDu4RKXLHfBGRiBMTHcX9155OamIsj767h6q6Ju69eiIxEToYcSA9tHw3ReW1/OHWacTHRHtdTlAp6EQkrERFGXd9fjypibH8YtlOquub+MXcKSTE9q5f3l1RVFbDb97azRWTBnDOqN53e5b+DBKRsGNm/PtFo7nr8+N5Y/MhvvL4ao41NJ98x17IOcf3X9pEXHQU379ivNfleEJBJyJh65azh/PzL57Oqj0V3PDISo7UNHpdUshZsvEg7+ws4/9cPJqcCB/T8kQUdCIS1q6eOoiHbzyDrQePcu3DKzhQVed1SSHjaH0Td7+ymdMGpnDTjKFel+MZBZ2IhL0Lx+ew6JZpHKyq5wu/eZ9tB6u9LikkPLB0J4ePNnDPVb27w07v/c5FJKLMHJHB0wtm0OIc1/x2BW/v6N2Dv2/eX8Xj7+/h+mlDmDw4zetyPKWgE5GIcdrAVF782tkM6pfILY+v7rVDhrW2Ov77xU2k94njW5eM9boczynoRCSiDEhN5NnbZ/KZkZl894WN3Pf6Nlpbe9eN5U+t3sfavZV873PjSE2K9boczynoRCTiJCfE8rt5+Vw/fQi/fWs333hqLfVNLV6XFRRlxxq47/VtTB+ezlVTetcIKCeiG8ZFJCLFREdxz5UTGJqexL2vbeNAVT2PfDmf9D5xXpfWY5pbWvnWcxuoaWjmx1dO0FigfjqjE5GIZWYsmDWCB6+fysaSKq588D0276/yuqwe4Zzjhy9v5u/bDnPX58czKifZ65JChoJORCLe5ZMG8NT8GTQ0t3D1b97nmdX7vC4p4Ba+XciTK/ey4Nw8bpo5zOtyQoqCTkR6halD+vHqN84hf1g/vvX8Bv7z2fXUNUZGu90rG/Zz72vbuHzSAL59qXpZtqegE5FeI7NvPE98ZTrfOH8kz35YzFW/eY89ZTVel9UtBUUVfPOZ9eQP7cf9157eq+aZ66xuBZ2ZpZvZUjPb6f/a7wTbvW5mlWb2Srvlj5vZHjNb539M7k49IiInEx1lfPPiMfz+ljM5WF3P7F+9y+ubDnhd1ikpLD3GbU8UMCgtkUe+nK8ZHE6gu2d03wGWOedGAcv8rzvyM+CmE6z7T+fcZP9jXTfrERHplM+OyebVb5xDXnZfbn9yDT9+ZQsNzeFzKbPsWAM3/3410Wb8/pYz6RfBvUm7q7tBNwdY5H++CLiyo42cc8uAo938LBGRgMpNS+TZBTO5+axhPPruHj73i3dYVVjudVknVdfYwm2LCjh8tJ5H5+UzNKOP1yWFtO4GXY5z7vg5/0Eg5xTe4x4z22BmD5hZfDfrERHpkriYKH44+zQWfWUajS2tfGnhSr793AYqa0Nzyp+KmkYWPPkh64sr+cXcKUwZ0mGLkbRx0qAzs7+Z2aYOHnPabuecc0BXx9n5LjAWOBNIB779KXXMN7MCMysoLe3dg7WKSODNGp3FX/9tFgtm5fHcmmIu/PlyXlpXgu9XW2h4f1cZl/3ibVbuLufeqyZyyWn9vS4pLFh3/hHNbDtwnnPugJkNAN5yzo05wbbnAf/hnLviVNa3lZ+f7woKCk61bBGRT7V5fxXfe2Ej64urOHd0FvdcOYHB6Ume1dPU0srPl+7goeW7GZ7Zh1/OncKE3FTP6glFZvahcy6/o3XdvXS5GJjnfz4PeKmLhQ3wfzV87XubulmPiEi3nTYwlRf+5Wzu+vx4Piyq4KIHlnPva1s5WFUf9FqKymq45rfv89u3djP3zCG88vXPKOS6qLtndBnAM8AQ4CPgi865CjPLB253zt3m3+4dfJco+wLlwK3OuTfM7O9AFmDAOv8+x072uTqjE5Fg2V9Zx72vbePVDfuJjjJmn57LV88dztj+KT36uc45nl9Twl0vbSImOoqfXD2RyyYO6NHPDGefdkbXraDzioJORIJtb3ktj723h6dX76OuqYVzR2cx/5w8zh6ZEdDBk1tbHauLKli0ooglGw8yfXg6D3xpMgPTEgP2GZFIQSciEiCVtY08ufIjHn//I8qONTBuQApfzB/EmcPSGTcghehTGJnEOcfm/dUsXr+fl9fv50BVPQmxUdz52ZHccd7IU3rP3kZBJyISYPVNLby0roRH39nDzsO+Fpe+8TFMGZLGmcPSyR/Wj8mD00iK++RsaC2tjvqmFg5U1fPqhgMsXl/C7tIaYqKMWaOzmD15IBeOy6FPvGZS6ywFnYhIDyqprKOgqILVRRUUFB1h+6GjOAcxUcagfok0tfiCraG5lfqmFprbzXg+fXg6sycP5HMTBmiEk1P0aUGnPxdERLopNy2R3Mm5zJnsm9G7qraJNXuPsLqogn1H6oiPiSI+JoqE2Oh/+pqcEMt5Y7LU/tbDFHQiIgGWmhTLZ8dm89mx2V6XImiaHhERiXAKOhERiWgKOhERiWgKOhERiWgKOhERiWgKOhERiWgKOhERiWgKOhERiWhhOQSYmZXimxaoJ2UCZT38GeFGx6RjOi6fpGPSMR2XjgXiuAx1zmV1tCIsgy4YzKzgROOm9VY6Jh3TcfkkHZOO6bh0rKePiy5diohIRFPQiYhIRFPQndhCrwsIQTomHdNx+SQdk47puHSsR4+L2uhERCSi6YxOREQimoJOREQimoLOz8yuNbPNZtZqZifs5mpml5rZdjPbZWbfCWaNwWZm6Wa21Mx2+r/2O8F2LWa2zv9YHOw6g+Vk//ZmFm9mT/vXrzKzYR6UGVSdOCY3m1lpm5+P27yoM5jM7DEzO2xmm06w3szsl/5jtsHMpga7Ri904ricZ2ZVbX5WfhCoz1bQ/cMm4Grg7RNtYGbRwIPAZcB44DozGx+c8jzxHWCZc24UsMz/uiN1zrnJ/sfs4JUXPJ38t78VOOKcGwk8ANwX3CqDqwv/H55u8/PxaFCL9MbjwKWfsv4yYJT/MR/4bRBqCgWP8+nHBeCdNj8rdwfqgxV0fs65rc657SfZbBqwyzlX6JxrBJ4C5vR8dZ6ZAyzyP18EXOldKZ7rzL992+P1HHCBmVkQawy23vb/oVOcc28DFZ+yyRzgCeezEkgzswHBqc47nTguPUZB1zW5wL42r4v9yyJVjnPugP/5QSDnBNslmFmBma00syuDU1rQdebf/uNtnHPNQBWQEZTqvNHZ/w9f8F+ie87MBgentJDW236PdMVMM1tvZq+Z2WmBetOYQL1RODCzvwH9O1j1X865l4JdTyj4tGPS9oVzzpnZie5FGeqcKzGzPODvZrbRObc70LVKWHoZ+LNzrsHMFuA74z3f45okNK3B97vkmJl9DngR3+XdbutVQeecu7Cbb1ECtP2LdJB/Wdj6tGNiZofMbIBz7oD/0srhE7xHif9roZm9BUwBIi3oOvNvf3ybYjOLAVKB8uCU54mTHhPnXNvv/1Hgp0GoK9RF3O+RQHDOVbd5vsTMfmNmmc65bg+CrUuXXbMaGGVmw80sDpgLRGwvQ3zf2zz/83nAJ856zayfmcX7n2cCZwNbglZh8HTm377t8boG+LuL7BEZTnpM2rU9zQa2BrG+ULUY+LK/9+UMoKpNE0GvZWb9j7dpm9k0fPkUkD8Ue9UZ3acxs6uAXwFZwKtmts45d4mZDQQedc59zjnXbGZ3Am8A0cBjzrnNHpbd034CPGNmt+KbFumLAP7bL253zt0GjAMeNrNWfD+YP3HORVzQnejf3szuBgqcc4uB3wF/MLNd+Brd53pXcc/r5DH5hpnNBprxHZObPSs4SMzsz8B5QKaZFQN3AbEAzrmHgCXA54BdQC1wizeVBlcnjss1wB1m1gzUAXMD9YeihgATEZGIpkuXIiIS0RR0IiIS0RR0IiIS0RR0IiIS0RR0IiIS0RR0IiIS0RR0IiIS0RR0ImHIzO4zs381s03+AbXHeV2TSKhS0ImEGTM7C7gYWIdvjMQfAr/wsCSRkKagEwk/0/DNCmBAE/A6cIanFYmEMAWdSPjpaNy+lqBXIRImFHQi4ecd4HL8A+ICX/AvE5EOaPYCkTDjnFtjZs8BjwDp+Oa9u9HbqkRCl2YvEAlTZnYe8B/OuSs8LkUkpOnSpYiIRDSd0YmISETTGZ2IiEQ0BZ2IiEQ0BZ2IiEQ0BZ2IiEQ0BZ2IiES0/w9pkx2rB4DCPAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "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": [ "$\\displaystyle \\bar{\\epsilon} \\left(δ \\cos{\\left(j \\left(- θ_{0} + \\operatorname{atan}_{2}{\\left({\\partial_{1} {φ}_{(0,0)}},{\\partial_{0} {φ}_{(0,0)}} \\right)}\\right) \\right)} + 1\\right)$" ], "text/plain": [ "\\bar{\\epsilon}⋅(δ⋅cos(j⋅(-θ₀ + atan2(D(φ[0,0]), D(φ[0,0])))) + 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": [ "$\\displaystyle {φ}_{(0,0)}^{3} - \\frac{{φ}_{(0,0)}^{2} α \\operatorname{atan}{\\left({T}_{(0,0)} γ - T_{eq} γ \\right)}}{\\pi} - \\frac{3 {φ}_{(0,0)}^{2}}{2} + \\frac{{φ}_{(0,0)} α \\operatorname{atan}{\\left({T}_{(0,0)} γ - T_{eq} γ \\right)}}{\\pi} + \\frac{{φ}_{(0,0)}}{2} - \\bar{\\epsilon}^{2} δ^{2} \\cos^{2}{\\left(j θ_{0} - j \\operatorname{atan}_{2}{\\left({\\partial_{1} {φ}_{(0,0)}},{\\partial_{0} {φ}_{(0,0)}} \\right)} \\right)} {\\partial_{0} {\\partial_{0} {φ}_{(0,0)}}} - \\bar{\\epsilon}^{2} δ^{2} \\cos^{2}{\\left(j θ_{0} - j \\operatorname{atan}_{2}{\\left({\\partial_{1} {φ}_{(0,0)}},{\\partial_{0} {φ}_{(0,0)}} \\right)} \\right)} {\\partial_{1} {\\partial_{1} {φ}_{(0,0)}}} - 2 \\bar{\\epsilon}^{2} δ \\cos{\\left(j θ_{0} - j \\operatorname{atan}_{2}{\\left({\\partial_{1} {φ}_{(0,0)}},{\\partial_{0} {φ}_{(0,0)}} \\right)} \\right)} {\\partial_{0} {\\partial_{0} {φ}_{(0,0)}}} - 2 \\bar{\\epsilon}^{2} δ \\cos{\\left(j θ_{0} - j \\operatorname{atan}_{2}{\\left({\\partial_{1} {φ}_{(0,0)}},{\\partial_{0} {φ}_{(0,0)}} \\right)} \\right)} {\\partial_{1} {\\partial_{1} {φ}_{(0,0)}}} - \\bar{\\epsilon}^{2} {\\partial_{0} {\\partial_{0} {φ}_{(0,0)}}} - \\bar{\\epsilon}^{2} {\\partial_{1} {\\partial_{1} {φ}_{(0,0)}}}$" ], "text/plain": [ " 2 2 \n", " 3 φ_C ⋅α⋅atan(T_C⋅γ - T_eq⋅γ) 3⋅φ_C φ_C⋅α⋅atan(T_C⋅γ - T_eq⋅γ) φ_C 2 2 2 \n", "φ_C - ─────────────────────────── - ────── + ────────────────────────── + ─── - \\bar{\\epsilon} ⋅δ ⋅cos (j⋅θ₀ - j⋅atan2(D\n", " π 2 π 2 \n", "\n", " \n", " 2 2 2 \n", "(φ[0,0]), D(φ[0,0])))⋅D(D(φ[0,0])) - \\bar{\\epsilon} ⋅δ ⋅cos (j⋅θ₀ - j⋅atan2(D(φ[0,0]), D(φ[0,0])))⋅D(D(φ[0,0])) - 2⋅\\bar{\n", " \n", "\n", " \n", " 2 2 \n", "\\epsilon} ⋅δ⋅cos(j⋅θ₀ - j⋅atan2(D(φ[0,0]), D(φ[0,0])))⋅D(D(φ[0,0])) - 2⋅\\bar{\\epsilon} ⋅δ⋅cos(j⋅θ₀ - j⋅atan2(D(φ[0,0]), D\n", " \n", "\n", " \n", " 2 2 \n", "(φ[0,0])))⋅D(D(φ[0,0])) - \\bar{\\epsilon} ⋅D(D(φ[0,0])) - \\bar{\\epsilon} ⋅D(D(φ[0,0]))\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": { "text/latex": [ "$\\displaystyle \\left\\{ \\pi : 3.14159265358979, \\ T_{eq} : 1.0, \\ \\bar{\\epsilon} : 0.01, \\ j : 6, \\ α : 0.9, \\ γ : 10, \\ δ : 0.02, \\ θ_{0} : 0.2, \\ κ : 1.8, \\ τ : 0.0003\\right\\}$" ], "text/plain": [ "{π: 3.14159265358979, T_eq: 1.0, \\bar{\\epsilon}: 0.01, j: 6, α: 0.9, γ: 10, δ: 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(φ_tmp, 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 'Target.CPU' or 'Target.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", " dh.swap('phi', 'phi_temp')\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": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
FUNC_PREFIX void kernel(double * RESTRICT const _data_T, double * RESTRICT const _data_phi, double * RESTRICT  _data_phi_temp, double * RESTRICT  _data_phidelta)\n",
       "{\n",
       "   #pragma omp parallel num_threads(4)\n",
       "   {\n",
       "      #pragma omp for schedule(static)\n",
       "      for (int64_t ctr_1 = 1; ctr_1 < 301; ctr_1 += 1)\n",
       "      {\n",
       "         double * RESTRICT _data_phi_10 = _data_phi + 302*ctr_1;\n",
       "         double * RESTRICT _data_T_10 = _data_T + 302*ctr_1;\n",
       "         double * RESTRICT _data_phi_11 = _data_phi + 302*ctr_1 + 302;\n",
       "         double * RESTRICT _data_phi_1m1 = _data_phi + 302*ctr_1 - 302;\n",
       "         double * RESTRICT  _data_phidelta_10 = _data_phidelta + 302*ctr_1;\n",
       "         double * RESTRICT  _data_phi_temp_10 = _data_phi_temp + 302*ctr_1;\n",
       "         for (int64_t ctr_0 = 1; ctr_0 < 301; ctr_0 += 1)\n",
       "         {\n",
       "            const double xi_0 = _data_phi_10[ctr_0]*_data_phi_10[ctr_0];\n",
       "            const double xi_1 = 954.92965855137209*atan(-10.0 + 10.0*_data_T_10[ctr_0]);\n",
       "            const double xi_2 = -2222.2222222222222*_data_phi_10[ctr_0];\n",
       "            const double xi_3 = xi_2 + 1111.1111111111111*_data_phi_11[ctr_0] + 1111.1111111111111*_data_phi_1m1[ctr_0];\n",
       "            const double xi_4 = cos(-1.2000000000000002 + 6.0*atan2(-16.666666666666668*_data_phi_1m1[ctr_0] + 16.666666666666668*_data_phi_11[ctr_0], -16.666666666666668*_data_phi_10[ctr_0 - 1] + 16.666666666666668*_data_phi_10[ctr_0 + 1]));\n",
       "            const double xi_5 = xi_4*0.013333333333333336;\n",
       "            const double xi_6 = xi_2 + 1111.1111111111111*_data_phi_10[ctr_0 + 1] + 1111.1111111111111*_data_phi_10[ctr_0 - 1];\n",
       "            const double xi_7 = 0.00013333333333333334*(xi_4*xi_4);\n",
       "            _data_phidelta_10[ctr_0] = xi_0*xi_1 + xi_0*5000.0 + xi_1*-1.0*_data_phi_10[ctr_0] + xi_3*xi_5 + xi_3*xi_7 + xi_5*xi_6 + xi_6*xi_7 - 3148.1481481481483*_data_phi_10[ctr_0] - 3333.3333333333335*(_data_phi_10[ctr_0]*_data_phi_10[ctr_0]*_data_phi_10[ctr_0]) + 370.37037037037038*_data_phi_10[ctr_0 + 1] + 370.37037037037038*_data_phi_10[ctr_0 - 1] + 370.37037037037038*_data_phi_11[ctr_0] + 370.37037037037038*_data_phi_1m1[ctr_0];\n",
       "            _data_phi_temp_10[ctr_0] = 1.0000000000000001e-5*_data_phidelta_10[ctr_0] + _data_phi_10[ctr_0];\n",
       "         }\n",
       "      }\n",
       "   }\n",
       "}\n",
       "
\n" ], "text/plain": [ "FUNC_PREFIX void kernel(double * RESTRICT const _data_T, double * RESTRICT const _data_phi, double * RESTRICT _data_phi_temp, double * RESTRICT _data_phidelta)\n", "{\n", " #pragma omp parallel num_threads(4)\n", " {\n", " #pragma omp for schedule(static)\n", " for (int64_t ctr_1 = 1; ctr_1 < 301; ctr_1 += 1)\n", " {\n", " double * RESTRICT _data_phi_10 = _data_phi + 302*ctr_1;\n", " double * RESTRICT _data_T_10 = _data_T + 302*ctr_1;\n", " double * RESTRICT _data_phi_11 = _data_phi + 302*ctr_1 + 302;\n", " double * RESTRICT _data_phi_1m1 = _data_phi + 302*ctr_1 - 302;\n", " double * RESTRICT _data_phidelta_10 = _data_phidelta + 302*ctr_1;\n", " double * RESTRICT _data_phi_temp_10 = _data_phi_temp + 302*ctr_1;\n", " for (int64_t ctr_0 = 1; ctr_0 < 301; ctr_0 += 1)\n", " {\n", " const double xi_0 = _data_phi_10[ctr_0]*_data_phi_10[ctr_0];\n", " const double xi_1 = 954.92965855137209*atan(-10.0 + 10.0*_data_T_10[ctr_0]);\n", " const double xi_2 = -2222.2222222222222*_data_phi_10[ctr_0];\n", " const double xi_3 = xi_2 + 1111.1111111111111*_data_phi_11[ctr_0] + 1111.1111111111111*_data_phi_1m1[ctr_0];\n", " const double xi_4 = cos(-1.2000000000000002 + 6.0*atan2(-16.666666666666668*_data_phi_1m1[ctr_0] + 16.666666666666668*_data_phi_11[ctr_0], -16.666666666666668*_data_phi_10[ctr_0 - 1] + 16.666666666666668*_data_phi_10[ctr_0 + 1]));\n", " const double xi_5 = xi_4*0.013333333333333336;\n", " const double xi_6 = xi_2 + 1111.1111111111111*_data_phi_10[ctr_0 + 1] + 1111.1111111111111*_data_phi_10[ctr_0 - 1];\n", " const double xi_7 = 0.00013333333333333334*(xi_4*xi_4);\n", " _data_phidelta_10[ctr_0] = xi_0*xi_1 + xi_0*5000.0 + xi_1*-1.0*_data_phi_10[ctr_0] + xi_3*xi_5 + xi_3*xi_7 + xi_5*xi_6 + xi_6*xi_7 - 3148.1481481481483*_data_phi_10[ctr_0] - 3333.3333333333335*(_data_phi_10[ctr_0]*_data_phi_10[ctr_0]*_data_phi_10[ctr_0]) + 370.37037037037038*_data_phi_10[ctr_0 + 1] + 370.37037037037038*_data_phi_10[ctr_0 - 1] + 370.37037037037038*_data_phi_11[ctr_0] + 370.37037037037038*_data_phi_1m1[ctr_0];\n", " _data_phi_temp_10[ctr_0] = 1.0000000000000001e-5*_data_phidelta_10[ctr_0] + _data_phi_10[ctr_0];\n", " }\n", " }\n", " }\n", "}" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ps.show_code(φ_kernel)" ] }, { "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", "phi_temp| ( 0, 0)| ( 0, 0)\n", "phidelta| ( 0, 0)| ( 0, 0)\n", "\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7kAAAF1CAYAAAAtEi0mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAB9t0lEQVR4nO29ffSuV1nf+b04geALQkIUs5LYxBKmE3WKkAl0rJQxvARWJbQCBqnETih1FVod7ayGYQSNugbsWCsz1DZASnAoL0Itp21sGt7GtqvEnCAiCY05RjHnNBLzQkQsYHKu+eO573P27/pd17X3/bzc9/N77u9nrd/6Pc9+vfa+7/Pd33vv5zw/UVUQQgghhBBCCCG7wKOmDoAQQgghhBBCCFkXfMglhBBCCCGEELIz8CGXEEIIIYQQQsjOwIdcQgghhBBCCCE7Ax9yCSGEEEIIIYTsDHzIJYQQQgghhBCyM/AhlxBCCCGEEELIzsCHXEIIIbNGRP6k+DkhIv+teP+KqeMjhBBCyDD4kEsIIWTWqOrX9z8A/gDA9xZp7546PkIIWQUReZqI3CUifyAiPzZ1PISMAR9ySYqI/E8i8lsi8iwR+SMRuVVELpo6LkIIIYQQ0sQdAC4CcDmA14vId4vI2SLy70TkRSJyu4jcLSJ/feI4CVkbfMglISLytQB+BcBbAbwCwLsBvBfAu0VEpoyNEEIIIYTUUdUvqeqXVfU3sfByzwfwTwD8EYDHA/gqgJcAuF5EzpkuUkLWBx9yScYlAATA2wB8HYAvAPi/AXwHgG+ZLixCCCGEENKCiFwlIjeLyPUAHgHwWADfC+BnAZwO4AuqejOA3wLwnOkiJWR98CGXZDwJwHFV1T5BVb8M4EEA3zxZVIQQQgghpIqIfAeAnwHwYgD/DMCrsDjBPQTgblP8HtDfkR2BD7kk4ziAc8uPJovI1wA4A8CxyaIihBBCCCEtPAfAv1bVewDcAuBLAN6DxUeU7afyzgH9HdkR+JBLMm7GQgxfg8XHlg8B+CkA/0lVj08ZGCGEEEIIqfLHWPg3APgJAL+uqn+AxXes/ASAxwCAiFwO4C8AuGGKIAlZN6dNHQDZXlT1zzrRuw7A07D4fxz/EcAPThoYIYQQQghp4d0A/pqIfBLA5wD8L136j2DxnSs/B+BrAHwjgL+uqg9OEiUha0aK/25JSIiI/L8AjqrqT04dCyGEEEIIWR0ReRWAv6Gqz546FkLWCT+uTAghhBBCCCFkZ+BDLiFko4jIdSJyr4h8JsgXEXmLiBwVkU+LyNPGjpEQQlZFRC4TkTs6LbvayX+WiHxSRB4WkZeYvCtF5M7u58oi/eki8ttdm2/h36gnhEzNQdE6flyZELJRRORZAP4EwLtU9dud/BcC+LsAXgjgGQB+UVWfMW6UhBCyPCJyCMDvAHguFt9OewuAl6vq7UWZ8wF8A4C/D+Cwqn6gSz8TwBEAFwNQALcCeLqqPigivwHg72HxRZA3AHiLqv7aWOMihJCSg6R1PMklhGwUVf11AA8kRS7H4gFYVfUTAJ4gImePEx0hhKyFS7D43oq7VPWrWHxz7eVlAVX9fVX9NIATpu7zAdykqg90X/pzE4DLOh38BlX9RPf36t+Fxd86JYSQqTgwWseHXELI1JyDvX+Q/liXRgghB4VVdCyqa/9mKbWREDI1B0brtuJPCJ111ll6/vnnTx0GIQeaW2+99T5V/cYhdZ7/P3+d3v/AI6v1++mv3Abgy0XStap67UqN7ijUOkJWh1q3/VDrCFkdat1qbMVD7vnnn48jR45MHQYhBxoR+dzQOvc98AhuvvHclfp99Nm/+2VVvXiFJo4DOK94f26XtnNQ6whZnS3VulV07DiAZ5u6H+/SzzXpB0IbqXWErA61bjX4cWVCZo3iET2x0s8aOAzgld23LD8TwEOqes86GiaEkAUb17pbAFwoIheIyGMAXIGFtrVwI4DnicgZInIGgOcBuLHTwT8WkWd23zT6SgAfWm78hJB5QK3r2YqTXELINCiAE9jsN6yLyHuw2Lk7S0SOAXgjgEcDgKr+Uyy+Re+FAI4C+FMAf3OjARFCZsemtU5VHxaR12Jh4g4BuE5VbxORawAcUdXDIvI/AvhVAGcA+F4R+SlV/TZVfUBEfhoL8wgA16hq/2V9fwfAOwF8DYBf634IIcSFWncKPuQSQjaKqr68kq8AXjNSOIQQshFU9QYsNu3KtDcUr2/B3o/kleWuA3Cdk34EwL4/vUYIIVNxULSOD7mEzJwT+77hnRBCdg9qHSFkDlDrFvAhl5AZo1A8opv9uDIhhEwNtY4QMgeodafgQy4hM2fT/yeXEEK2AWodIWQOUOsW8NuVCSGEEEIIIYTsDDzJJWTGKIBHuONHCNlxqHWEkDlArTsFH3IJmTn8WAshZA5Q6wghc4Bat4APuYTMGAX4BQWEkJ2HWkcImQPUulPwIZeQmcMvmieEzAFqHSFkDlDrFvCLpwghhBBCCCGE7Aw8ySVkxiiUX1BACNl5qHWEkDlArTsFH3IJmTMKPEItJITsOtQ6QsgcoNadhA+5hMwYBf/vBiFk96HWEULmALXuFHzIJWTWCB6BTB0EIYRsGGodIWQOUOt6+MVThBBCCCGEEEJ2Bp7kEjJjFMAJ/t8NQsiOQ60jhMwBat0p+JBLyMzhx1oIIXOAWkcImQPUugV8yCVkxigohoSQ3YdaRwiZA9S6U/D/5BJCCCGEEEII2Rl4kkvIzDmh3PEjhOw+1DpCyByg1i3gQy4hM4YfayGEzAFqHSFkDlDrTsGHXEJmjELwCP/XAiFkx6HWEULmALXuFHzIJWTm8GMthJA5QK0jhMwBat0CPuoTQgghhBBCCNkZeJJLyIzh/90ghMwBah0hZA5Q607Bh1xCZo3gEeUHOgghuw61jhAyB6h1PXzIJWTGKIAT/F8LhJAdh1pHCJkD1LpT8CGXkJnDj7UQQuYAtY4QMgeodQv4qE8IIYQQQgghZGfgSS4hM0aV/3eDELL7UOsIIXOAWncKPuQSMnNO8GMthJAZQK0jhMwBat0CPuQSMmMWXzXPHT9CyG5DrSOEzAFq3Sk4C4QQQgghhBBCdgY+5BIyaxb/d2OVH0II2X42r3UicpmI3CEiR0Xkaif/dBF5X5d/s4ic36W/QkQ+VfycEJGndnkf79rs875pzRNDCNkpqHU9/LgyITOGf0+NEDIHNq11InIIwFsBPBfAMQC3iMhhVb29KHYVgAdV9ckicgWANwP4flV9N4B3d+18B4B/paqfKuq9QlWPbCx4QsjOQK07Bd0tITPnEZWVfggh5CCwYa27BMBRVb1LVb8K4L0ALjdlLgdwfff6AwAuFRHb8Mu7uoQQshTUugU8ySVkxiiEX1BACNl51qR1Z4lIecpwrape270+B8DdRd4xAM8w9U+WUdWHReQhAE8EcF9R5vux3zD+cxF5BMAHAfyMqupqwyCE7CrUulPwIZcQQgghpM59qnrxphoXkWcA+FNV/UyR/ApVPS4ij8PC+P0ggHdtKgZCCMGOaB2PcAiZOSf0USv9EELIQWDDWnccwHnF+3O7NLeMiJwG4PEA7i/yrwDwnrKCqh7vfn8RwL/A4qOChBASQq1bwJNcQmYM/54aIWQOjKB1twC4UEQuwMLgXQHgB0yZwwCuBPCfAbwEwEf7j+OJyKMAvAzAd/eFO3P4BFW9T0QeDeCvAvjwJgdBCDnYUOtOwYdcQmaMgl8eRQjZfTatdd3/O3stgBsBHAJwnareJiLXADiiqocBvAPAL4vIUQAPYGEOe54F4G5VvatIOx3AjZ3pO4SF6XvbxgZBCDnwUOtOwYdcQgghhJAVUdUbANxg0t5QvP4ygJcGdT8O4Jkm7UsAnr72QAkhZAUOitbxIZeQmcO/k0sImQPUOkLIHKDWLeBDLiEzRhV4hF8eRQjZcah1hJA5QK07BR9yCZk1ghPg/8klhOw61DpCyByg1vXwIZeQGaPgjh8hZPeh1hFC5gC17hScBUIIIYQQQgghOwNPcgmZOfw7uYSQOUCtI4TMAWrdAj7kEjJjFIIT/Du5hJAdh1pHCJkD1LpT8CGXkJnDHT9CyByg1hFC5gC1bkF1FkTksSLyGyLyWyJym4j8VJd+gYjcLCJHReR9IvKYLv307v3RLv/8DY+BELIkCuCEPmqln12BWkfI7kKtOwW1jpDdhVp3ipaRfAXA96jqXwTwVACXicgzAbwZwC+o6pMBPAjgqq78VQAe7NJ/oStHCCHbDrWOEDIHqHWEkJ2n+pCrC/6ke/vo7kcBfA+AD3Tp1wN4cff68u49uvxLRYQfDidkKxE8suLPrkCtI2SXodb1UOsI2WWodT1N/ydXRA4BuBXAkwG8FcDvAviCqj7cFTkG4Jzu9TkA7gYAVX1YRB4C8EQA95k2Xw3g1QDwWHwtnvuol642EkJmzuNwxtOH1uk/1kIWbFrrDp1xBi74xZ9fZCiwby3x0mxea72sfEZLDLVyLW0PbWvZcbSOp0zDwL6i9oZcs6hsSzxDr/8q91IWSzSmZeMNXj/mvHOpdStCrWuMoVaupW1qHbWuFhO1biM0PeSq6iMAnioiTwDwqwD+wqodq+q1AK4FgG+QM7VSnBCyIXZp125VNq11p5933imtKxe5/ne2CGeXyVtsW+rV2sryWhZ4j9I8ROO25Wpxef3X6kUxDzHjGd61K1c6a3KiPlr6zsbRek9512Kd91JtHmuxrEGmqHWnoNYN6Idal0Oto9ZtKYMe9VX1CwA+BuAvAXiCiPQPyecCON69Pg7gPADo8h8P4P51BEsIIWOwMa2LzMnJjr1gnDpq8tSUa902jOrXYrB5LeupZzhajEA2J1ncQ+egbK/88fr0fgO5UfHKRacYtTHVxp+N3euzTG+5l6J+s2tWM3Hl+GvXntvia4Na1xBXmUetW0Cti+tQ67aGlm9X/sZupw8i8jUAngvgs1iI4ku6YlcC+FD3+nD3Hl3+R1WVl4mQLURV+C18HVuhdXbB6xdBz/zZ3/Z1i5Hz6rfE1ZI3xFR6JqqMKxt/y3hb8HbyvZ39FvMezVPN8NjrMuS6e+Vrc2LnLZqDaDyRMctOlKK4aobXzs1AqHWnoNYNiKslj1q3vxy1Lo6LWjcaLR9XPhvA9d3/33gUgPer6r8RkdsBvFdEfgbAbwJ4R1f+HQB+WUSOAngAwBUbiJsQsiYe2bCgichlAH4RwCEAb1fVN5n8b8HiS02e0JW5WlVv2GhQPtujdeXCnxkGSX63tN/KMm1GMUXlyjy7627NjZ2flnhKaictnpmrGb6+XK3vrJ1sDr1Yomtejqnl2lm8drM4yrSMWlzRfd+XX9H8bVrrDhDUughqHbWOWrczVB9yVfXTAL7TSb8LwCVO+pcB8FukCDkAKIATy6poA52JeisWJwXHANwiIodV9fai2P+Bhcn6JRG5CMANAM7fWFABo2idt7tb7iC37JxHi6eXbnfOMwOZpXnmqLZD7RnWKAav3dZYszZsXDUjFBmdzHiVfQyNuaVszTijMa9vv6Xv6B7zrk12vct2vXvS9puZy2VNbNHUJrXuIEGtS9KodXGZsg9qXdwutW5raPriKULIriKb3vG7BMDRzjxBRN6LxZ+jKB9yFcA3dK8fD+C/bjKgSYkW6BYzok5ayzqWLeJRGzUDU1uAPQM41HjZtMhoRv3b9ry+I4Prte3NV8s4Wq61bdvOVTTvQ3xMzcC1GKuaKc7q2N/eQ0nNePb1aqcoQSA83RgRah21jlq3PwZq3ahwFgghm+Tkn57oKP8sRc9PAvgbInIMi1PcvztOaFuEXQztwlYujJFxsGW9BdOWaSEzCTUz5r2OxhqdGpR5kSH13mfj9/KzNhG05bWZjQNOXm1HP4rB9ukRzU3flmJv/SGmO+qvdoJS9u+NU82PvYdaH3jIdkKto9Z5cVHrqHUbgCe5hMwYBXBCV1bRs0TkSPH+2u5PSbTycgDvVNWfF5G/hMX//fp2VT2xamBbR8182F1doG3xz3aIs11524dnGmsnCD3Lnl5kY64Z1MxYZeOv9TF0x75WLxtHOWfWNNbiqBmyqK3atfXuJRtjbd6ie8m7r1vunRVZk9aRVqh11DqvLWodtW5E+JBLyMx5ZPUPdNynqhcHeSf/9ERH+Wcpeq4CcBkAqOp/FpHHAjgLwL2rBrZ1RKYqW0Rri6K3+EfGLTIAUZxemcgA2P68sdX69sxaq9GIsOPv07ITE48WkzzEYGVzkZ1QtVw/79pE17OGdy2idqJrVKvnjT/LX5I1aB1phVqX902ti9vO2vdipdbtg1q3gA+5hMwYhWx6x+8WABeKyAVYPNxeAeAHTJk/AHApgHeKyH8P4LEA/miTQU1KtPB5i2a2YEa7xGVaZhRqu9Veu1H8tm9bvtV0ZEan1Rxl7fXla3Md1UXwPjJHrbSa0Oz6eeWW/acdXbNofmoGs3bP9GVqptJ73cgIWkcs1LoYal0OtW7/60aodafgQy4hM+fEBnf8VPVhEXktgBux+PNA16nqbSJyDYAjqnoYwI8DeJuI/K9YSPoP7fzf1m4xLTVadoltHzWz4O3+e+225mWxWmo74BZbtmVnvRy/TYd53WK07ftlr6ed+5Z5KMuU4yrzIvMamVwvhha8dr17qeXhIrqfWoxkhU1qHQmg1u2HWnfqNbVuf2zUurXBh1xCyEbp/ubtDSbtDcXr2wF819hxTYJd3GqLoFdv6EnBkAWztpPdcgLgbU94px2RQfDi9AxNFi+C9NoYhpyulPnZGKJTqGj8kYmz/Xp50fvWk5poDmsPKp7Ja5l3m6fIxw1Thmwv1DpqnddnWZdaV4+FWrcSfMglZMaoAo/wYy3bRe10oPVyeYZEgvctbXnvo8U921X36vV59rXXT2TUvPHUDImX581TFKf9Hc1py/WsmTM45bwxZfdPZiozc5dd02wMGV5M9j7xynnXsAFq3RZCraPWeXXglKPWNUOtOwUfcgmZOfy/GyPiTXWLSfIWe7uI2/azdiNTleGduLQYQq/PaKfbGjBg//ijvrP+vJi8uMv+ojqRObJjyMjKRCcftfqRsfXyvLlHkVeb78g4RrGVZW3dlrG15DVArRsRap1fl1rn51Hr2vMaoNYt4EMuITNm8QUF/L8bo5LtIntmDvCNSc24tRqlLJ6sfEu5Mo6sXmawvHxv7JnJW2ZsQ41xbWyZmYz6K+dgaGytZrklzWsnu0dbjKm97zJTavNs3w1Q6yaAWrcfah21jlo3GpwFQggZC89kRQtvuchli6zXR/nba9t7H+2KezF5+V55W7dmfryY7fhrdZcwBdUdd1t2GQMSmZuoj/6n1exlBtPeT+Vv26/33pv/bAw1M9lyPcv86KSObC/Uurgvr70+jVpHraPWrQ2e5BIycx6hik5HtKC1fBwr20EG9rdrDUu0Y213nKPdbNuuLRfFH8VrY/bee7FmZb04o/ysDzuval63/BPyDFl0HVrKDDmt8NrKThGysq3U5tDrJxtT6zwnUOsmhFpHraPWUetGhg+5hMwYBf/vxqh4C1jtvdcGgjKZUcqMYPTblm2pM8SAZca1RmQGojZr5sJro8UoDmnXlmsxulkMNRMY9dsyb0NlITLi1jBnDzfZv4Oy3hJQ60aGWketo9btLU+tGx0+5BIya/h/N0YnWujKtNrJRbbQlwujXVhtfiurGrXoxMPGItgfv+0rMwO18rZMLVYbV2RebJrXV2T67Tiy+LITnSGnHpn5ysYS3Zdemx5lG9lDio3X5i9l/qh1o0Oti2Oh1lHrqHUbhw+5hMycE4O3Mcna8D7mBLSbtdoi65Ut+6kZyNqtERkir8+o/czMlGWzGGtGpyxbSxtiSGrxtxq1oacTWVytJtLG05LWUr51/moPMtG1GXICZqDWTQi1jlpHraPWjQwf9QkhZCq8RTva/a4ZDgRly3ZsGWsCvdj6PA3Kegu6t+PtxZLtlpevM3NX9ldrp3wf9W3HWvbjtdVi2FrHn8VTxmDjap2H8n0Ulx13dH2i+fDShhrM8nXWz5YhIpeJyB0iclRErnbyTxeR93X5N4vI+V36+SLy30TkU93PPy3qPF1Efrur8xYROZjulVpHraPW+WnUur7O2rWOJ7mEzBj+0fAJqO3Oeot86y68LZvtGg/ZJbblh/Zd1vfaqC3qWazWlHl9trTVUnbotRtyUhDVqZ1W2PFG6dEpSBZXdt+1zmGGvV+yuMr3S5jATWudiBwC8FYAzwVwDMAtInJYVW8vil0F4EFVfbKIXAHgzQC+v8v7XVV9qtP0LwH4WwBuBnADgMsA/NpmRrFmqHX726DWUeuodaNpHU9yCZk5J/RRK/2QFYl2gmsGw6sbGQ/PcEWGwesjW9C9XfXohMLLL+PzdtMz7M5/ZlIyw+UZkGi334sh6887TWi55t7JR3RKY8tEpyBR3GVb1kRnfXoxekQnL5Ep9U5eaga4gQ1r3SUAjqrqXar6VQDvBXC5KXM5gOu71x8AcGl2WiEiZwP4BlX9hKoqgHcBePESQ98OqHXUOmrd3vrUOgCb0zo6VEJmzOKPhq/2QwZSLs6REfDqeO8jo1C7LJGps3HZ/qK8bOc825GP2oRJy3bSPdNg62fUzKI1yLb9yEi2mG/bjmfiyja96+3Vq825F0M53uxe8sqU2JjtPGSmN7rma5CZNWndWSJypPh5ddHFOQDuLt4f69LglVHVhwE8BOCJXd4FIvKbIvL/ich3F+WPVdrcXqh11Dpq3f44qHWjaR0/rkzIzOEXFEyAt+h6u/TZApu1WRLt+nu71zXDUOsj2qmuGZ2s32j8nuGKjGBt/iID7ZXLYrYm2iMyetkctl6X1nulZtpaHhyy8bfMYXT9amNYQa7WoHX3qerFqzbicA+Ab1HV+0Xk6QD+lYh82wb6GR9qHbWOWje8bktsCdS6BTzJJYSQMfF2x2u7yV79Wlot3+6ee4Yl6iva4a+ZnrJM/9ozhVH5skzNWA0ZR1mmtmNfjjE7BYrKReP3rr8319n4anMZMfQkITN15e/aNWipG12n7eM4gPOK9+d2aW4ZETkNwOMB3K+qX1HV+wFAVW8F8LsAntKVP7fS5vZCrdv7mlq3N51at7c8tW7tWseHXEJmjAL8uPKYZAavz4/qZYt77dRhyE52iwlY5tTFG0PrIu7VjU5mvJOKltMRz2DW5rd1/FFfLeOPTl680xE7H545tdfNM7URNdNtyebUthnVtfNnx93ICFp3C4ALReQCEXkMgCsAHDZlDgO4snv9EgAfVVUVkW/svswFIvKtAC4EcJeq3gPgj0Xkmd3/Z3slgA8NH/0EUOv2p7VArTtVnlpHrVsRflyZkJnDL4/aAmpm0FtAI0PT50ft2T6zcl7fWaxZe6UR8Xayh8Tg9ROdtLS2adMiw9XSZtZ2azkvP7sW6vz2qLVpKec1u24tc1Qz4La9ljoD2KTWqerDIvJaADcCOATgOlW9TUSuAXBEVQ8DeAeAXxaRowAewMIcAsCzAFwjIn8G4ASAH1bVB7q8vwPgnQC+BotvGj0Y36wcQa2j1tViyN5T65qg1i3gQy4hc4ansePSsqihKGMX1cjwRWWzNrN++9jKPoe0EZHt1Gf5dvxezNZYDjXMmYHx4ovKeHF79VY1k165qPyQa5bdS5nRi8bf5w29Z9YtSyNonaregMWfvijT3lC8/jKAlzr1Pgjgg0GbRwB8+3ojHQFqnd8Otc6HWrc+qHUn4REOIYSMTbkQZrvQmUG0i3G0u6/Ya4hsfdtmZFSi97atqF2bHq3B5ThaDHLWRtl3ZuiWbdtLy+au1nctvtqJRWT4o/SoTe9e8gyd10Zmemt9DzGn5GBAraPWUeuodRPBk1xCZoyC3648CbVdX7tLbY2MYL/BKxdn74TCEu38RybGK9diXmx6duoQ7dK3LPRZPN78LdMW4I+jpT2vjZY6rUayViZL9+4b736wMdVO2yJD2L/PzH8UwzInJaDWTQa1Lm+TWketo9ZtDD7kEjJz+HHlCSgXs2jXvdWg2PZq71sM21BjFMXbssNt+xtiqGz/LeZvSCxeTEPL2fI1E2vnbOh9kLVZu/ZejJEh9/pb9nTMtuG9X4NMUesmgFoX902to9Z576l1a4MPuYTMGAXFcHSi0wlLzayUi3G5IEevvcXdLubRaxtXFG+UlpmXzOh5449OX2p1vXED++cga6M2P1kbEdn8e21m44/GlM2vZ7qjMS5r4FvMa9TfMu0ETVPrRoZa54/Da4NaR62j1q0dPuQSMnMohltEtNBGCzTgm4+WBbJlt3pIrNkCnbWbmdrMwEanGbUTFS+uFtOXvY7wjGNLPN4Ofzb+aEyZecwMXtlWZEZbT0xq1Iz/GuWJWrdFUOv21qfWnXpPrVsZat0CfvEUIYSMSbkge3k9Wvy2i32Gmt9entdfFovXjmc4W9st46i9jsafmbsI22bU9zLt10xd2abtN+o7Mns2Ps8o2nq2fe+UxDNeZb92TrL5s/15ZPdqGWftnifbCbWOWket299+VI9at3Z4kkvIjFHwTwhNQsuUeycHXlrNzEUmJ9r9ru0o13a1ozLZLnYUu23PMx212O1YI0MZmZ7I6GTGtzbWrK4l6tMzgNFDRTSGMs2eWERxeSccWdlojofeR2UfSxhBat1EUOt8qHX7odZR69YMH3IJmTn8Fr6R8RbDmtmy9dBQPisXtRMt1HaBj+KNdtGtafP6s21kBi4ywV6bLSYrSrMmseyj1mdkdOCkZUarTPfMnm271ZBaWu/Bss0s9lbT57Vbi2tJyaLWjQy1zu/PtkGt2x8ftS4u1wC1bgEfcgmZM8r/uzEqnpGIjCCQL/Aoytj2PTwj45Xx2hlisjIzZ9vy2pWk/JC2LbZd257XTmRw+vc1E2r7j9Ja4o7Ktponmxc9dNh70KZHMdpy3rxm91XWbnT6MgRq3bhQ66h1Xhq1Lm+XWrdW+JBLyIxRUAxHp2V3OdvlzRZ4zwR6Rsvrx2vfMwc1o1czpi072VmexxBjGtXzjEVkyKO8MvYsJi+WZQy1Nw5L1G5maqP4swcP72HGu5eAvfPoxWnTaua3AWrdBFDr4vipdXHc1Dpq3ZrgF08RQshYRAtftqjb3WHbhhY/MK+zhTaLsWYQh2DNUmTsyvRoR9zru8VIRzFHc2tjaTEunmnx8rLxZzFa4wkn3ZuryPiW94mt491LtTmsXdOa58rusezfB9lOqHXUOmqdD7VuNHiSS8jM4Y7fBHg76CXRLnTUVl8nKlfbHc/K1E4IvDJefBk2dq+/obvtNq1W1jMe3s66jSGiZc6zul5fmRmttVHWs/nLXrvIDHqxtcRs+13Hg0cBtW4CqHX765dlqXV+X9Q6at0a4EMuITOG38I3If0iaHfNowUeiBdCb8d8qGHLFnEbi23XOyXIjKjXh52LrC0x9bzfkcnzytu2vLqRYcnMje0zwrsOmXnKxu/F4cWa4d1LNj8zw9npT8t9mN3LS8oVtW5CqHX726DWUeuodRuHD7mEzBylGI5PuTBG05+ViRb2VRbIaMH22s3MCZx8G09k1mrm08YSla3t+tu5tf1mZtzWabkWNi0zqbV5qBnQ7FoNIbqXovaifr38rK0Wc70k1LoJoNZR66h11LqJ4P/JJYSQKbC7332aOvlZG6UB8Wgxc3YB9naz7ftsES/zvbLWUNR2vL269jQkMg/ZTnsWT80wRW3YGKP2vTZq827b9uao5VpF/dXSbfvZw0nLvWVj0+B1mUbvdvCg1lHrqHXUugngSS4hM4d/T21kIsPl7XRHZDvLtbJeut01z2KomZmsbJluF/zW23Do4h/17e3A1wxSbee9Jb2MwatTG1+W740nO2WJWFUSsntJnDJlnvfvYMi/jQRq3chQ66h11Dpq3YTwIZeQGaPKLygYFbvYRa831e8yBtErZ3ehW9ocap5qbQ6Zr5aPzpXlbPutRmwZkxXF6fVbM+fRSYFtu2zDe/hojTMz1cvgjT/76N+Qpql140Kto9a1xun1S607VW5o09S6k/Ahl5CZw/+7MSKlQSjflwxZfL02vXY9s+CV99qs5bUYnczwjoFndHqieYjMXM28Z9fWq9OSXzNULadFtdOG8nVtTqK5seba66/FHHpxtDy8VKDWjQi1zo9v01DrqHWg1vXwIZeQWcNv4ZuE2qKXLeTZR5xQSYtiqLWZmZ1a39mC7RlLj5YTBGB/W974S1MSzVEZb2SUvHFlhsujdo2WPS2JzFyLIcvyagY+aiM6XcqwJxwI3jdDrZsEap1fLoJa117H9kmtO1mRWreAXzxFCCHbgCJf7DIi87JKLNYkRIttZCBt3ZqBse1EfXhG0Ro0z7DZ9uwYvT5bDEZmXFvNrPfe6yPrz54EeL+zNloeErwyUR9lvq1f9h+ZcDHlM8NJDhbUur15UR/UOr9tah1phCe5hMwcfqxlRFoWx57sZCHaqR66exwt5N4CHZmUFkPotZOZ1ZadcNt/q4nyDGk0Ps/8RXNTG6ON08urlanNadR/ZJqyeyQ7vfHyvHiyurZO+Tsb91InG11Vat14UOuoddQ6at3E8CGXkBmj4BcUjIqdamswMoNj0yKzmC3Ytl8vz4vRi7cWSxaf10/UR2s5mx6NsWWOh4wvq+uV88YfGZss7tb+LTUzVrZj8zND6fXvvY4Mdpnm5UXjbYRaNzLUOmodtW5/LNS6UeFDLiFzRgFdYbeQrEC0AAP7F+OW04q+XPnbqxudjETlW0xhNI4hZrbMy4xq2a/Nt/1nMURxeCc3mfnM8iNaxh/1l+VH1K6DZ/K8e2lI2y15XpnIVK7q2ah100Gt86HWUeuA/fWodWuDD7mEzBz+PbWJqO28txgur83aaUW2M23baD0VyQyr134L2c7+kPet5m5o3dZx2DlpqZcZ9SyuKM+OqcUM18q1mN6sjaxe64POQKh1E0Gty6HW7e2PWrcy1LoF/OIpQggZEy1+gHgn2Xtd1vF2uVt3k6OYyjZqhqNllz2LZUjd1p32mpmw5TLTNySurE/vmkTjt9fBGsaaySrHFd1LQJsRq5mv6B4s+xt6GhPNFzmYUOuoddQ6H2rdKPAkl5AZowC/oGBsshOIlhOHHm+BLBfZ2glAZISyhbq2g+6dfkQnJdFJiDh52S551L4XH4JymemJYqzF4F2vsm42X17MXr9e37VTh6y8N+f2tY3Lu8bRg0MUg00D/PmLDHsD1LoJoNbF5ft0at3+mL1+vb6pdS7UulPwIZeQWcO/pzYqQxYsuyB6i2BmYiIjlJmpFgMatWXzvTw7Fs8IWrzxe+mtt7EXY9RmND9R3y1la+WzvCzWDG882fXP6nn5LeOxfdfK92TzNwhq3ahQ66h1tfLUur1Q69YOH3IJmTn8goIRyUzKkFMJb/Gu7e5Hu+0tpyZZ/ZayPbV1t2YOo7ZrYy/LZOMdah6zufd2+lvii97X0oeMoWYaW+Ity7XEWuur5SGkNa6oK2rdeFDrcqh18ftaOrWuCrVuAf9PLiGEbAPlQmwXKO/jU+X76JTD1rPvJUgv+80WblvXthEttFl6ZorLPqO+y/LZ+Mt0r86qsVhzGLVhy3t1WvvMjI33wODVyU7MLF57LfXKflsehlraIwcHah21rsyn1lHrNgRPcgmZOfy/GyOSmYrMZLWcCmTt2EVWnddemhd73062UFtDYNutGYKW8dfmy/YbtWPxjJW36x7FUoutJY4hJwvRe6/N7D5Yps2Weq33VxTz0FiyMKl140Gto9a1xEGt2xvnkDpZmNQ6AHzIJWTWqFIMR8WaipqhQFLeK1Orn5W1JsmWL/MjY1NrNzO7UVlv/OV4svZsfa9MZhBbjIpXt8UU19ot8yNj66Vl90ptPFndWro3fu8EpcVcZm3a9hqh1o0MtS6Pm1rnx06to9atET7kEjJz+AUFIxKdPNRMTvk7arM0Cd7CaE1TtMi33A4thsgzdZ7JytovXw+N085DZC4i42pNxxDj2pIfjT+7R2qnJrafVoMXjdWar6wvr0w0n1EZG3M0/iUli1o3ItQ6ap1t30Kto9ZtGP6fXEJmzmLXb/kfMhBvV7pMtzvBNayhi/IV+xffFiNSUjMy1mjVzJ1Xvzb+FsPoGb6+7VYj2WoKMzxTE7Vv4/PaaZnT1v5tOc+kt9T1yO7DqEw5fu80xMsbyKa1TkQuE5E7ROSoiFzt5J8uIu/r8m8WkfO79OeKyK0i8tvd7+8p6ny8a/NT3c83LT8DI0OtO5VPrdubR62j1o2gdTzJJYSQsWjZye3T+7SyTPQaQfmo75a6Xgy1HfXMvNgYbB9e7Nn4LdGYor5t/zWsAbEss3HujT/rM+rbmqKobDa3luh+iOasdt28Ol5/Xlve7y1DRA4BeCuA5wI4BuAWETmsqrcXxa4C8KCqPllErgDwZgDfD+A+AN+rqv9VRL4dwI0AzinqvUJVj4wykHVBrdsLte5UXWodtW4kreNJLiEzR1VW+iEDyYxLVN4zVC3myp4WeH1mZSKjNmSX3OszKl/GVJaJdr6jWEvDEe3SZ2Yk250v4/FituleTGVay3zZvrPYvAcI217LfHrxeHNm76HW04fMPHtz2Npuwoa17hIAR1X1LlX9KoD3ArjclLkcwPXd6w8AuFRERFV/U1X/a5d+G4CvEZHTVx/xxFDr4vJlTGUZat3evrPYqHUh1LoFfMglZMYoVhNCPuQuSbaLa197Jwt2MWxdFKMTC7tzHJ0oRDvMa1iUT7brjTWKcSgt5jnb8Y/qlacTnlEfeoqSxdcy/trcZPNZM4LemOz4bbuewY/68u69lmtSYU1ad5aIHCl+Xl10cQ6Au4v3x7D3hGJPGVV9GMBDAJ5oynwfgE+q6leKtH/efXzvJ0TkYIkutc6HWketo9ZtXOv4cWVCZs661mzSgN21rUm4LROZIGs6op3+mrmsxVgzJVnstn4L0fhb+y/LZfNSptVMR2ZksvLL/kOLxpmZ1swsRWa05V6qGVhbxzuliK5HNo4ybQXBWoPW3aeqF6/ejI+IfBsWH+t7XpH8ClU9LiKPA/BBAD8I4F2bimFtUOuGQa2j1lHr1q511ZNcETlPRD4mIreLyG0i8iNd+pkicpOI3Nn9PqNLFxF5S/efjT8tIk9bJUBCyAZR8CS3YxSti4zHMobIa7tm+mq7ya3Y8pkJsmmeGfQMwjpoGX+5O+/VbRlrZk7s9cjM8Cq7+LbvmumOzGuf58XiPXSUbWfx18Ze1o/mbMX52bDWHQdwXvH+3C7NLSMipwF4PID7u/fnAvhVAK9U1d89Gbbq8e73FwH8Cyw+KrgS1LoBUOv2Qq3bG0PUJrVuK7Su5ePKDwP4cVW9CMAzAbxGRC4CcDWAj6jqhQA+0r0HgBcAuLD7eTWAX1o1SEIIGYHxtM4zZ0Pq2V3f/reXbutlsfTvbfkhpxoemZnz+o/MsDVVmWH05tgbTwu18tYEtcyzl56ZSu+9d+0jc5WdJGTtRTFE90R0WuQRmepsvtf5YLBebgFwoYhcICKPAXAFgMOmzGEAV3avXwLgo6qqIvIEAP8WwNWq+p/6wiJymoic1b1+NIC/CuAza4iVWte/p9bthVq3vw61znJgtK76kKuq96jqJ7vXXwTwWSw+a305Tv2n4usBvLh7fTmAd+mCTwB4goicvWqghJANoSv+7AijaF25kJbzl+28e/XLOtbUlOllH9G1ik48vDjsAu/dC3Zctb68dvs0r36LgYriqrURGcnI2HjXsXXMtf7LvmvGu9XYZ/9+a/dSeT/ZE4eW8dv7wzN5mcnNxtjKBrWu+39nr8Xi20I/C+D9qnqbiFwjIi/qir0DwBNF5CiAH8Oph8jXAngygDeYP59xOoAbReTTAD6FxenI21aYgT5Wal2URq2Ly1DrqHU4WFo36P/kyuLvHH0ngJsBPElV7+my/hDAk7rX0X9IvqdIQ/efmF8NAI/F1w6NmxCyJjb9kWMRuQzALwI4BODtqvomp8zLAPwkFvL6W6r6AxsNqsKmtO7QGWd0iabDmhnxdvVrO8cnAyh+2wU3aiPaZbd1sz6tGVDsjyHrw6a1GGObbs2Mbd++z04avPlqGUfLtfbMlDVdtXhreGbLttVyXW1btRg8M2ljsOOsXbMhZrqsumGtU9UbANxg0t5QvP4ygJc69X4GwM8EzT59nTFaqHVJHLV4y77619S6uE+vbWpdHBO1bmWaH3JF5Oux+I/AP6qqf1x+6VV3BD3oUqjqtQCuBYBvkDOXvIyEkFVp+cPfyyINf09NRC4E8DoA36WqD8oa/gD4KmxS607/lvP8utGCD5NuX2f5Xhte+eoA4BsezzzY997rKM5sjEPMYsv4h7YJDOsju1ZDrnUWQ9SnzcvatqcHmeluuV+yWLz+4by21ym7hwaySa07iFDrgrapde19UOvifGrd5DT9CaHu89EfBPBuVf2XXfLn+4+rdL/v7dJb/kMyIWQetPw9tb8F4K2q+iAAqOq9mIiNa51nPvYEUPzud3E907UvcOxfdL1dZC/PGkIbV7TI2ri8upbshMKOGfDHP6S9skw0j0NOTaI+yvTMnEVtlXOXxTWk3b7tsl0bt8CP2buXbIzRvNn47b3knU7YtqO4yNqg1jlxUevqfZTp1Lq9bVDrto6Wb1cWLD5b/VlV/UdFVvmfiq8E8KEi/ZXdt/E9E8BDxcdfCCFbhALb8PfUngLgKSLyn0TkE93Hm0dnNK0rFz67ENq8IQtfZtD6tsvfZX5LWmQyIpNgzWXL6Yo1DutY+GtGKjJqrTvhNVPnxWDH5pUr741VTeAQM1mW99KGnjDYe9wzmdE8WYbMh1eV3yQPgFpXTaPW+VDr6m1R67aKlo8rfxcWf6vot0XkU13a/w7gTQDeLyJXAfgcgJd1eTcAeCGAowD+FMDfXGfAhJA1ogBWF7RV/57aaVh8a+ezsTgh+HUR+Q5V/cKqgQ1k81rn7f726XbXuOWEoS9X2/3v+/VOMrJ2okU5OiXIThZaY/d2+Gumtqzbp9v5E/jteWmlWYmMYWZ0o9dRO5kRzsrUTiVsuve+1kcUb9nPkGtm+7T1bV+2Xq3tiPVo3a5AraPW7a1PrYvjLfuh1h0oqg+5qvofEU/zpU55BfCaFeMihIzEhv/vRsvH3I4BuFlV/wzA74nI72Dx0HvLRiMzTKZ10cJbMxxIykV5dtH2DED/ehmzZeO0RtaOLTI7UV7NiGVGotZfLT0ye5ERbnkdmcfITNeMWav5XNZQReWiB4H+ffZw0L/2xmqxxnMg/H9qC6h1oNZl6dQ6at2O0PR/cgkhO4yu+JPT8vfU/hWAZwNA93fSngLgrpXGtM2Uu+bRQtqyQNVOKuzryPCV/ZdmrXWR9Exi1kcLQ/qvtW3Hb18P6SMzsVm/tf6ifGumW6kZzjKm2rxk48jGFt1L9rV3r2TGvWZ+MzardcRCratDrdvfNrWOWrcm+JBLCNkYjX9P7UYA94vI7QA+BuB/U9X7p4l4BLLFvGYm7AlE1ke2ULUYTi9OzyzUqC3SkSmonUS0jC/qr9VoR6YnMiteH9lJkS07ZEyZmcr68tqK+ujT7IOKVz6bj9pcZKcvO2S4Zge1Lu6TWhdDrSNrYtDfySWE7Bqb/5KBhr+nplj8sfAf22gg20S2I14zd94Or2cKot34VuNY1otOL6J+bRnbvzfGqG4Wo9eujcdr28YbGS8bj1fGli1fZ3Nkx++NtWYEy/jtaxsHgjxv7N4ctVzvKE6vjv0djT+by0Hs1heqHBioddQ6ah21biL4kEvI3OHu4biUC69nqPrf3o51zSS0GKdWw+eVjcxmi0G0ZWpm0rafxV4zJVH5bHyWmgHx0jOTnhk9z5QN8Sy1614zacvcS9E4ozFl17lPj+ZtWc2i1o0LtS5um1rnl6HWUevWCB9yCZkzCu74jUl26pC9jwxAtvNcMxaZmcgMk1e2xTCV4671X/bnjb/1dMdrL+pjmd3z1uvnjbv1mg3p1/tty9Tm0osjartl3oeY18xQt5ykZFDrxoVa19Z/2R+1rr1fal0Mte4k/D+5hBAyFtGOsTUG1qDYfJj3Q9azmiGKzFFpHKKdZi/2ss2+brZbXqvfYgCGGIxWM2F36T1TF5keb86y/hT5XNqYs9+ZmfLqRFij6eV7RtDWH3LC4MXJE4qDAbWOWketa4datxF4kkvI3KGQToO36+vtWvdlS+OVGSlrULw+vTa9uMr3Nq5oZzsyrVndWp43/ppxqo2/Fa/f1p39yPjUxh8ZyP517VSlFpdnorIHgiGxR/20zL33b8G7t5eFWjcN1Lr9UOuoddS6jcOHXEJmzzKrIdkI2YKflfWMwTL99URmNDOK5W+70NeMS9T/srdmrV5mtmtEY6sZ1ixt1RhbTLj3QFAzj9m9lF3jWrmWsbWWGwS1bmug1sX9DIFaR61zodYBfMglhHDHb1qiXeXaDrdXbi2LYxJPyw549HrZE4nWNrNTnVpb2SlEVNaaK+86ZKchMOW9WGqnHbXTgBZjZuPIrtcQ09jSR9+ON/7a3C8DtW5aqHVt+dS6OEavTxsTtY5a18GHXELmDsVwGryPR8GktZSvLfi2jcw02NdeW9nC7eVnBipqOxu/jSlLr8VekhlIO7fRXNeMWvnea6N2jb0xtBre1no2vWVehjyIeNe+1bzX7u8a1LppoNZR66h19T5s3Sy/BrUOAL94ihBCxkWRGwTPCNo8Ly0zS9HucZ/nmTLbf7YoZ0bRMzFZuQjB/jmQIN2ePAw5+aiVs23bdK8da4y9PiKDnMUWzUkUs23H3gPZvRSleebU5kenEtlJTlnH1qeBOxhQ6/JyEdS6/WnUOrIEPMklZM4oAH7V/LjUDFfLDq63250ZjyjNq9vatjeGlng3ebtFO/OtdXvj0XLak82LpWbMvPQhpzNZvLbtWszZuFrTor6GGEsvjiFz7sVCrRsXat3moNZR67JYqHUA+JBLyOxR7hSOS7ag10xHj7e7vsyi6BkCr++aCbJttPRTvu/LtBpOr62hJxhefN51iGKvzUUUry0btdlieLy6tTiiubdlvOuUmc6We2nI/TPEfDZCrRsZah21jlrn90utGwU+5BIydyiG41IzS0PMhLdge23XTFSf75kw21b5G8jbHrLT7pmeIeP32vHGX5sH77VFi98t5ihqM9u9V5MelY9ibTGKUV0vPbu3au17cUb3RmaOlz3ZKNsh40Gt8/OpddQ6at0o8CGXkLnDj7WMS7bYD6nfYjiihbSsB5OvJj9btG3bmeGqjaF1HiITEL3PYvRizvq0sVoi81nmlfnRNYv6iAxs6/1TKz/kPqzF3rfntd1iviMDvIp5o9aNC7WOWketo9ZNCL94ihBCxiIyUmV+VC8yEFFbZX9ZX9kudrY4L2M2PLPVglfXMwXl775ezWBF76M5z+pEZaK+WsYfmVs7foVvkmwfnkH20j3U/K6RzaltM6obPSCQ7YZatz+tBWrdqfLUOmrdivAkl5CZIxTR6amZwZbTAZsftWf7zMp5fWexZu2VRsTbBR8Sg9dPdIrR2qZ3gtESg1cma7u1nJefXQt76lEzVNF7Szmv2XVrmaOaAa+dCmXtNECt2wKoddS6WgzZe2pdE9S6BXzIJWTOeLuNZHMMWbBadua9neq+TnZS0rJwt55eRO1HptGLKzO8WV5mDloMo2eyvBOYrHz0fgjLnnJYUxbF6pnCyCTavFYTmd07tXspuyeWvQ+9Pqh140Gto9Z5UOvitql1a4cPuYTMGgH/78YEROao5VSiJDMfQy9rbWe8NdbMKGRx2TpZ/N4412XErOGIxuyNPzIpteuSzYnXV9Rebcy2ftZerd3aqVvtlMWWz/LXArVuEqh1fvst5Wweta7ejs2n1s0a/p9cQgiZgpaF0atTKx+ZNo9sB3wZI+WdDPS/1aRFsdSMwtCYvHp2HltPXLxThLI9gX+NvHmIYizza6cMtfFHsWTpXkz2te3bG793L3mnKLZ9O34PnlIcLKh1fvvUur39Uev2Q61bCZ7kEjJ3KKLj07I7ne2Cq/Pa21mu9e2lt5iJaEc6G0vLyU3L7nxkGFrq9+UzU2tPMLw+bVrrdWs5xchii8a/7LX20srxZw8CWYxZP1572fiRpA+FWjc+1Lphsdl2qHX763pQ6/a3Q3iSS8js0RV/yHCy3e+oTERkTGwZr/3aiUPWX0l2opKdeGT1s/77OtZ8RW16hrmlj6hPr4xHNlYvvqxNbx5Lw2TbtHVrcUcGuqV8K6WZLNvIxm8fcrKYalDrxodaR62j1lHrJoInuYTMnR0StANBuVgP3cn2ytn30QLute/11WKOhuzSRzEsexLR2laUZo1E60lFFovXphdf1I81cC0nKl770SlNZHyzWMs2bHoUv+23Ni9eeZvfYkJbodaNC7Uub8NrrxYHtc5vs0+j1u3ta+bwIZeQOaMAv6BgZGoGwyvrYU1L/3uIKSvTM0OStYfg9zJGsnZKE9Xp30fteGktRsTrL4qlnMdWw27Ta2P38of2Zct4pq3P89qIHli8+2/oPV6LufV0KqpLrRsXat1+qHX786l1bX23Qq07CT+uTAghY7LKDqsWv1vMldeXt8hHfUlQbpUTkLKNodg6npmr9eGZLFs2GnetXWsay3LqvG9ptxZH2W/0uywbtZudSJXlo7nJ5ixKz66B956+7WBBrWvrv6UOtY5aRwbDk1xCZg7/aPjIZDvpdqe5tsPunSbU6nq78N6O/pC2l6F2cuOVGXI6UOsvOkmptZWNf+iJSGtabc6t6aqNNYsnut5DTipqp1StJzJrvueodSNDrdvfl4VatzeNWketWyM8ySVNyMXfjgev/Es48d3fOXUoZN3oij9kOeyCOmRx83aS7Y56aeq83e6yHRsHgvJR3zVad71bDJxXvmzP7tJ74492+jd5P1vz1NpfVsZes+xUonxfO+nw7qUhePeS11fLCVutvSFsWOtE5DIRuUNEjorI1U7+6SLyvi7/ZhE5v8h7XZd+h4g8v7XNAwG1jlpHrcuh1m1E6/iQS6ocuugpeMY7fhO/8X/+Er7vn/17PugSsiqeyWvZGS8XoGjh9gyGJHkW77TDay9bCLOF0saaxRHV9cxM1l7t5CIzQTatNn/e2COzkxmsckw1wxhdV3HyvLEPuZdq2Pp2bEOu6zIPGRMhIocAvBXACwBcBODlInKRKXYVgAdV9ckAfgHAm7u6FwG4AsC3AbgMwD8RkUONbW431Dq/XdtHVJdatxdq3eQcJK3jQy6p8pWzH4ef+sbbAAA//ITjeOhbHztxRGSdiK72Q5Yg2q2N3kcmoNyptoYjMjnWXJSxDNlNHrrjXDNnXluRaWjZmY/ee3VqeHPktdXPZxZDZNBsP9519trITGP0eqh5jfrz5sLei9FJk60T9WPnrvGkwWPDWncJgKOqepeqfhXAewFcbspcDuD67vUHAFwqItKlv1dVv6KqvwfgaNdeS5vbDbUuL2PfU+uodVkbjVDrFvAhl1Q5/XMP4G/+wXcDAN74R9+GM2/74sQREbIDeAuJNVOZefNMRoupispHhqRlkfXMlz1RqcVmzU4ZU2ReI6PlxRQZGdufZy6s8Y7MR+0UopWaAbLjL3/bskNPJ6I5j/JbHjb6ct4pU62foXnTcQ6Au4v3x7o0t4yqPgzgIQBPTOq2tLn9UOv21qfW+e1Q64blTceB0Tp+8RSp8sjR38O9P3A+Lj3/Kpx+zx9DP/uZqUMi64RfNT8NpYGwu9r2fW1331tArZGqLcaRSfFi8XacbZutO+llG96ufjZ+a3rK35GJzMxSVtczpkAcr3cdI7zr4PVp2/XG7cXhxZpRM7DRfNg59dpsuQ+ze3kVuVpd684SkSPF+2tV9dpVG915qHX726DWUeuodRuHD7mkiYfv+n2cdtfv45GpAyHrJdqlJZulZac9KxMt7KsskNGC7bWbmRM4+TaeyKzVzKeNJSqblbH9eP1mZtzWabkWNi0zqbV5qBnQ7FoNIbqXovaifr38rK0Wc70M69G6+1T14iDvOIDzivfndmlemWMichqAxwO4v1K31uZ2Q62j1lHrqHUTaR0/rkzI3NEVf8hy2N3vPk2d/KyN0oB4tJg5uwB7u9n2fbaIl/leWWsoajveXl01ZSLzkO20Z/HUDFPUho0xat9rozbvtm1vjlquVdRfLd22nz2ctNxbNjYNXpdpQw2srb85rbsFwIUicoGIPAaLL1c5bMocBnBl9/olAD6qqtqlX9F9I+kFAC4E8BuNbW4/1DpqHbWOWjeB1vEklxBCpqLc3bbmz9utRlHG0mpQbJtZX9FCa3fny37KvDKG2phsrFlZr4+WNstyNqYIb6xZXtl+LX7vvW2zlhe1h2BsLeNpnZPsHmq5bll9Ly96MNkCVPVhEXktgBsBHAJwnareJiLXADiiqocBvAPAL4vIUQAPYGHk0JV7P4DbATwM4DWq+ggAeG2OPba1QK2LY6XWUetsHrVuLVrHh1xCZg6/IXlkMoPi7VjXTEOLYaj1bduqLfyrtFG732rjj/oYYliz1+V7azqy+YtMk3cNLNl4MgMUzXXt+kRjb72XvDQ7/lr89sEgirt2Lw5g01qnqjcAuMGkvaF4/WUALw3q/iyAn21p88BArcuh1vlp1LqVodYt4EMuIXOHD7njYRfSaEc3qldrF0G5aFHuX2eLeGaQajHZ9ltOE2qGzWvb1osMbzSOWl/evNXm3JZZ5hpnffdpNZObjce2YeeyxeAPuTez8rVyQ+6hLBYyDtQ6ah21jlo3MXzIJWTuUAzHI1qwaotozWwNMVM105T1V1t4I+NQM07eTncUTxZfy/zV4mgxtvbfTGu/UZ0hRtS2u4rpyx4IakQPB5Yo9tp1G2qyW6DWjQe1bm9Zal1eP8qj1i0HtQ4AH3IJmTWNf/ibrJuWRXKVXdyyjb5db3d9mV3mGna3PBvruu696OSg5RTAi6llXmw9r8+oTNlXmeeZ4ZZrUevTO93xykfptsyyRjoiK7+Ofweg1k0GtY5aR62L46z1twTUulPw25UJIWQbsIuSNUbZoqVBfc8MeXWi0wV1ytYWYa/NKPZsZ92Oyb6unYjY9iNqBq3lNKjWZlauZgC9dr3xZ4beM3xR+54ht223GigvJsX+e2/IvUXzdvCh1sXxee1R6/a3E/VPrSMFPMklZO6s/kfDSSvlYtdqOLx8W791d72sl9W3ZWo7zOXuubcwW+Phxef177XRv64Z4Wz8rSc4q55oeO14ZaN7IbsuQwzQMnNV1ovmspxPr1yt7ZbTmmVNpwe1bjyodX6fXltRG/1ral071LquXWodwIdcQgh3C8ejxVCVJqp/D+RGwevD1s3MTmZwonKRQbRxRGbTlvfiz0xFLV4bSy1eG6vtIzNvQ/yEd31q+dYAZeWy+LyTp/59zWxHfUXlMry5LNuyxrKss4p3o9aNB7Vuf9tZ/NS6U2nUOmrdmuBDLiEzh/93Y0Ts6UbLYuuZj9ouvdd+VM7mZf1HbUd1vX7L9FpfXv1lTVZksKPyQ9vq07N+svHXTGDG0Gthr73Xt3evRe1mpxG1eDMDaP+9rAi1bkSodXvTqXX722iNp1aOWre/W2odAD7kEkIohuOxzG54Zggi42JNYbRrHy36tv+onaHlPbNaI4rPM72APzfZyc6QNoecAnh52fiHzGv2vrWcF1cZx5B7ZOj19OLMDPMaTN/JNsk4UOuoddQ6at3E8IunCCFkTFqMVolg70Id7URb+sU1Mzmo5EX5tfJZvy0Gx5b3Tg3KfqITlyiObNc9O4GK6g0xFKu0Ec3zUNNXtmVPJqJ7KcvPYveup5eXlSMHE2pde1t9eWrd/jLUOrIkPMklZM4o+LGWsbE75eXC65kUuzBb7PWzu8LR9fV2tG2MXp2WtlpOEiK8kxmvXRtPZPCy+bHte7+z9r04IjKj3YI3lmy8tfnyaDkFajmZqPU59J5YtV5Xl1o3MtS6HGpd3Ae1jlq3BviQS8jcoRhOgzUYmUErqS3G1lja92U7HtHi32qYbBtDF+qaqfDiaTFlyxi32jiiubF1lzEsteu8zOmRV9eb69q94z28DDF4qzwYrAK1bhqodW39U+v2xkmtWx5qHQA+5BJCKIbbQbYYZgYkWni9BT5qZ2h6S9/eOLx2h6a3mJ/I4Hj3emRiVo0lm4Nszrz81j4zE7VKGzXDXLsPsnhajHVLey1Q67YDah21Lsqn1tXba4FaB4D/J5cQQsZDzW8vDxi2gNp8rx3F/t32srw1RVF8WuRFC7Xt32u3ZggsLQbMlrfjtUYjwpoiG0+LQcmuZZQ2JD/qM6vnjafWRqtRyup593zricYysZDtgFoX1836pta19UmtIw3wJJeQmcP/uzEiLbvhltpuOxrztFLWmkJbPjKRtTEM2eX3yg45qcj6iMbvGWI7tmx+vbqROQTivj0is+qNP7u3vFij/Civlu6Nv8X0tVwzOwdLaha1bkSodXnc1Do/dmodtW6N8CSXEELGxC7QPdGOsGewvDa91+V7r53sJEEQmw7bhncKYPto2c22MbSMv+w/MzuemfPa9+YnMpq2TDa3kSG0bZVz6RmiqP8yL4ojmyevbnZfeSY5u/ci0+7dS/3raM6G3EtkOqh1dah1+2Px2izTqHWkEZ7kEjJ3uOM3HaUJtItsbRe4toBn5co+7fWP+sl2070+UeStYzc9OgmpGZLWnf7aKUpLjF567WTE6y+al2z8fR1vnjITtex9tey9Z/PLtqL7scxbxfRR66aDWketo9ZR60aGD7mEzBkFP9YyBS3GzquTvbdtRIv9EKMVxVqLt8X02LJl+ezUJevf7uRHZsfLt79bDGQ0X7ZMy1y1li/7i2L1yrS25Y3XXpeoDODHn5nR7IEiqrcM1LppoNbtL1uWp9ZR66h1G4MfVyaEkKkoF6Ihu7Zqfntkpsgr5xmQbKGNzKMdU22x9eYgM2xe/+WufvkT9VPWt2YxuiZD4vHybNu2z6H3QhZbi8HxzGIUY3aPeO3UHlK8E53oPohiIQcLah21jlpHrRsZnuQSMne44zcNyxiJbJe41mZLGS8/MhWti3I2jtpvr55ts6wz5BTBS/d287PTDa+tbGfeizk7BcjGa997cUbX3DuVsH21Gqxl7sfoHovGlt1zQ6DWTQO1jlpHrdv7nlo3CjzJJWTu6Io/ZDi10wm7497/znaxh1yLqJ6NK9rpliDPq2fLeePIjIodf83cDjlN8d57ca2T7DrW+htiPmtl7Gtvzpfp37aR3V/962z8nrFdFmrd+FDr9rfp/bZtUuuoddS6leFJLiEzRgD+342xKXdsWxf6FnNid9bLvKiOjaF2wlCmRXk21uy0xDtNyNI9ItMYzW/t9GZIWkY0P5np8q7H0LiGmLMojuhear1vvTaicWUnN+u4DkU1at3IUOviOKh11Dpq3cbhSS4hc4c7fuPhGSdv9zdbNLM2vPRsN9/b0c7ia13wo53olvvFmg1v/BYvra8bxeDNe/966AlJLX8ZM9YSt4fXT3YdaicN3r1Xlsuuh2cko5hqMa8Dat14UOva6lPrqHXUuo3Bh1xCCBkLu1OvQVpftqzj7fLbBdguxtaEeYtvmW7rtxqWWjlr5Cy1XfVo/F7f1qhG44/MaZnntenl27zsVCIzEN4Jj5dnY6sZ7OgBIerTu5ei62IR+PdSFoONJbsfycGAWketo9btj8HGQq3bKPy4MiFzRsGPtYyNZypaTysiQ2HNT7So2j4jUxUZseh3RBR3GW/0uqUdm5YZ0NZ0bye+xDOanrHLxuSNuaXP2vjLNK+PbA6iMtk95+XXTHLLvA69hi1Q68aHWketo9bth1o3GjzJJWTu6Io/FUTkMhG5Q0SOisjVSbnvExEVkYtXGs9BwFs07Xx6C2o0554By3aG7eIdxZItxrVF2ovF68ur37KL770eOg8epRHzdvlbjJEdX8v4PTMdnQzY1zVsfduONYneHNq5KB8cbF/ZfZ3F2PIAswob1jriQK3zX0fly5i8MtS6HGrd3lhmrnU8ySVk7mxQ0ETkEIC3AngugGMAbhGRw6p6uyn3OAA/AuDmzUWzZZSLJ+Av7pYhJskupN7JRHaqEbXdusM+JN+bi9a+y/JlG9Y01U4QbB+1E5eSmmH3Xrf0vczJQFbWmwvvpMbGHN0nSPJaHwyihwuv31XZIfN2oKDW7Y+XWuf/jtou26XW1aHWAWg4yRWR60TkXhH5TJF2pojcJCJ3dr/P6NJFRN7Sndh8WkSetsngCSFbzyUAjqrqXar6VQDvBXC5U+6nAbwZwJfHDK5kFK3LDEKZ7u0O2zRb3jMcnhlo2UGOdvI9E9N6alDGFO32Z6ca0VhsG+X7Wp0s9mgOonJZ2dLYZGZmyFzWxlLro2YMbR0v9ug6ev1F8WT3tTqvyVrYuN5R66h11Lq4DrVuFFo+rvxOAJeZtKsBfERVLwTwke49ALwAwIXdz6sB/NJ6wiSEbArR1X4AnCUiR4qfVxfNnwPg7uL9sS7tVP8Lw3Seqv7bTY+1wjsxptb1hsamwUm3eWLe23p2lz9aUBX762cxRDvfZZ7tz0uzsUUxlHgxewa5tT07hto8AfvbiQy5V65s3xt/mR61YetEBs3Li9rw5qws650YRebPu+5evaic9wDgjWdJE7gGrdsl3omx9I5aR62j1lHrJqD6kKuqvw7gAZN8OYDru9fXA3hxkf4uXfAJAE8QkbPXFCshZBNki0DLD3Cfql5c/Fzb2rWIPArAPwLw4+sazrKMonXlIubtvttyMGW9BbBsp7YTX7a/TFy2bJQfxZ+dbNj69kTA9hvtzreMwTM12WlIbX6za2Rfe+P3Tmq8WLzftn+bZ42VN69Zf2U72bUpxxHdSza+Mq7IdHu03ueW1bVuZ9i43lHrqHXUOmrdxCz7xVNPUtV7utd/COBJ3evqqU2PiLy6P/n5M3xlyTAIISuxqhDWxfA4gPOK9+d2aT2PA/DtAD4uIr8P4JkADm/Rl0+tVese+ZMvLW/OyjRNypT5fRlrBLy8aDEuy5eLemRCvfit+bBlohgstl87FttXnz7EKKj5Hc1dZH49w1Smt5hCb/xlnh1/7d9hZupqhqs2D15f3r3k3U/2/vOuZ9ln63g9Nq91u8BKeketM+1Q6/yy1Dpq3Uis/O3KqrrUlKjqtf3Jz6Nx+qphEEKWZMMfa7kFwIUicoGIPAbAFQAO95mq+pCqnqWq56vq+QA+AeBFqnpkQ8NdmnVo3aGv+7rYTAHxrm7NcHhmJTNPXptevrc7nZkGi7fwZ2UiQ1vrBya/pcwqeaX5rVGLx2sjM5RRXu16276iua6dPESGN+vb3gdlWs1ElvdhNgct4WxW63aKZfSOWlfUo9bth1qX16fWrZ1lH3I/339Upft9b5deO7UhhMwIVX0YwGsB3AjgswDer6q3icg1IvKiaaNrYv1a5xmBaJfbppWLbnZKELXTk5nHMr+2KLcarJbFvSwflYnmxWOpLYmibouZynbgW2kxMVkcmdn35ivrNzLY1vBGKPbPgZ1Le62HmvQdMmBbyHr1jlq3vwy1rh5PVJdaRway7EPuYQBXdq+vBPChIv2V3TfxPRPAQ8VHXwgh24iu+FNrXvUGVX2Kqv55Vf3ZLu0NqnrYKfvsLTvFXa/WeQamZbe6ZozKMuX7zEjZ994C27qLnN0H5cJfiyczcWU5r+2sXZtWu28zU1Iz06salSHlyzn15teb15aHgchA1u4/ewrhtV27p2q60nqqlLW9Ia3LiL652Cl3ZVfmThG5skv7WhH5tyLyX0TkNhF5U1H+h0Tkj0TkU93Pq1YIc316R63L46HWUeuodRvXupY/IfQeAP8ZwH8nIsdE5CoAbwLwXBG5E8BzuvcAcAOAuwAcBfA2AH+n1j4hZFr4sZYFo2idXZyjBda+7xfTKD9b9G35sk5EdHIRLYbZYl4u+K39l/2V8UR1W02rt4APiSeKL3sfjdvrr+XEqNaP7cOmWXMXxWHft8Sfjallfj3j6PXX+kBim59W66JvLj4Vn8iZAN4I4BlY/Pm1NxYG8f9S1b8A4DsBfJeIvKCo+j5VfWr38/aWYDaud9Q6ah21LoZaN4rWnVYroKovD7IudcoqgNfU2iSEbBGrC9pOMIrWWRNkd6ajBdLmR+0MicFr374u823d0mxYMxEZsayPvl7UfzR+j9p8iPO6ZdfdG4NN9/qx/dWuWWZyvPHXfkdGMYrRI3tgiep7/Qy5X6P6yzKt1l0O4Nnd6+sBfBzAPzBlng/gJlV9AABE5CYAl6nqewB8DABU9asi8kksPjK8NBvXO2pd3Edfj1pHravVXxZqHYA1fPEUIeQA4+1WD/0hy+GZgf61t8sL7F80renyFvisT6/NzMC07lBnMXh50Y55bfxeH9Hrvr1+riODHc2h7TcyQvZ1y/izWG1fkWmK+vGMacscZPeSbS8axzL3Qknr9W9hPVqX/U3wGtE3F5e0/F3xJwD4XixOSHq+T0Q+LSIfEJHy/85uB9Q6ah21Ls+j1m1E66onuYQQQtaEXSBbd73LvLJM7USiZSe9dXEtzWlUPzONZdtRXNlpRTaOzNBm/UTXwJvjiMiIZ7Fn4/HMfS2O7Ppn5SMjlt1LWfnswSQr5/WVvZ6O+1Q1/PNmIvJhAN/sZL2+fKOqKjL8Q4EichqA9wB4i6re1SX/awDvUdWviMjfxuLk5HuGtr12qHXUuigtaoNaV68zHjuhdXzIJWTG2E1KMiKemcsW71bzZtsaeoFb+6kZnojSyETjtuVqcWXmrsVgePVq5WpEBr3sp+X6tvSdjaP1nvKuxTrvpSEPINGJzAqMoXWq+pywf5HPi8jZqnqP+ebikuMAnl28PxeLj/r1XAvgTlX9x0Wf9xf5bwfwc8Mj3zDUOmodtS4vS63biNbx48qEzJ3VP9ZCWonMSY83nzat3GUuf9uTkxai+rUYbF7LiuoZjhYjkM1JFvfQOSjb8+5vb85t/cyoeOWs2enTamOqjT8bu9dnmd5yL0X9ZtesZuLK8deu/bK6M63WRd9cXHIjgOeJyBndl7A8r0uDiPwMgMcD+NGyQv8nfzpehMWfapseat0Cat3+MtS6U3nUuo1qHU9yCZk5u/QNyQeecsGzC6C3cHoGIqpTYtuvGbeszSivf1+7v6zhaRlLtju/LGX9bJyt5j2bJ69f77RnyHX3ytfmxM556/x7dVv6LO+J1vb6en2ZFa7xxFr3JgDv777F+HMAXgYAInIxgB9W1Vep6gMi8tMAbunqXNOlnYvFxwD/C4BPiggA/D/dt4v+ve5vjj8M4AEAPzTmoJaGWketo9b59foy1LqVtY4PuYQQMgXebj8a33ttWDIDkC2e2WIeLcrWuES/vb6HGD2v76hMqwEq47K7+pGRrZmz1tOFsrw3vzWTY8efGcPIdNlxZKbXvs/aQ5Dn3Rs2fjt/6zT4E9F91M775uIjAF5VvL8OwHWmzDEEo1bV1wF43VqDXTfUuv3vqXV5WQu17sCwTVrHh1xC5g5Pcqch28XNyvTlvIXSMyO1U4fMiERki7fXRmubXjwtRGP10iKT1XK6UTODNYOMIC+7Hq31o7htenTP2Fha76UyPTN6UX0vnqxMzVhnUOumgVq3v00vnhaodXHcNp1aN3v4kEvI3KEYjk+0QNcW8Khca72+busOcdZ+thuOIK0lrqF1a0Yv6sebw8y4eX14ps+WsWVbx9FaZojpi9Jqhq0lpqGnDp7h9OKxca0CtW58qHVxXEPrUuv29hfFmqVR62YFH3IJmTOKqf/vxrwojY1ncjxDEbXjGbJosaydaGSnFCXeIhzFukxaZkYiovplnC2nD2X9mgFtmYdW4zLEsHksa3RtGZvXEn9mFrMTt1rsqxo8D2rduFDrqHUWap0f87qh1p2ED7mEzB2K4XjUTIBnTjwT4i3A2cLa0q99nRmKVmMa9ZGlt4w/azM6uYgYYjxajfgQWkxf9G+0dl1t3tBTs4jodKTFUEZteHnrNobUuvGg1uV9l3WpdXvf1+pS6+pQ6wDwTwgRQsj42AUoel/bZS+pLZBqXtv3ZXl72lHrr6zn1fH6ivK9tlp3+aO4on5baY25pX4URzbeFlNZG1/Lacky8x2Vr5lIa+xbrjeN28GDWhfnU+uoddS6jcKTXEJmDj/WMgHZzn/rDrNXJvuolpde7kqX6bXTi4iynm2vVi+L2WKNaTa2Wlte21lMQ9sb0ueQ8nb82TWz18T7jSAPGD6XQ+6l7AFj6H1UgVo3AdS6/fWymC3UOmrdElDrFvAhl5C5QzEcH29x7dOB9gWuxUC24hmBktZYh5yK2PZbytm8yBQsaxI8U4wirS/jjb/FJA4xkZHBjdprNdjR75Y+bGxeWu1e8spn+euCWjc+1Dq//ZZyNo9aV2/H5lPrZg0fcgmZOdzxm4jI/NXqoFI+W5AtnmlZZTffLvzZ7rnXvib5cPKGzFtmWr0+szmw42o1hNn1i06FVj1lyNJqRr6172j8Q07fhpyKLalZ1LqJoNbtb59aR62j1m0cPuQSQsjYtOxOtxiP8rW3s1zr20tvMRPRjnQ2lhbT2rI7H5m4lvp9+czkecbFMz2RQcnmpWZoa7HVzFTUb0t6eT/17+19FTHkHo5OjbL7dt2nHGQ8qHXDYrPtUOv21/Wg1hEHfvEUIXNG1/BDhmPnze6Se2UiImNiy3jte/229leS7aJ7JwleLEP67+t4JwFem9GOeq2PqE+vjEc2Vi++rE1vHqNTFq9uLe7IQLeUb6U0k2Ub2fjtQ04WUwa1bhqoddQ6ah21biJ4kkvI3NkhQTsQlIt1tKOrpmyNIScS5WsvnqyNsj9rurJ+ywW7xaja+l6aVz9q06tTi9crv6yBtHPQWt/rKxp/VK6VFvPZ2m523fo+Wu47a5i9B5ghUOvGhVpHraPWUesmhA+5hMwYAfh/N8bEW+y8Ra+2m14zaNmCGpmZ2s6yt+hmZsOanTJ/qNGrxVnWHWp4vLha5te+rtF6/fp2bdlabNkce9eoZs6isi1zVLvvopOUqK3WB4VKM9S6EaHWUetq7fft2rLUuvY6QTPUugV8yCVk7lAMpyNa0DKD1LpwZ7vftQXeOw3x7hOvjj3JyKgt6EPGXzMimUmL5jEyPsucVHgGvuXkosWMlbSaVs+IZSZ+KLU59PrJxrSMofdiItNAraPWUeuodSPD/5NLCCFj4S1gdjGqLXDZYhwZJZtnjZxtMzONLb8zg5XtbA9ZmDMz6rUp5n2tba+tMt+WaY09G/8QI1eOv6Vva7Bq/Qw1WbZN776wJ2O2n8zwDp1nMi3UOmodtW5veWrd6PAkl5CZI0olHZXajny0k9uyC9znlW2tYq7KdsrfXkwZ0YmHZxYiY+Tt/FszUCtvy9RitXHZaxOleX1lpr/lunr53nhbTj1sf5HhsmOJ7kuvTY+yjeghJTPd2UlbA9S6kaHWxbFQ66h11LqNw4dcQuaMgruFU1AzbtnCnV2vFiMWtZMZy8zMtLQRGYshprRl5zvKb2EZQ1IzWF6cNXOTEV2fzIiW1OYt68OWL+crM4gtfbay7Lz1cVDrxodaF9ePoNZR66h1a4EPuYTMHH5BwYh4i0+5OzzEuJX1a7vuXpnI2Hmx1Ha/bTltKOvRuqB7sUvltU2z8XkGpiyzzM5+zRhGebX7wLse2bUuY/XmInqdxZIZT1suO5VoMYXZg8gAqHUjQq2j1lHr6u1YqHVrhQ+5hMwdiuH0tJo+i2fOvLRsV7h2/TMzExnC2niy3e8WSnOSGQdv/NFJjWdMa6cFWdw1w+iNPzthsfX7/luucZnujd8zfF5sUXtDDHtmLKOHGC+eZaDWTQ+1bhjUulP9U+vaodYB4BdPEULI+GSGqBWvbHRyUiuXGRNrBKI2I4PgGQqvv1r7ti9r4lrmo+zbm4Myr2b6vH768tH82Rha071y5fi9eYvmp/UhIJsfey2905GyHVsmMo3ZQ0xtTsl2Qq3zy1Dr2qDWkRXgSS4hM4cfaxkRb9FrXcyihdsuqNHuvY2jfF0zULZetmOf7bZ7O+feODwiA1m2k53EWGwsUZ1o5z+KMTJUtfnK6kf5LebOm+eozyFtZ/GVhrRsz5I9cGzA6FHrRoRaF5en1lHrqHWjwIdcQuYOxXAaWhbiMg1BumdUWrALq92VX4WWdiID7I01G3/WZkuM3imGnY/W+WmJOYt16Ny3mLHo3rKva3Vb5mCIifXKeaZ9Xfdk3xYZH2rd3rLUOmodtW4U+JBLyJxR7vhNgmc2ykU8WqAt3o7+kP7L13ZXuqVtb9G2O9xenJnZaTG3Wb9eG15dz+xZamVaTi1a4rdzFhnSWv0y34vXKzvErLXck94JR3Rv2Osw5BRpKNS6aaDWUeuoddS6ieD/ySWEkDGJFs1oUfUW/aie9zurG8VQW3hbzIEUP0NiiMYX9eGVi/po3d3PTjwiMrPllfPKCPaPLTO2mXlufViwBr1Mi+4lby7s+IcYc2/8No8cPKh1OdS6/Wk2XmodWQGe5BIyd7jjNx7WUNhFsbZjXGsXQTnbvjVM3u643Xn24qjFZNtvOYnJ4o/atvU8s+XNdxS3Nxe23ZbTlOgUqaXfWt99WjSWaM69axidLrQY/CH3Zla+Vm7IPZTFQsaBWketo9ZR6yaGD7mEzBgB+LGWMYkWrNoiWjNbQ8xUzTRl/dUW3toOd7SwezvsUTxZfC3zV4ujxdhGpyG1fqM6Q4yobXcV05c9ENSIHg4sUey16zbUZFeg1o0MtW5vWWpdXj/Ko9YNhlp3Cj7kEjJ3lGo4Oi2L5Cq7uGUbfbve7voyu8w17G55NtZ13XrRyUHLKYAXU8u82Hpen1GZsq8yzzPDLdei1qd3uuOVj9JtmWWNdERWfh3/Dk62Ra0bHWodtY5aF8dZ629ZqHUA+JBLyOzhjt+EtJiJltOIiJad+qxOy/sez+xFC7pnRDLz0LfVYqqinfChu+8lQ08i1EnzymXXxMuLjOmq72vpPbVrm50OtZQv+2mNaQDUugmh1uX5/eu+LWodtW4FqHUL+MVThBAyJgp/V1uL12XZ0iBYM1TW9RY1r02b5+Xb9loWzCEmJtuF79/3MdjxDzEK2fgjvLKtu/WlMfOusRdjS/9l+tp2+pP0KK928mBNr9deS3lB27Um2w21jlpHrcvLU+s2Ck9yCZkzmciTzRDtWtde19LsdawZhKgfz2Bkxifa0c52sWsGyvZfo3Yi4dEylqwfe5pj5631BMYrk50URbSOpzRW0UlJ7V7yTme8+6ZmVrPTkKjsssaXWjc+1DpqHbXOj49aNwp8yCVk5siJqSOYEdmOsre4wkm3RAu4NQ2ZqfEMRmY+MkPiLdC2XhTHEOOS5du2vPHXjFg0/9GJk1c2onYdsusVjX+ZebL9RNen1aBFZjcy/lH/ZTlgb0wrmDdq3YhQ66h1WTlqHbVuBPiQS8jc4Y7feNQWzGghr5kAL69cVFtNnxdTZh48sjK1mGxeNC+tBrI2/haj45me6FTAa9/mZTHVxp6ViWKw6fa6thq6mjmuXfe+fPRQUPYTxdJ6D0ZQ68aDWketo9ZR6yaGD7mEzBx+QcEEZLu9y7QRLZblQur1a+tYU1IzXOrk1eKNjO0QbIxluheX14833944vL68flvNYGbsvX5q5rf1BMhrOyrvxT7UlEX9tsboxYTgfSPUugmg1lHrqHXUuongF08RQshUREbC4i2A0Y6xZ3RaDWXNOGbGJuu/FoM1YVE/XnuRgYuMSm0uPBPV/7bGdUi7Q6mZSM9AemVtuWwu+9fRtfTGHMURXZMWooeibMxbjoicKSI3icid3e8zgnJXdmXuFJEri/SPi8gdIvKp7uebuvTTReR9InJURG4WkfNHGtIwqHWn8qh1eRzUOmrdmrSOD7mEzBkFoLraD2mn37XtseZEg3IZ0S51nxftRJe/o1ht+7ZO7USj3KmujSnbzc6wO+Et44+MjRdjZoSikwVbtkZkzlrabJlf776zbXt9tRislpMHe19796mNq0yz419GdqbXuqsBfERVLwTwke79HkTkTABvBPAMAJcAeKMxiK9Q1ad2P/d2aVcBeFBVnwzgFwC8edVA1wK1Lu+PWrf3d0ub1Lo2qHUn4UMuITNHdLUfMpDMjGU79n1+tsDbBd0zMd5vYP+i6i24to53KhPtapfmKTMoXht9XmREIiPTEpt97c2x17c9QbGxZzFl5iiKuzb+2n1lr1U03ih+71qX5bL7pow/GkM2J3b8S55uTKx1lwO4vnt9PYAXO2WeD+AmVX1AVR8EcBOAywa0+wEAl4rIdpz/UOuoddQ6P84ojVrX2m6T1vEhl5C5oyv+kHai3ftsUY8MSLm4qvM+668WY4tBbMWajWjxt+bQq58ZD6/f2vijubWxWOPREkdmlqPxZzFmxssz6J6Jrj0Y2HbsfVWLz4ut1axl91j272MIq2vdWSJypPh59YDen6Sq93Sv/xDAk5wy5wC4u3h/rEvr+efdx/d+ojB3J+uo6sMAHgLwxAFxbQZqHbWOWudDrQNG0jp+8RQhhIxNtBiXWKMR7cxnu8tlXc/ARDFlRHHY/r34svEsk+7leQbMG3dWrxxLzbgMTff6tjEOaTMzmLWxR3VtG2V+rU6W7o3Pu2+ie6gvP/ThY33cp6oXR5ki8mEA3+xkvb58o6oqMvi85BWqelxEHgfggwB+EMC7BrYxPtQ6ah21bm971Loaa9M6PuQSMmME4EeOp6BmQoD9i17NFLUuwqsaq8jA1cxsZrKyGGvjHxKX10erEa6Z2SyerO3M4HgmaYhpj+It01tOaqL+h2rHMvdN2Y9g/5wM6HrTWqeqzwn7F/m8iJytqveIyNkA7nWKHQfw7OL9uQA+3rV9vPv9RRH5F1j8P7Z3dXXOA3BMRE4D8HgA968+mjVBrcv7ptZR62y/fT61bmWt48eVCZkzq345Ab94ahjeTnGLoQH8BTKa/nKH2NaJypfxeHnle+8kpfxt28jyamOJ5sy2Jdg/n9n4vXnNYorKeWPzDJKdp6gNr74t641/2X+K3r3i5QP5vHjxRPFF1yJqV1C/T2pMr3WHAVzZvb4SwIecMjcCeJ6InNF9CcvzANwoIqeJyFkAICKPBvBXAXzGafclAD6qugXCTK2j1lHrqHUTax1PcgmZOTzJHZGhu+h9vmL/4pztfGemxuurrGv7az2pyE4TbFve62yHPorHpnsmNcLLs/NWM+Q1U+T1lc1FFpuXV3toyOpGeZ7RLsvYOfeMmXcv2QeGIXF5fS7BxFr3JgDvF5GrAHwOwMsAQEQuBvDDqvoqVX1ARH4awC1dnWu6tK/DwgA+GsAhAB8G8LauzDsA/LKIHAXwAIArxhtSArXOj5VaV4/Ny6PWDYJat4APuYTMHT7kToe32LYYJyBfBGtmMms3Sx+CZ2Brr716Nh47V3ZnPDJ0LemRMY1i8NqJrkl0jW2frX2LSW+95lFctk0bn20zMuTZtauNNUprebhoYUKtU9X7AVzqpB8B8Kri/XUArjNlvgTg6UG7Xwbw0rUGuwmoddQ6al09jVq3Vq3jx5UJIWRqsoXMLv62jsBf+FvatwthZsiGUi7e/U8tBs/0eObHMwZDx2/7teU9o2bn3CsfjTOKcZW+h5rz6CHBM2te3T427xQoM+D2XsgMXMvYuDF3cKHWnXpNraPWAdS6DcKTXEJmDj+uPCG1U43s2mSnHXaXOGJVc9eCt+te5mVmMzol8IyDtyvuUZrF7ETEazdLj2Iq0zyT543Blo/GOfT6Re1Eed5pTmn+auPP+veoja3lemXdU+umg1pHraPW+TFR6zYGH3IJmTMK4ATVcDI8o5OZl5a8KL/Wh2ckazG2xGyJTJBXf8gYs3pR/9nYrHmO+svGPNQERnFE9b22MgMa5VuyeV7GLEcPIy1Gr/W61qDWTQu1Lm+bWketo9atHT7kEjJ3qIXj0WJ2ovzs1CIzCmWaZRmDNMSMti70LWbS/m4xurXTAZsmpk407sywtPaZGZ3oOkfjbzFELWbdxp5dg8zMeX1n17gW17qg1o0HtY5a57UXxUGtWy/UOgB8yCVk9vBjLSNiF0Bv8bQLu11wS9R5bc2Ll+/V98xV6+Lbal4zQxPFkI3f6zszl7avMt2LJ+qjxTBHabXTEM/YlXmWFqMZ3Vten72ps33aaxDFE93bWQw2FhvX0PsxgFo3ItQ6ah21jlo3MfziKUIIGZPMVHiLL5x8OO+j115ZzwRF/ampE/328AyYrWNfZ8YpGletflTXS48MTdlH7X00Pi+tZR7L65Vd1zIteijwsIbLu5e8e9RrM7p3o/c2xlqc5OBAraPWUev2Q60bDZ7kEjJ3Vv/D32QZolMOwF+4a4axTyvxdoS9vmwb1ghFv20b3ilEZuSi19mJRf++Vr6McYhhsNcli9Gjpa/a/Nr5y+6Tlj5KWk+ivHvJe5+d/GR9RWPM5mZVqHXTQK2LX1Pr9uZT69YDtQ4AT3IJmT2iq/1U2xe5TETuEJGjInK1k/9jInK7iHxaRD4iIn9uE+PcOryFs7YA9tRMR9ROVMbmRzvunhnty5eGIWrbiyeitkOezU82ftu2mh/bbjT3EuRbWoywV6ecWxtzFFeLGfVMnVe2dp948+W1lcWUtb0BNq11JIBal0Oto9atGWrdAj7kEjJn7KK3zE+CiBwC8FYALwBwEYCXi8hFpthvArhYVf8HAB8A8HOrD2yLsXMXLe7ea2s6PGNW2032rpmafDG/e4YYCFs+6q+lbhZ7S3+2rjhptfpZXFmf3jWJxm+vQ/RQEMVVjiu6lwC/HS/Glv6icZRttJqmFjO5LBvWOuJAraPWUet8qHWjwIdcQsgmuQTAUVW9S1W/CuC9AC4vC6jqx1T1T7u3nwBw7sgxjkt2epGdCNgF0Vsgy5OGqN+yfrnQRmbPq9MSozVrnvmw5sQzw17bNiavfS++qJw1VZkBrxmdMi0zY974o5OobN6jeW59ALBm08Zny3kxe31n91I2X7aM936HTNhOQ63bW95rm1q3P2ZqXZxHBsH/k0vIjBEAstn/u3EOgLuL98cAPCMpfxWAX9tkQJMyZKrLRd8zJpEZsvWidr36LQY0asvme3l2LC2nJd74vfTWHXEvxqjNaH6ivlvK1spneVmsGd54suuf1fPyW8Zj+66V78nmbwAjaB0podZR62rlqXV7odatHT7kEjJ3TqzcwlkicqR4f62qXju0ERH5GwAuBvBXVo5oW4kWLHsq4Jkk244tV1vEo/J2F9q25eXZ9iNj0Of15b047Y66fd1iQr1TlZYykdGJ3mcnBy2nJjaO8n1tDlti9U4Uspijfr17KTtdiky7dx/X2inrlm2sg9W1jrRCraPWUeuodRPDh1xCZs4advzuU9WLg7zjAM4r3p/bpe2NQeQ5AF4P4K+o6ldWDWhr8cxVZPRqC2WL+ct2lDMjEeW1xBLlefUzs5OVqZljL17PyGR9eXOaxR6ZRM9Ye21k1y+KyYs3Gk9EdA9m46vNjW2/JZbsPvZiWAKebowItY5aR61rjyNqa0modQv4kEvInFGsb+fQ5xYAF4rIBVg83F4B4AfKAiLynQD+GYDLVPXejUYzNdkimC381jTYa1YzexavTMuus9d/ZkRsDLavmlHJTBNMnpeWjScyh1ldr59oDJnR9sa/jCmuxdtiqst2amYyum5Rm1ks3j1dK7MKm9c6UkKt29smtW5/WvQ+iyGKl1q3ty1qHYANffFU7U+GEELmgao+DOC1AG4E8FkA71fV20TkGhF5UVfsHwL4egC/IiKfEpHDE4W7FEvpXbTYZWbQK+eVb1nsh5gEu/CX8XjmplysvXQbSxRjWS8av+27xSBkZti24xkfry3PCNs6Q8bvxZmd1LS057Vp82wfkXkr7wNrqGrXYcj16svYcdHATQK1zvym1lHrqHVbzdpPcos/GfJcLL5k5hYROayqt6+7L0LIqig2/UfDVfUGADeYtDcUr5+z0QA2yMp6Fy3u3k5zKy11s9MRry372mvHmrzMNGS76NYo1XbSS3PQMmdRv5nRsmMoy/TptROjbPxlG1F/3pii8beSXVsba1m+NHy1NmrYMWRxrtLPCFq3y1DrnDaoddS6IVDrRmcTJ7nVPxlCCNke+EfDV2I5vastcNEi7JW1ZWrmx/bh7dBH9bxde2/RzuKvmd3yx6sftZXNXcmy5sTG4PXb8u8hG1efX7v+dvzeqYQ9dfCIjGvZZmTMbJxlftRvNn9eu9GYl9Qdat1KUOuodad+U+v2xxqNoaxn31PrNsomHnK9Pxlyji0kIq8WkSMicuTPsLvfM0PI1qO62s+8qepdqXWP/MmXukSnpWxRLBfhqI5nvFoMZmTIbFoUX21B7g1IZirK39E4slijecgMiGeOy/6zE4oyPlsnu0ZlmjVb0euyL+911G/txCB6kGi9Bl5btYcOr98yvTbm7F5tgVq3CtQ6ah21jlp3oJjsi6e6PzFyLQBcfPHFetORX5kqFEJ2AhG5dXAlBYRfNb9RrNYd+ZEfnzgiQg428qN/n1q3hVDrCFkv1LrV2MRJbtOfDCGEkB2AekcImQPUOkLIgWITD7kn/2SIiDwGiz8ZcqC+LZWQWcGPtawC9Y6QgwK1bhWodYQcFKh1ADbwcWVVfVhE+j8ZcgjAdap627r7IYSsid3Rs9Gh3hFygKDWLQ21jpADBLUOwIb+T673J0MIIduJ7NCu3RRQ7wg5GFDrVoNaR8jBgFq3YBMfVyaEEEIIIYQQQiZhsm9XJoRsCdzxI4TMAWodIWQOUOsA8CGXkHmjAPhV84SQXYdaRwiZA9S6k/Ahl5AZI1D+3w1CyM5DrSOEzAFq3Sn4kEvI3KEYEkLmALWOEDIHqHUA+MVThBBCCCGEEEJ2CJ7kEjJ3uONHCJkD1DpCyByg1gHgQy4h84ZfUEAImQPUOkLIHKDWnYQPuYTMHH5BASFkDlDrCCFzgFq3gP8nl5C5o7raDyGEHAQm1DoROVNEbhKRO7vfZwTlruzK3CkiV3ZpjxORTxU/94nIP+7yfkhE/qjIe9VKgRJCDj7UOgB8yCWEEEII2TRXA/iIql4I4CPd+z2IyJkA3gjgGQAuAfBGETlDVb+oqk/tfwB8DsC/LKq+r8h/+8ZHQgghMVujdXzIJWTWrLjbx5NcQsiBYHKtuxzA9d3r6wG82CnzfAA3qeoDqvoggJsAXFYWEJGnAPgmAP9h1YAIIbsIta6HD7mEzBnF1GJICCGbZz1ad5aIHCl+Xj0ggiep6j3d6z8E8CSnzDkA7i7eH+vSSq7A4jSjFN/vE5FPi8gHROS8ATERQnYNat1J+MVThMwdfgsfIWQOrK5196nqxVGmiHwYwDc7Wa8v36iqisiyO4RXAPjB4v2/BvAeVf2KiPxtLE5OvmfJtgkhuwC1DgAfcgkhhBBCVkZVnxPlicjnReRsVb1HRM4GcK9T7DiAZxfvzwXw8aKNvwjgNFW9tejz/qL82wH83HLRE0JIGwdF6/hxZUJmjqiu9EMIIQeBibXuMIAru9dXAviQU+ZGAM8TkTO6byR9XpfW83IA79kzpoWJ7HkRgM+uGigh5GBDrVvAk1xC5g4fVAkhc2BarXsTgPeLyFVYfGPoywBARC4G8MOq+ipVfUBEfhrALV2da1T1gaKNlwF4oWn374nIiwA8DOABAD+0wTEQQg4C1DoAfMglZN4ogBN8yCWE7DgTa133UbtLnfQjAF5VvL8OwHVBG9/qpL0OwOvWFykh5EBDrTsJH3IJmTX8hmRCyByg1hFC5gC1rof/J5cQQgghhBBCyM7Ak1xC5g53/Aghc4BaRwiZA9Q6AHzIJYRQDAkhc4BaRwiZA9Q6AHzIJWTe8IunCCFzgFpHCJkD1LqT8CGXkFmjgJ6YOghCCNkw1DpCyByg1vXwi6cIIYQQQgghhOwMPMklZO7w/24QQuYAtY4QMgeodQD4kEvIvOH/3SCEzAFqHSFkDlDrTsKHXELmDnf8CCFzgFpHCJkD1DoA/D+5hBBCCCGEEEJ2CJ7kEjJ3uONHCJkD1DpCyByg1gHgQy4hM0cphoSQGUCtI4TMAWpdDx9yCZkzCuAE/54aIWTHodYRQuYAte4kfMglZO5wx48QMgeodYSQOUCtA8AvniKEEEIIIYQQskPwJJeQucMdP0LIHKDWEULmALUOAB9yCZk5yj8aTgiZAdQ6QsgcoNb18CGXkDmjgCq/oIAQsuNQ6wghc4BadxL+n1xCCCGEEEIIITsDT3IJmTv8WAshZA5Q6wghc4BaB4APuYQQfkEBIWQOUOsIIXOAWgeAD7mEzBtV/tFwQsjuQ60jhMwBat1J+JBLyNzhjh8hZA5Q6wghc4BaB4BfPEUIIYQQQgghZIfgSS4hM0f5sRZCyAyg1hFC5gC1bgEfcgmZNcqPtRBCZgC1jhAyB6h1PXzIJWTOKPhV84SQ3YdaRwiZA9S6k/Ahl5C5o/xYCyFkBlDrCCFzgFoHgF88RQghhBBCCCFkh+BJLiEzRgEoP9ZCCNlxqHWEkDlArTsFT3IJmTOqi4+1rPJTQUQuE5E7ROSoiFzt5J8uIu/r8m8WkfM3MVRCyIwZQesyRORMEblJRO7sfp8RlPt3IvIFEfk3Jv2CTh+Pdnr5mC6d+kkIOQW17iR8yCVk5ugJXeknQ0QOAXgrgBcAuAjAy0XkIlPsKgAPquqTAfwCgDdvYJiEkJmzSa1r4GoAH1HVCwF8pHvv8Q8B/KCT/mYAv9Dp5INY6CZA/SSEGKh1C/iQSwjZJJcAOKqqd6nqVwG8F8DlpszlAK7vXn8AwKUiIiPGSAghm6bUuesBvNgrpKofAfDFMq3Tw+/BQh9tfeonIWSb2Bqt4//JJWTubPZb+M4BcHfx/hiAZ0RlVPVhEXkIwBMB3LfJwAghM2Pabxx9kqre073+QwBPGlD3iQC+oKoPd++PYaGbAPWTEGKh1gHYkofcW2+99U9E5I6p4+g4C9uzODAWH8bi8+eGVvgiHrzxw/qBs1bs97EicqR4f62qXrtimzvJrbfeep+IfAnbc89s0/3LWHwYy362UutE5MMAvtmp9/ryjaqqiOz0N8NQ61IYiw9j2Q+1bgW24iEXwB2qevHUQQCAiBxhLPthLD7bFMsyqOplG+7iOIDzivfndmlemWMichqAxwO4f8NxTYKqfuM23TOMxYex+GxTLEMZQeugqs+J8kTk8yJytqreIyJnA7h3QNP3A3iCiJzWnXCUOrqV+kmti2EsPoxlPVDrTsH/k0sI2SS3ALiw+7a8xwC4AsBhU+YwgCu71y8B8FFV3elTDkLI7Ch17koAH2qt2Onhx7DQR1uf+kkI2Sa2Ruv4kEsI2RjdTtxrAdwI4LMA3q+qt4nINSLyoq7YOwA8UUSOAvgxxN/ERwghB5U3AXiuiNwJ4Dnde4jIxSLy9r6QiPwHAL+CxZeqHBOR53dZ/wDAj3U6+UQsdBOgfhJCtout0bpt+bjyNv3/Pcbiw1h8timWrURVbwBwg0l7Q/H6ywBeOnZcE7JN9wxj8WEsPtsUy4FCVe8HcKmTfgTAq4r33x3UvwuLb6u36dusn9t0vzAWH8bis02xHCi2SeuEn2ohhBBCCCGEELIr8OPKhBBCCCGEEEJ2hskfckXkMhG5Q0SOisjo/5dERH5fRH5bRD7Vf122iJwpIjeJyJ3d7zM21Pd1InKviHymSHP7lgVv6ebp0yLytBFi+UkROd7NzadE5IVF3uu6WO4oPke/rljOE5GPicjtInKbiPxIlz763CSxTDI35OBCraPWObFQ68jOQa2j1jmxUOvI+KjqZD8ADgH4XQDfCuAxAH4LwEUjx/D7AM4yaT8H4Oru9dUA3ryhvp8F4GkAPlPrG8ALAfwaAAHwTAA3jxDLTwL4+07Zi7prdTqAC7preGiNsZwN4Gnd68cB+J2uz9HnJollkrnhz8H8odZR64JYqHX82akfah21LoiFWsef0X+mPsm9BMBRVb1LVb8K4L0ALp84JmARw/Xd6+sBvHgTnajqrwN4oLHvywG8Sxd8Aou/I3X2hmOJuBzAe1X1K6r6ewCOwvlP4ivEco+qfrJ7/UUsvpX3HEwwN0ksERudG3JgodZR67xYqHVk16DWUeu8WKh1ZHSmfsg9B8DdxftjyG+0TaAA/r2I3Coir+7SnqSq93Sv/xDAk0aMJ+p7qrl6bfdRkeuKj/eMFouInA/gOwHcjInnxsQCTDw35ECxDfcFtS6HWufHAlDrSDvbcF9Q63KodX4sALVup5j6IXcb+Muq+jQALwDwGhF5VpmpqoqFYI7OlH13/BKAPw/gqQDuAfDzY3YuIl8P4IMAflRV/7jMG3tunFgmnRtCloBaF0Oti2Oh1pGDBrUuhloXx0Kt2zGmfsg9DuC84v25XdpoqOrx7ve9AH4Vi48gfL7/WET3+94RQ4r6Hn2uVPXzqvqIqp4A8Dac+njGxmMRkUdjIT7vVtV/2SVPMjdeLFPODTmQTH5fUOtiqHVxLNQ6MpDJ7wtqXQy1Lo6FWrd7TP2QewuAC0XkAhF5DIArABweq3MR+ToReVz/GsDzAHymi+HKrtiVAD40VkxJ34cBvLL7xrlnAnio+IjHRjD//+GvYTE3fSxXiMjpInIBgAsB/MYa+xUA7wDwWVX9R0XW6HMTxTLV3JADC7VuP9Q6ah3ZPah1+6HWUevIFOjE33yFxTeo/Q4W31b2+pH7/lYsvjHttwDc1vcP4IkAPgLgTgAfBnDmhvp/DxYfifgzLD7jf1XUNxbfMPfWbp5+G8DFI8Tyy11fn8biH/nZRfnXd7HcAeAFa47lL2PxkZVPA/hU9/PCKeYmiWWSueHPwf2h1lHrnFiodfzZuR9qHbXOiYVax5/Rf6S7eIQQQgghhBBCyIFn6o8rE0IIIYQQQggha4MPuYQQQgghhBBCdgY+5BJCCCGEEEII2Rn4kEsIIYQQQgghZGfgQy4hhBBCCCGEkJ2BD7mEEEIIIYQQQnYGPuQSQgghhBBCCNkZ+JBLCCGEEEIIIWRn+P8BJxB1PyMaYGIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "timeloop(10)\n", "init()\n", "plot()\n", "print(dh)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result = None\n", "if 'is_test_run' not in globals():\n", " ani = ps.plot.scalar_field_animation(timeloop, rescale=True, frames=600)\n", " result = ps.jupyter.display_as_html_video(ani)\n", "result" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "assert np.isfinite(dh.max('phi'))\n", "assert np.isfinite(dh.max('T'))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.5" } }, "nbformat": 4, "nbformat_minor": 2 }