{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pystencils as ps\n", "import pystencils.plot as psplot\n", "import numpy as np\n", "import sympy as sp" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial 05: Phase-field simulation of spinodal decomposition\n", "\n", "In this series of demos, we show how to implement simple phase field models using finite differences.\n", "We implement examples from the book **Programming Phase-Field Modelling** by S. Bulent Biner. \n", "Specifically, the model for spinodal decomposition implemented in this notebook can be found in Section 4.4 of the book.\n", "\n", "First we create a DataHandling instance, that manages the numpy arrays and their corresponding symbolic *sympy* fields. We create two arrays, one for the concentration $c$ and one for the chemical potential $\\mu$, on a 2D periodic domain." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "language": "python3", "language_info": { "name": "python3", "pygments_lexer": "python3" } }, "outputs": [], "source": [ "dh = ps.create_data_handling(domain_size=(256, 256), periodicity=True)\n", "μ_field = dh.add_array('mu', latex_name='μ')\n", "c_field = dh.add_array('c')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the next cell we build up the free energy density, consisting of a bulk and an interface component.\n", "The bulk free energy is minimal in regions where only either phase 0 or phase 1 is present. Areas of mixture are penalized. The interfacial free energy penalized regions where the gradient of the phase field is large, i.e. it tends to smear out the interface. The strength of these counteracting contributions is balanced by the parameters $A$ for the bulk- and $\\kappa$ for the interface part. The ratio of these parameters determines the interface width." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle {c}_{(0,0)}^{2} A \\left(1 - {c}_{(0,0)}\\right)^{2} + \\frac{κ \\left({\\partial_{0} {c}_{(0,0)}}^{2} + {\\partial_{1} {c}_{(0,0)}}^{2}\\right)}{2}$" ], "text/plain": [ "c_C**2*A*(1 - c_C)**2 + κ*(Diff(c_C, 0, -1)**2 + Diff(c_C, 1, -1)**2)/2" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "κ, A = sp.symbols(\"κ A\")\n", "\n", "c = c_field.center\n", "μ = μ_field.center\n", "\n", "def f(c):\n", " return A * c**2 * (1-c)**2\n", "\n", "bulk_free_energy_density = f(c)\n", "grad_sq = sum(ps.fd.diff(c, i)**2 for i in range(dh.dim))\n", "interfacial_free_energy_density = κ/2 * grad_sq\n", "\n", "free_energy_density = bulk_free_energy_density + interfacial_free_energy_density\n", "free_energy_density" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In case you wonder what the index $C$ of the concentration means, it just indicates that the concentration is a field (array) and the $C$ indices indicates that we use the center value of the field when iterating over it. This gets important when we apply a finite difference discretization on the equation.\n", "\n", "The bulk free energy $c^2 (1-c)^2$ is just the simplest polynomial with minima at $c=0$ and $c=1$. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAGICAYAAAB2hyVtAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbKxJREFUeJzt3Xd4lFX6N/DvtEx6IwnpkITQawwEMEgVsVCkCYKgIrbdtS3ryvpzfdVdWdu6u/YVRSkK2AUEFZAOIfReEkjvdVInU877x2QGsiSQSTJ5pnw/15XLyydzJveEycw959znPjIhhAARERER2YRc6gCIiIiInBmTLSIiIiIbYrJFREREZENMtoiIiIhsiMkWERERkQ0x2SIiIiKyISZbRERERDbEZIuIiIjIhphsEVGHqaurw/Lly3HmzBnLtcrKSixfvhwXL15s8/2WlpZi48aN+OSTT7Bnz56OCJWIqNMopQ6AiDqXwWDAihUrmlzz8PBAZGQkhg8fDrVa3eb7rqysxOLFi/HOO++gb9++AID8/HwsXrwYK1asQHx8vNX3efjwYdx6663o168fevXqBaVSieTk5DbHSETU2ZhsEbkYnU6HxYsXIy4uDmPGjAEA1NTU4MCBA6irq8O7776LmTNnShvkVZYtW4bQ0FDs3r1b6lCIiNqEyRaRixo+fDiWL19u+X+9Xo/Ro0djwYIFmDBhAvz9/aUL7ioXLlxAXFyc1GEQEbUZky0iAgAolUrcfvvt2LdvH9LS0pCYmAgAyMvLw08//YQ77rgD4eHhTcYsX74cgwYNwtChQ63+ednZ2fj5558RGxuLcePGXfP9zMxM/PrrrygoKIBCobAkhtOmTYNKpcJXX32F0aNHIz4+Hvv27cOFCxcwYcIEREZGAgDKy8uxd+9elJWVITo6GsnJyVAqr33Ja+3tmnOjsYcPH8bRo0exaNEi1NTU4Ndff0VDQwNuueUWhIWFtfs+q6ur8dtvv6Gqqgrz5s0DYJq53LFjB4qKijB48GD069cPO3bsQH5+PubOnQsA2LBhA7RabbMzmEeOHMGRI0cwd+5ceHl5ter3QEQ3IIjIpdTV1QkAYt68edd8b/bs2cLNzU2UlJRYrv36668CgPj111+vuT0A8ec//9ny//n5+QKAeOeddyzXzp49KwCIFStWWK4dPHhQhIaGiuTkZFFcXNxsnEePHhWLFi0Svr6+olu3bmLRokVi0aJFIiMjw3Kf7777rrj99tvFxIkTRWJioti8ebMQQojXXntNeHh4iMTERHHPPfeImJgY0bNnT3H27NkmP6O1t2tOa8Y+//zzAoBISUkRCQkJYtasWaJfv37C3d1dbNq0qV33uXv3bjFw4EAxffp00adPHyGEEOnp6aJnz54iKChITJ8+XYwcOVIsWbJEzJs3T0RERFju4+WXXxYymUykpaVdE8PQoUNF7969b/j4iaj1mGwRuRhzsjV8+HDx8ccfi48//lj8+9//FrNmzRJdunQRy5cvb3L7jk62vv76a+Hh4SHmzZsn6uvrbxhvXFycmDp1apNr5vuMi4sTBw4cEEIIYTAYRG5urvjkk08EAPHhhx82ecy33nqr6Nmzp9DpdEII0erbNae1Y82J0fz580VVVZUQQgidTieSk5NFjx49hMFgaPN9zpo1S2g0GiGEEFlZWUIIIRITE0VsbKzIz8+33Mdbb70loqKimiRb+fn5QqVSiSVLljR5XIcOHRIAxBtvvNHiYyci67H1A5GLKi4uxoEDB3DgwAGkpqbi4sWLCA0NRUBAgM1+5j/+8Q/Mnj0bf/7zn7F69ep27XwEgKSkJCQlJQEA5HI5wsPD8eqrryIpKQmPPPKI5Xbu7u546aWXcOHCBWzduhUAWn275lg7dtGiRfD29gZgWq5dsGAB0tLSkJmZ2eb7fOCBB+Dj4wMAiIqKwqFDh3Do0CEsWbIEoaGhlts9+eST18QfGhqK6dOnY8WKFaivr7dcf//996FSqbBgwYIWHzsRWY81W0Qu6n8L5AHg+eefx8yZM7F169Zm66ja429/+xvS09OxatUqzJ8/v0Puc/DgwU3+v6ysDOnp6ejWrRs+++wzCCEAAEIIVFVVAQDOnDmDYcOGtep2kyZNuuZntvZnXD120KBBTe7DXFeWnZ2NmJiYNt3n/z72kydPNntdoVCgf//+OHHiRJPrjz/+ONatW4f169djwYIFqKiowNq1a3HXXXchJCTkmsdNRG3HZIuILB5//HG8+uqr+PTTTy3JlkKhAAAYjcYmt716RqQ1+vXrh6ysLGzYsAH33HMPVCpVu+MNDg5u8v9arRaAqci8ueanixYtQq9evVp9u+a0Zax5BsrM/NgbGhrafJ//+9j1ej0AwM3N7ZrxzV275ZZbMGDAALz//vtYsGABVqxYgdraWixatOia2xJR+zDZIiILd3d3AKaO7WbmJani4uImt7W2I/zdd9+NRx55BDNnzsTkyZPx7bffwtPTs50RNxUSEoKAgADEx8dfM2t3NYPB0KrbtedndPZ9mttjXLx4ETfddFOT76WlpTU75rHHHsPjjz+OI0eO4MMPP0R4eHizs3lE1D6s2SIii02bNgEARowYYbkWGxuLwMBAbN68ucltP/zwQ6vv/4477sAvv/yCAwcOYOLEiaioqGhXvP9LoVDg4Ycfxvfff49Dhw5d8/3s7GxUV1e3+nbt+Rm2iPt6Ro0ahYiICLz33nuWWS4A2Lp1a5PasKvdd9998PX1xf33348LFy7g/vvvt8xkElHH4cwWkYtKT0+3zKJotVqcOnUKK1aswNixY/HMM89YbqdWq/HCCy/g6aefhkqlwoABA7B9+3bMmzcP77//vtU/Nzk5GTt27MBtt92GMWPG4JdffunQGqFXXnkFmZmZSE5Oxvz58zFgwABUVlbi+PHjOH78OPbu3Qtvb+9W3649P8MWcbdEpVJh5cqVuOuuuzB69GjcfffdyM/Px8WLF3HXXXc1uzzp7e2NBQsW4N1334VMJsODDz5oVcxE1Dqc2SJyMUqlEosWLUK/fv0suxFPnTqFyMhI/PTTT9i+ffs1b+pPPfUUNm/eDD8/PxQUFODFF1/E3LlzsWjRIgwbNsxyO09PT8t9m/n7+2PRokXo2bOn5drgwYOxZ88eDBs2rElBeHNmz56NW2+9tcm15u7TTKVS4csvv8SuXbvQrVs3nD9/HkIILFy4EOfPn0fXrl2tul1zWjs2MTERixYtglze9KU2IiICixYtQkRERIfdJwCMGzcOp0+fxm233YaMjAz07dsXX3/9NTQaTYuJ2sKFCwEAo0ePZqd+IhuRieu9yhERkcOLj49Hnz598OOPP17zvY8++giPPvpoh+4SJaKmOLNFROQkiouLLTsczb766iukpaVh1qxZzY5Zs2YNgoOD7erwcSJnw5otIiInkZGRgYULF+Kuu+5CREQETpw4gc8//xwzZsxoMmslhMAnn3yCI0eOYPfu3fjvf/9r2YlKRB2Py4hERE4kIyMD33//PdLT0+Hl5YWxY8fitttua3Ibo9GIhx9+GJ6enhg/fjymTp0qUbREroHJFhEREZENsWaLiIiIyIaYbBERERHZkNMUyBuNRuTl5cHHxwcymUzqcIiIiMiJmQ+KDw8Pb7bv3dWcJtnKy8tDVFSU1GEQERGRC8nOzkZkZOR1b+M0yZaPjw8A04P29fWVOBoiIiJyZhqNBlFRUZb843qcJtkyLx36+voy2SIiIqJO0ZrSJRbIExEREdkQky0iIiIiG2KyRURERGRDTLaIiIiIbIjJFhEREZENMdkiIiIisiEmW0REREQ2xGSLiIiIyIaYbBERERHZkNN0kCciupF6nQGaOl2rbuvnqYJaqbBxRETkCphsEZFT0xuM2J1Wgu+P5uLn0wWo1xlbNc7LTYFJ/cMwPSECw2O7QCG/8ZEcRETNYbJFRE5HCIHTeRp8eyQXPx7PQ0m11vI9uezGZ5kJIVDTYMA3R3LwzZEchPq6Y+rgcNydEIHeoTx7lYisIxNCCKmD6AgajQZ+fn6orKzkQdRELiq3og7fH83F90dzcbGo2nI90MsNUwaF4+4hERgY6deqZOtwZjm+PZqLjcfzoKnXW77XJ8wX04dEYOrgcIT4utvssRCRfbMm72CyRUQOTQiBH47lYW1qFlIul8H8iuamlOPWvl0xfUgEbukZDJWibfuBtHoDfjtXhG+P5OK380XQGUw/QC4Dbu4RhHuHRWNS/9AbJnBE5FyYbDHZInIJJdVa/Omr4/jtfLHl2vDYQNw9JAK3DwiDr7uqQ39eeU0DNp3Mx3dHc3E4s9xy/c4BYXj17gHw8+zYn0dE9ovJFpMtIqe380Ix/rj+OEqqtXBTyvH4mDjMSoxChL9Hp/z8zNIarE3Nxse7LkFvFAj3c8e/5gzBsJjATvn5RCQtJltMtoicllZvwBtbzmP5nssAgJ5dvfGfuUMkK1w/nl2BJ9ceRUZpLeQy4Hdje+CJ8fFtXrYkIsdgTd5h9avBihUrEBsbC6VSid69e+P777/vkDHZ2dm499574e/vj8DAQPzud79DdXX1tXdGRC4rragKd7+3z5JoLRjRDT/+PlnSHYKDovyx8YlRmHlTJIwCeGd7GmZ/tB9ZpbWSxURE9sWqZOunn37Cww8/jL/97W8oLy/H448/jlmzZiE1NbVdY0pLS3HzzTdDq9Xi1KlTyM/Px+DBg7Ft27a2PzIichpCCHyRkoW73tmDM/kaBHq5YfmCRLw8tT/cVdI3HvVWK/HmrEF4Z+4Q+LgrcTSrAnf8Zze+P5ordWhEZAesWkYcP348/P398c0331iuJSUlIT4+HqtXr27zmCVLluDLL79Eeno63N3btpWay4hEzqm8pgHPfXsCP58uBACMig/CW7MG2W3bhZzyWjy19hgONRbQ3z0kAi9P7QefDi7WJyJp2WQZUQiBAwcOYPTo0U2ujxs3Dvv27WvXmB9++AHTpk1rc6JFRM5pX3oJbv/3bvx8uhAqhQzP39EHnz8wzG4TLQCIDPDE2oeH4+kJPaGQy/Dd0Vzc8Z/dTXYvEpFraXWyVVVVhdraWgQHBze5HhISgsLCwnaNycjIgJ+fHyZOnAg3NzdERETgySefvG7NllarhUajafJFRM5jw/E8zF+eggJNPWKDvfDd4zdj8S2xkDvAsTlKhRxPTojH+keGIzLAA9lldbjno/3Yfq7510oicm7t3i5jNBqtbub3v2OEEHjrrbfw6KOPoqqqChs3bsQPP/yAJ554osX7WLZsGfz8/CxfUVFRbX4MRGRftp8rxNPrjsEogGmDw7HxD8noH+EndVhWu6lbIH56chQm9u0KvVHgsdVHsD+9VOqwiKiTtTrZ8vHxgZeXF4qKippcLy4uRmhoaLvGhIaGYtKkSZg+fTrUajWGDBmCJ598Et9++22L8SxduhSVlZWWr+zs7NY+FCKyY/vSS/Do6iPQGwWmDg7HW7MHw9PNcY9x9XVX4b15CZjQJwRavREPfZ6KY9kVUodFRJ2o1cmWTCbDyJEj8dtvvzW5vm3bNowcOdLy/0ajEQaDwaoxo0aNgtFobHIbg8EAubzl8NRqNXx9fZt8EZFjO5pVjsWfH0KD3ogJfbrizVmDoHCAZcMbUSnkePfeBIyM64KaBgMWfnoQ5wpY+kDkKqxaRlyyZAk2bdqETz/9FMXFxXjjjTdw/PhxPP3005bbvPzyy+jSpYtVY5YsWYKtW7fiiy++gEajwf79+/Hvf/8b8+fP74CHSESO4Gy+BvevSEVNgwE39+iCd+8d4lSNQd1VCny8IBFDov1RWafD/OUHcbmkRuqwiKgTWPVKNnHiRHz++ed444030L17d6xZswbff/89hgwZcuUO5XIolUqrxtx000347rvv8NZbbyE8PBwLFizAAw88gDfeeKMDHiIR2bvLJTW475ODqKzTYUi0P/57X6Jd9M/qaF5qJT67fxh6h/qgpFqL+ctTkFdRJ3VYRGRjPK6HiCSVW1GH2R/uR25FHfqE+WLt4uFOf6BzcZUW93y0H5dKahAb5IV1j4xAsI9a6rCIyAo2Pa6HiKijFFeZZndyK+oQG+yFVYuGOX2iBQDBPmqsfigJEf4euFRSg/s+SUFlrU7qsIjIRphsEZEkKmobcN8nKbhcUoMIfw+sXpSEIG/Xmd0J9/fA6odMj/lcQRXu/+wgarR6qcMiIhtgskVEna5aq8f9K1JxrqAKwT5qrHkoCeH+HlKH1eligryw+qFh8PNQ4WhWBRavPIR6nUHqsIiogzHZaqVvj+TgL9+d5JEbRO2k1Ruw+PNDOJZdAX9PFVYvSkL3IC+pw5JM71BffP7gMHi5KbAvvRS//+II9AbjjQcSUYv2pZXg/74/iY0n8qQOBQCTrVbberYQX6Rk4TibERK1yz82n8P+S6XwVivx+QPD0CvUR+qQJDc4yh/LFw6FWinH1rNF+M/2NKlDInJoqRnlWH0gC7svlEgdCgAmW60WFeAJAMgqq5U4EiLHtf1cIVbszQAA/GfuYAyK8pc0HnsyIq4LXp85EADw7vaLSLnEY32I2sr8Xh0VaB/lCUy2Wiky0JRs5ZQz2SJqiyJNPZZ8dQIA8ODNMRjXu6vEEdmfqYMjMPOmSBgF8NS6Y6iobZA6JCKHlF1uTrY8JY7EhMlWK0U3/oNll7EBIZG1jEaBp9cfQ1lNA/qG+eLPt/eSOiS79dKUfogJ8kJ+ZT3+/M0JOEkrRKJOlVPGZMshRQWYpiKzy2v54kdkpY92XcLetFJ4qBT4z9whUCudrzt8R/FSK/HO3CFQKWT4+XQh1qRkSR0SkUNp0BuRr6kHcKUESGpMtlopIsADMhlQ22BAaQ2n9ola61h2Bd765TwA06xNjxBviSOyf/0j/PDnSb0BAK9sPIMLhVUSR0TkOHIr6iAE4KFSIMjbTepwADDZajW1UoGuPu4AgGwWyRO1SlW9Dk98eRR6o8CdA8MwKzFS6pAcxoM3x2B0z2Bo9Ub84Yuj7L9F1Erm9+jIAA/IZDKJozFhsmUFS91WOeu2iG5ECIH/+/4UsspqEeHvgVfvHmA3L3yOQC6X4c1ZgxDkrcb5wir8fdNZqUMicgjm4vhoO6nXAphsWSWycQspZ7aIbuzbI7n44VgeFHIZ/jN3MPw8nP/Mw44W7KPGW7MHAQBWHcjEL6cLJI6IyP6ZN7LZS3E8wGTLKuZCOyZbRNd3qbgaL/xwCgDw1Ph43NQtUOKIHNfonsF4+JZYAMCz35xAfiVn1omu5+plRHvBZMsKUZZlRCZbRC1p0Bvx5NpjqG0wICkmEI+P7SF1SA5vycReGBDhh4paHZ5aewwGI3dEE7XE3npsAUy2rGJp/8BeW0QtevOX8ziZWwl/TxX+NWcwFHLWabWXm1KO/8wdAk83BVIul+GDHTzOh6gl5pkt1mw5qOgupn+4vIo6HhRL1IydF4rx312XAACvzRiIMD/7mcZ3dDFBXnhlan8AwNtbL+JwZpnEERHZn6p6HcprdQA4s+Wwuvq4w00hh94okF9ZL3U4RHalpFqLP64/BgC4b3g33NYvVNqAnND0hAhMGxwOg1HgiS+PQVOvkzokIrtiXnkK8FTBW62UOJormGxZQS6XIeKqTvJEdMWrP51FSXUDenb1xvN39pE6HKckk8nwyrT+iA70RG5FHf75ywWpQyKyK/ZYrwUw2bKaeXdDDuu2iCxSM8rw7ZFcAKblQ3cVj+OxFR93FV69ewAAYOX+DJzJ00gcEZH9yLazMxHNmGxZKZo7Eoma0BuMeOF7U5uHOUOjMCQ6QOKInF9yfBDuGBAKowBe/PEUz2slapTT2HTcXs5ENGOyZSVztpzFXltEAIDVBzJxrqAKfh4qPNt4nh/Z3v/d2RceKgVSM8rx3dFcqcMhsgtZlpkt+9qcw2TLSmxsSnRFcZUWb/1qqhtaclsvBHrZx6GvriDc3wN/GG/qYfbqT+dYLE+Eq5YRObPl2MzZMs9HJAJe23IOVfV69I/wxb3DoqUOx+U8lByL2CAvlFRr8favLJYn1yaEsCwj2lOPLYDJltXM/4DFVVrU6wwSR0MkncOZZfj6cA4A4OWp/dm8VAJuSjn+35R+AICV+zNxroDF8uS6SqobUKczQCYzzfzaEyZbVvLzUMGnsXdHDovkyUUZjAIvfH8aADA7MRIJLIqXzC09g3F7/1AYjAJ//f40i+XJZZk3roX5usNNaV/pjX1F4wBkMhkiWSRPLm5NSibO5Gvg667En1kUL7n/u8tULH8wowzfH2OxPLkmywHUdraECDDZahOekUiurKRaizd/Pg8A+NNtvdDFWy1xRBTh74Hfj7tSLF/FYnlyQfZ4JqIZk602sPTa4swWuaDXNp+Dpl6PfuG+uDepm9ThUKOHRsUgJsgLxVVa/GvrRanDIep05gkQe9uJCDDZapMoNjYlF3U4sxxfsSjeLqmVCkux/Gf7MlgsTy7nylE99lUcDzDZahPzP2QWlxHJhRiMAi/+aOoUP+umSNzUjUXx9mZ0z2Dc1q+rqVj+BxbLk2vJstOjegAmW21inqLMKavlixm5jC8OZuFUbmNR/O0sirdXL9zVF+4qOQ5eLsOPx/OkDoeoU+gNRuRX1gNgzZbTiGxMtqq0elTWsRCVnF9ptRZvbDkHwNQpPohF8XYrMsATvx9rKpb/+6azLJYnl5BfWQ+DUcBNKUewHb4+tSnZqq+vR3Z2NhoaGjpkTHl5OdLS0pp8Xb58uS2hdQoPNwWCfUz/mNyRSK7g9S3noanXo2+YL+axKN7uLb4lFt27eKKoSot/s1ieXICl7UOAB+R2WEtqdbL14osvIjAwEIMHD0ZQUBDeeeeddo9577330L9/f0yaNMnyNWPGDGtD61Tm9g/stUXO7lh2BdYdygYAvDKtH4viHcDVxfIr9mXgYmGVxBER2VaWnZ6JaGZVsrVq1Sq88cYb2Lp1K0pLS7F69Wo8/fTT+PXXX9s9ZvDgwU1mto4cOdK2R9RJuCORXMXrjcuHMxIicVO3QImjodYa0ysEt/Y1Fcu/+ct5qcMhsinze7E91msBViZbH3zwAWbMmIGRI0cCAKZMmYJRo0bhgw8+6JAxhYWFqKystCYkyZizZ/baIme2L60E+9JL4aaQ45mJPaUOh6z07G29IJMBP58uxImcCqnDIbIZS48tO2z7AFiRbBmNRhw5cgQjRoxocj05ORmHDh1q95iUlBT069cPoaGh6NWrF7Zs2dLa0CRhaWxazpotck5CCLzROCNyb1I0IuzsYFe6sfiuPrh7cAQA4M1fLkgcDZHtWHpsOfoyYlVVFbRaLbp06dLkelBQEEpKSto1pk+fPjh48CBKSkqg0Wgwffp0TJ06FSdPnmwxHq1WC41G0+SrM0UGmo/s4cwWOaft54pwNKsC7io5Hh8bJ3U41EZPTegJpVyGXReKcfBymdThENlEth332AKsSLbkctNNdbqm24gbGhqgUCjaNWbGjBkYOnQoAEClUuHVV19FZGQkVq9e3WI8y5Ytg5+fn+UrKiqqtQ+lQ5iz59zyOhiN7LVFzsVoFJaZkPtHxiDEx13iiKitort4YvZQ0+vjmz+fZ29Acjq1DXqUVJs6HTh8suXj4wM/Pz8UFBQ0uV5QUIDIyMgOGwMAMpkM0dHRyMzMbPE2S5cuRWVlpeUrOzu7tQ+lQ4T5uUMhl6HBYERhVX2n/mwiW/vpVD7O5mvgo1bi0dGxUodD7fSHcT3gppTjYEYZdl1sfiWCyFHlNJbz+Lor4eehkjia5llVID9mzBj8/PPPTa5t3rwZY8aMsfx/WVlZkx5ZrRmj1+ubfF+j0eDEiRPo0aNHi7Go1Wr4+vo2+epMSoXcUsPCXlvkTPQGI/75q2lW66FRsfD3dJM4ImqvMD8P3Dfc1B/trV84u0XOxd6XEAErk63nn38eu3btwl//+lekpqbid7/7HbKzs/HHP/7Rcpv//Oc/GDJkiFVjxo4di08//RRHjx7Fr7/+irvuugsqlQqPP/54BzxE27lyRiLrtsh5fHc0F5eKaxDgqcKDyd2lDoc6yGNj4uDppsCJnEr8fLpQ6nCIOoy999gCrEy2hg4dil9++QWpqam4//77kZubi507dzaZgQoMDERsbKxVY1avXo3Dhw9j8eLFeOWVVzBs2DCcPHkS4eHhHfAQbYftH8jZNOiN+Pc2U8fxx8bEwcfdPqfkyXpB3mo8eHMMAOCfv56HgbWm5CTMq0vRXew32VJaO2D06NEYPXp0i99/4okn8MQTT1g1plu3bnjvvfesDUVybGxKzmZdahZyyusQ4qPGghHdpQ6HOtjiW2Kxcn8GLhRWY8PxPEwbEiF1SETtdqXtg/22p+FB1O1gTrZyWLNFTqCuwYB3tqcBMBVUu6ua32VMjsvPQ4VHRpvaeLy99QJ0BqPEERG1n+VcRGep2aKmzFk0Z7bIGaw6kIGiKi0iAzxwz9BoqcMhG7l/ZHcEebshs7QWXx/OkToconYRQlh2IzpNzRY1ZZ7ZKtDUQ6s3SBwNUdtV1evwwY50AMCT4+PhpuRLg7PyUivx2BhTzex/tl1EvY6vXeS4ymt1qNaaOhpEchnROXXxcoOHSgEhTM1NiRzVp3syUF6rQ2ywF+5mHY/Tm5cUjTA/d+RX1uOLlCypwyFqM/MSYldftV2XPjDZageZTMYzEsnhVdQ2YPnuSwCAZ27tCaWCLwvOzl2lwB/GxQMA3t+RhtoG/Q1GENknez8T0Yyvqu0UxTMSycF9uPMSqrR69AnzxR39w6QOhzrJrMRIdOviiZLqBqzYmyF1OERtYm77YM8NTQEmW+0WyV5b5MCKqurx2T7TiQ9LJvaEXC6TOCLqLCqFHE9NMM1ufbQzHZV1uhuMILI/WQ7QPR5gstVu7LVFjuz939JRrzNiSLQ/xvUOkToc6mRTBkWgZ1dvaOr1lqVkIkeS4wA9tgAmW+1maf/AXlvkYHIr6izF0X+a2AsyGWe1XI1CLsMzt/YCAHy65zJKq7USR0RkHUc4FxFgstVu5uMBOLNFjuaDHWloMBgxMq4LRvYIkjocksht/bpiYKQfahoM+Hj3ZanDIWo1g1Egt4I1Wy7BvAOiolYHTT1rHsgxFFdpsf6QqaHlE+PjJY6GpCSTySw7E9ccyOTrGDmMAk09dAYBlUKGUF93qcO5LiZb7eSlViLQyw0Ai+TJcazYexkNelOtVlJMoNThkMTG9w5BfIg3qrR6rDnAvlvkGMzvuRH+HlDY+eYeJlsdgHVb5Eiq6nVYdSATAPDY6DjWahHkchkebTwz8dO9l9lVnhyCo9RrAUy2OoTlQGrWbZED+CIlC1X1evQI8caEPl2lDofsxJTB4Qj3c0dxlRbfHsmVOhyiGzI3E4+084amAJOtDmFOtrK4jEh2rl5nwPI9piLoR26JZV8tslAp5HhoVCwA4KNd6TAYhcQREV2feWYrmjNbriGKjU3JQXx3NBfFVVqE+blj6mCegUhNzRkWBX9PFTJLa7H5VL7U4RBd15VlRPvusQUw2eoQliN7eD4i2TGDUeC/u0yNKx8aFQs3Jf/8qSlPNyUWjugOAPhwZzqE4OwW2S9HORcRYLLVISyHUZfV8sWJ7NbPpwtwuaQGfh4qzBkaJXU4ZKcWjuwOD5UCp3I12JNWInU4RM2q1xlQqDE14WWBvIsI9/eAXAZo9UYUV7EDM9kfIQQ+2JEOwPRm6qVWShwR2atALzfMGWZKxs3PGSJ7k9O4kuStViLAUyVxNDfGZKsDqBRyhPmZlxJZt0X2Z29aKU7mVsJdJcf9I7tLHQ7ZuYdGxUIpl2FfeimOZ1dIHQ7RNczvtZEBHg7RvobJVgeJZK8tsmMf7jTNUMwZGm1pwkvUkgh/D0wZHA7gynOHyJ7kOFCPLYDJVoe5um6LyJ6cyKnAnrQSKOQyPDQqRupwyEGYm5xuOV2A9OJqiaMhasq8Ic0RiuMBJlsdhr22yF6ZZyamDgp3iOZ/ZB96dvXBhD4hEAL4785LUodD1ERWqbnHlv23fQCYbHWYK+0fmGyR/bhUXI3NpwoAAI80zlQQtdZjY0zPmW+P5qCgsl7iaIiusLR94DKia7nS2JQ1W2Q//rvrEoQAJvQJQa9QH6nDIQdzU7dADOseCJ1B4JM9nN0i++FI5yICTLY6jPkfPL+yDjqDUeJoiIBCTb3ljLtHOatFbWSe3foiJQuVtTqJoyECKut00NTrAVzZnGbvmGx1kGBvNdRKOYwCyK/gdDtJ79M9l9FgMGJo9wAkdg+UOhxyUGN6BaN3qA9qGgxYdSBD6nCILLNaQd5qeLo5Rs9AJlsdRC6XWTJsFsmT1CprdVh9IBPAlZkJoraQyWSWmdEVezNQ12CQOCJydY50JqIZk60OZF5KZJE8SW11SiZqGgzo1dUHY3uFSB0OObi7BoYhMsADpTUN+OpwttThkItzpDMRzZhsdaArRfJMtkg69ToDVuy9DAB4dEysQ3RXJvumVMixeFQsANOmCz3rUklC5o1onNlyUZbGpuXckUjS+epwDkqqGxDh74G7BoZLHQ45idmJUQj0ckNOeR02ncyXOhxyYeaZrWgH2YkIMNnqUOYsmzVbJBWjUWDFHtOs1kOjYqBS8E+cOoaHmwILR3QHAHyy5zKEENIGRC7L/B7LZUQXZe7OncNkiySy82IxLpXUwEetxKzEKKnDISczb3g03JRynMipxJGscqnDIRdkNArkmI/qcfaZrczMTOzevRv5+a2fSm7tmNraWmzduhUnT55sS2iSMv/Dl9Y0oEarlzgackUr9mYAAGYPjYK32jG2RJPjCPJWY+og09L0p43PNaLOVFytRYPeCIVchjA/d6nDaTWrki2DwYCFCxeib9++eOqppxAbG4s//vGPHTrm8ccfx6RJk/DSSy9ZE5pd8PNQwc9DBQCWzJuos6QVVWHXhWLIZLAs9xB1tAduNh1mvuVUAfIq+DpHncu8AS3c3x1KByqTsCrSd999Fz/++COOHz+Ow4cPY/fu3Xj33Xexfv36DhmzZs0anDlzBhMmTLD+kdgJ1m2RVMyzWhP6dEV0F8eZXifH0jfcF8NjA2EwCqzcnyl1OORiHLFeC7Ay2frss88we/Zs9OjRAwCQmJiIiRMn4rPPPmv3mLS0NCxZsgRr1qyBUum4yx9s/0BSqKzVWY7mebBx5oHIVsyzW18ezGKTU+pUlrYPzpps6fV6nDp1CgkJCU2uJyQk4NixY+0a09DQgDlz5uCVV15BfHx8q+LRarXQaDRNvuxB9yAvAMClkmqJIyFXsjY1C3U6A3qH+mB4LI/mIdua0KcrogI9UFmnw3dHc6UOh1yI+b21W5CTJlvV1dXQ6/UIDGz6Qt6lSxeUlze/K6W1Y/785z8jKioKDz30UKsDX7ZsGfz8/CxfUVH2sfMqPsQbAHCxkMkWdQ69wWhZznnw5hg2MSWbU8hllrrAFXvZBoI6j/m9tWeIj8SRWKfVyZabmxsAoK6uaUFkbW2t5XttGXPw4EF88MEHuOeee7B161Zs3boVJSUlKCoqwtatW1Ff3/yhzkuXLkVlZaXlKzvbPo6QiG98AqQVMdmizvHLmULkVtQh0MsNUwaziSl1jtlDo+DlpsDFomrsSSuROhxyAQajQHqx6b01vqu3xNFYp9XFUZ6enggODkZOTk6T6zk5OejevXubx+j1eiQnJ2P58uWW76enp0Mul+Mf//gH1qxZA3f3a7d3qtVqqNXq1obfaeJCTMuIpTUNKK3Woou3/cVIzsV8NM+8pGi4qxQSR0OuwtddhZk3ReLz/ZlYsTcDo+KDpQ6JnFxOeS20eiPUSrmlr6WjsKpAfuLEifjhhx8s/6/T6bBx40ZMnDjRcu3SpUvYuXNnq8eMHDnSMqNl/kpKSsKoUaOwdetWdO3atc0PTgqebkpEBph2JHJ2i2ztZE4lUjPKoZTLMH94N6nDIRezcGR3AMD2c0W4XFIjbTDk9MxLiHHB3lDIHatcwqpk64UXXsC5c+ewcOFCfPXVV5g1axa0Wm2TvlkrV67E1KlTrRrjbCx1W0y2yMbMs1p3DgxDV1/HafBHziE22BvjeocAAD5rfC4S2Yr5PdXRlhABK5OtXr16ISUlBZ6envjss8/QvXt3pKSkIDQ01HKb2NhYjBkzxqox/2vw4MEYOHCg9Y/GTsR3Zd0W2V5RVT02nMgDcGUrPlFne+Dm7gCArw/nQFOvkzYYcmoXi6oAXJnQcCRWN7Tq3bs3Pvjggxa/v2DBAixYsMCqMf/rb3/7m7Vh2ZUelpmtKokjIWe25kAWdAaBhGh/DI7ylzocclHJPYIQH+KNi0XVWJ+ajYdGxUodEjkp8wRGDwfbiQjwIGqbYPsHsjWt3oA1KaZ2D5zVIinJZDLLc/Dz/RkwGNkGgjqe0SgsyZbTLyNS65hntoqqtKis5bQ6dbwNx/NRUt2AUF93TOrf8pI8UWe4e0gE/D1VyC6rw9azhVKHQ04or7IOtQ0GqBQydAt0rJ2IAJMtm/BxV1lOI08r5lIidSwhhKUwfsHIblA50GGs5Jw83BSYMzQawJVNG0QdyVwcHxvk7VAHUJs5XsQOogeXEslGUjPKcTpPA3eVHHMb3+CIpLZgRDco5DIcuFSGs/n2cXwaOY+0xvfSHg64hAgw2bIZcyd5tn+gjmaeObh7SAQCvJo/vYGos4X7e1iWtDm7RR3NkXciAky2bMZcwMdkizpSdlktfj5dAAC4fyQL48m+PNjYBuL7Y3kordZKGww5FUuPLQfciQgw2bIZc/adVsiaLeo4qw5kwihM2+17hTrmiw45r4ToAAyK9EOD3ogvUrKkDoechBDCsozoiDsRASZbNmOu2cqrrEcVG/1RB6ht0GPtQdMbmLmRJJE9uboNxKoDmWjQGyWOiJxBoUaLKq0eCrkM3bt4SR1OmzDZshF/TzcE+5gOoU4v5plh1H7fHMmFpl6P7l08MbZXiNThEDXrjgFhCPFRo6hKi59O5ksdDjkBc71W9y6ecFM6ZtrimFE7iCvNTbmUSO0jhMDKfRkATIf/yh3sEFZyHW5KueVQ9JX7M6QNhpyCeVe/o9ZrAUy2bMpSt8UieWqnlMtluFhUDU83BWbcFCl1OETXNWdYFFQKGY5kVeBUbqXU4ZCDc+QDqM2YbNlQj65s/0AdY9V+09E804ZEwNddJXE0RNcX4uOOSf3DAACrD2RKHA05urTGZcQeDtr2AWCyZVPxPJCaOkChpt7S7uG+xuUZIntnfq5+fywXlXXcJERtI4TABS4j0vWYk62c8jrUNugljoYc1ZcHs6A3CgztHoA+Yb5Sh0PUKkO7B6B3qA/qdUZ8czhH6nDIQZVUN6CyTge5DIgNdsydiACTLZvq4q1GoJcbhAAucUcitYHOYMSXje0e5nNWixyITCazPGdXH8iE0SgkjogckXllKDrQE+4qhcTRtB2TLRvrwaVEaoetZwpRqNEiyNvNchQKkaOYNiQC3molLpXUYF96qdThkAMybzDr4cBLiACTLZuL54HU1A4rGwvj5wyNhlrpuJ/qyDV5q5WYkRABgG0gqG0uOnjneDMmWzZ2pUieyRZZJ62oCvsvlUIuA+YmRUsdDlGbmJcSt54tRF5FncTRkKNx9AOozZhs2Vh8Y/sH9toia5nbPYzv0xUR/h4SR0PUNvFdfTA8NhBGAUv9IVFrpTn4AdRmTLZszJyNZ5bWoF5nkDgachQ1Wj2+OZILAFgwgoXx5NgWjOgOAPjyYDbPS6RWK6tpQEl1AwAgLsRxdyICTLZsLthHDV93JYwCuFzCHYnUOt8fy0W1Vo+YIC/cHBckdThE7XJr367o6qtGSbUWWxp7xhHdiHlWKzLAA55uSomjaR8mWzYmk8ksS4ms26LWEEJYlhDnD+/GcxDJ4akUcswdZqo7XL2fHeWpdZylXgtgstUpLGck8kBqaoVDmeU4V1AFd5UcMxN4DiI5h7nDoqGQy3AwowznCjRSh0MO4MpORMeu1wKYbHWKHtyRSFYwt3uYOigCfp48B5GcQ1dfd9zWryuAK5s/iK7nSo8tzmxRK3AZkVqrqKoeW07lAwDuY2E8OZn7hncHAHx3NBdV9Twvka6Py4hkFfMTJaOkhjtx6LrWp2ZDZxAYEu2P/hF+UodD1KGGxwaiR4g3ahsM+LZxty1RcyrrdCjUaAFwZotaKczPHV5uCuiNApml3JFIzdMbjFiTYupDdB/PQSQnJJPJLM/tVQcyIQTPS6TmmZcQw/zc4ePu+OUUTLY6gUwmQw8uJdINbDtXhPzKegR6ueGOAWFSh0NkE9MTIuDppkBaUTUOXCqTOhyyU2mNS4jOMKsFMNnqNDwjkW7EXDQ8OzHKoU+3J7oeH3cV7h5iOi9x1YEMaYMhu2XZiejgnePNmGx1kitnJLL9A10rvbgae9JKIJMB83gOIjk58+aPn08XolBTL3E0ZI/Mq0COfgC1GZOtTmJ+wvCMRGrOmgOmWq1xvUIQFegpcTREttU71BfDugfCYBQ8L5GadeVMRCZbZAXzVOil4hroDdyRSFfUNujx1eFsAMB8tnsgF2F+rn+RkgUdXxPpKtVaPXIr6gCwZousFOHvAXeVHA0GI7LKaqUOh+zIj8fyUFWvR3SgJ0bHB0sdDlGnmNQvFEHeahRVafHrmUKpwyE7kt44qxXso4a/p5vE0XQMq5Mto9GI7du3Y8WKFdizZ0+HjamursaWLVuwevVqHDx40Nqw7J5cLmMnebqGEAKrDpgK4+clRfMcRHIZbko55g6LAsCO8tTURSdbQgSsTLZqa2sxZswY3H///diyZQtmzpyJadOmQa/Xt2vMF198gcGDB+O9997D5s2bMWXKFNxyyy2oqnKuYnLzUiLrtsjseE4lTudp4KaUY1ZilNThEHWqOcOiIZcB+y+V8nWRLJypc7yZVcnW66+/jvT0dBw5cgTr1q3D/v37sW3bNnzyySftGhMWFobjx49jw4YNWLNmDU6dOoWjR49ixYoVbX9kdsgys8UDqanR6sZZrTsHhCHQyzmmy4laK8LfA+N6hwAw1W4RAUBaY9uHHk5wALWZVcnW2rVrcc899yAoKAgAEBMTgzvvvBNr165t15ixY8fCy8vL8v9+fn7w8PCAwWCw6sHYu3guI9JVKmt12HA8DwAwfzjbPZBrmtfYUf7rw9moa3Cu13xqG2dcRlS29oY6nQ4XL15E3759m1zv27cvtm/f3u4x5eXlWLduHWpqarBhwwYkJSVh8eLFLcaj1Wqh1Wot/6/RaFr7UCRjPpA6ragaBqOAgvU5Lu3rIznQ6o3oHeqDhOgAqcMhksTo+GBEBnggp7wOG0/kcTndxdU1GJBdbtpE5iw7EQErZraqq6thNBrh7+/f5HpAQECLiY41Y+rq6nDs2DGkpqYiLS0NISEhkMtbDm/ZsmXw8/OzfEVF2f8faFSAB9yUcmj1RuSW10kdDklICIE1KY2F8cO7QSZj4k2uSS6X4d7GRr6ruZTo8tKLqyEEEOCpQhcnKq1odbLl4eEBANcUrWs0Gnh6Nt+E0Zox4eHh+PDDD7F27VqcOHECv/zyC15++eUW41m6dCkqKystX9nZ2a19KJJRKuSIDTItl7KTvGvbf6kUl4pr4OWmsBxdQuSqZidGQaWQ4Xh2BU7lVkodDknoSjNTH6f6ENrqZMvd3R2RkZG4fPlyk+uXL19Gjx49OmwMAAQGBmL06NFISUlp8TZqtRq+vr5NvhxBPA+kJlzpGD91SAS81a1ezSdySkHeakzqbzp83TzjS67JPBHRw0mO6TGzqkB+ypQp+Prrr9HQ0ADANEP1448/YsqUKZbbHDp0qMkuwtaMSU9Pb/JztFotDh06hLi4OOsfkZ3jgdRUVFWPn08XAADmJ7FjPBEAzG9cSvzhWB409TqJoyGpXDmA2oWTrf/7v/9DdXU1Jk6ciGXLlmH8+PEIDg7GE088YbnNxo0b8fTTT1s15v7778f8+fPxxhtv4OWXX8bQoUPR0NCAF154oQMeon0xP4HSuIzostanZkNvFEiI9kffcMeYkSWytWExgYgP8UZtgwHfH82VOhySyNXLiM7EqmQrLCwMR48exeTJk1FQUICFCxciJSWlyRJeYmIiHnzwQavG7Ny5EzNmzEBpaSl0Oh3+8pe/4Pz58+jWzfk+9ZsPpL5YVA0hhMTRUGczHbxrqi+cx1ktIguZTIZ55kL5A5l8fXRBWr0BGaU1AK68VzoLmXCSZ7RGo4Gfnx8qKyvtun5LZzCizwtboDcK7H1uHCL8PaQOiTrRtrOFWPT5Ifh7qnBg6Xi4qxRSh0RkNzT1OiT9fRvqdAZ89egIDO0eKHVI1InOFWgw6V+74eOuxIkXJ9p9gbw1eQcPou5kKoUcMeYdiewk73LWNG5tn5kQyUSL6H/4uqswZVA4gCunK5DruLpey94TLWsx2ZKAeXqUZ4G5luyyWvx2vggALH2FiKip+Y0d5TefLEBptfYGtyZnctFJ67UAJluS6NH4ROKORNeyNjULQgA39+iC2GDnqkcg6igDIv0wMNIPDQYjvjqcI3U41InMG8ecrV4LYLIliStnJHIZ0VU06I1Yl2oqjGe7B6LrM/+NfJGSBaPRKcqKqRXMExDOdEyPGZMtCXBHouv55UwBSqobEOKjxoS+XaUOh8iuTR4UDh93JbLKarE7rUTqcKgT6AxGXC4x70TkMiJ1gJggL8hlQFW9HkVVrElwBeZi3zlDo6BS8M+O6Ho83BSYkRAJgIXyriKztAZ6o4CXmwLhfu5Sh9Ph+KovAbVSge5dzDsSWbfl7NKKqnDgUhnkMmDOMBbGE7XG/OGmv5VtZwuRX1kncTRka1cvITrbTkSAyZZkerBuy2WY2z2M690V4eyrRtQqPUJ8kBQTCKOApREwOS/zTsQeTrgTEWCyJZmejWvSZ/M1EkdCtlTXYMA3jTuqzJ/Uiah1zG0g1h7Mgs5glDgasiXze2FPJ9yJCDDZksygKH8AwLHsCknjINvacCIPmno9ogI9cEt8sNThEDmU2/qFIsjbDUVVWmw7Wyh1OGRD5vdC83ujs2GyJZHBjU+oi0XVqOIJ905rTWNx773DukEud746BCJbclPKMTsxCgCw+kCWxNGQrRRq6pFfWQ+5DBgQ4Sd1ODbBZEsiwT5qRPh7QAjgZE6l1OGQDZzMqcTxnEqoFDLMToyUOhwihzR3WDRkMmBPWomlNQA5l6NZFQBM5TVeaqW0wdgIky0JDY72BwAc5VKiU1qTYprVur1/GLp4qyWOhsgxRQV6YkxP0xL8FylsA+GMzEuIQxrfE50Rky0JDWHdltOqrNPhh2N5AK4U+RJR25j/hr46nIN6nUHiaKijHcsuB3ClvMYZMdmS0OCrki12kncu3x7JQZ3OgF5dfTC0e4DU4RA5tDG9QhDh74GKWh02nsiXOhzqQAajsJTSDI5y3tdKJlsS6h/hB6VchuIqLfIq66UOhzqIEAKrGgvj54/o5pQN+og6k0Iuw7zG1imr2FHeqVwsqkJNgwFebgqnPBPRjMmWhNxVCvQOM/XbOtZYIEiOb396KS4V18DLTYG7h0RIHQ6RU5idGAU3hRzHsytwIqdC6nCog5jf+wZG+kPhxDu2mWxJ7MpSYrm0gVCHMX/ynp4QCW8n3VlD1NmCvNW4fUAoAJ6X6EzMNcuDnbg4HmCyJTnzGjWL5J1DQWU9fjljar7IwniijnVf49/Uj8fzUFnL/oTOwJJsOXFxPMBkS3LmJ9jJ3EoeR+EE1qZmwWAUGNY9EL1CnfOMLyKp3NQtAL1DfVCvM+LrIzlSh0PtVKPV40Kh6XzgIUy2yJZig7zg465Evc6I8wU8lNqR6QxGfHnQ1OV6/gjOahF1NJlMhvsa/7ZWH8iE0chd3I7sRE4ljAII93NHiK+71OHYFJMticnlsiYtIMhxbT1TiEKNFkHebpjUL1TqcIic0rTBEfBWK3G5pAb70kulDofawVXqtQAmW3aByZZzMBfGzxkaDTcl/7SIbMFLrcSMBNMu31UHMqQNhtrFFZqZmvEdwQ4w2XJ8aUVV2JdeCrkMmJsULXU4RE5tXmOh/K9nCpFfWSdxNNRWV4rjnbeZqRmTLTtgTrbSi6uhqecOG0e0+oCpVmt8n66I8PeQOBoi59azqw+SYgJhFMCXKVlSh0NtkF9Zh0KNFgq5DAMi/KQOx+aYbNmBLt5qRAV6QAjgRHal1OGQlWob9PjmsGln1H1s90DUKcyF8l+mZqNBz53cjsbczLRXVx94uCmkDaYTMNmyE1f6bbG5qaP54VgeqrR6dOviieQeQVKHQ+QSJvYNRbCPGsVVWvxypkDqcMhKrlQcDzDZshus23JMQgis2t94DmJSN8id+LgJInvippRj7tAoAOwo74iOukgzUzMmW3bi6mRLCPaOcRRHsytwJl8DtVKOmTdFSh0OkUuZmxQNhVyGA5fKcLGQfQodhd5gxMkcU8mMszczNWOyZSf6hftCpZChpLoBOeXcXeMoVjfOak0eFI4ALzeJoyFyLWF+HhjfOwQAZ7ccyYXCatTpDPBRKxEX7C11OJ2CyZadcFcp0CfMFwCXEh1FWU0DNp7IB8DCeCKpmAvlvzmSixqtXuJoqDXM73EDo/xcpvSCyZYdYd2WY1l/KBsNBiMGRvphkItMhRPZm5vjghAT5IVqrR7fH8uVOhxqBVdqZmrGZMuOMNlyHAajwJqUxsJ4zmoRSUYul2FeYyPhVfszWfPqAFypmamZ0toB1dXVWL9+PTIzMxEfH4/Zs2fDze36tSqtGXP06FHs3LkTOp0OQ4cOxZgxY6wNzeGZk61TuZXQGYxQKZgL26tdF4qRXVYHX3clJg8MlzocIpc266YovPnLeZwrqMLhzHIkdg+UOiRqQVW9DheLqgFwZqtFpaWlSExMxPvvv4/6+nq8+uqrSE5ORm1tbbvGTJ8+HYsXL0ZGRgYKCgowc+ZMzJkzx+U+ocQEecHPQwWt3ohz+dxZY8/MxbizEqNcoiEfkT3z81RZPvSwUN6+ncyphBBAhL8Hgn3UUofTaaxKtv7+97/DYDBg165deO2117B7926kp6fjvffea9eYP/3pTzh06BD+9a9/4a233sL27duxbt06/PLLL21/ZA5IJpNZan/Y3NR+ZZfVYvv5IgCwLF8QkbTMhfI/nSxASbVW4mioJUddrJmpmVXJ1rfffotZs2bB09MTANClSxdMnjwZ3377bbvGjBgxosmY/v37w83NDdnZ2daE5xTM06pHWbdlt744mAUhgFHxQYh1kW3LRPZuYKQ/BkX6ocFgxPpDrvfe4SjM9Vqu0l/LrNXJllarRWZmJnr06NHkeo8ePXDhwoUOGwMA69evh06nw8iRI68bj0ajafLlDIawSN6u1esMWHvQdPDtvCQWxhPZE/NmldX7M6E38LxEeyOEuKo43l/SWDpbq5Mtc42Vr69vk+t+fn6oqanpsDFnz57FY489hqeffhp9+/ZtMZ5ly5bBz8/P8hUVFdXah2LXzMuIl4prUFmrkzYYusYPx3JRXqtDhL8Hbu3bVepwiOgqkweFI9DLDXmV9fjlTKHU4dD/yKusR3GVFkq5DP0j/KQOp1O1Otny8vKCTCZDZWVlk+sVFRXw9m5+KcXaMWlpaZgwYQLuvPNOvPHGG9eNZ+nSpaisrLR8OcuSY6CXG7p1MS25Hs+pkDYYakIIgRV7MwAAC0d2g8JFmvEROQp3lQL3DjPVUX7W+LdK9uNYVgUAoHeYD9xVrrWxqNXJlpubG2JjY3H+/Pkm18+fP48+ffq0e0x6ejrGjh2L0aNH4/PPP4dcfv3Q1Go1fH19m3w5C/bbsk8HLpXhXEEVPFQK3JPIwngiezR/eDco5TIczCjDqdzKGw+gTuOKzUzNrCqQnzVrFtavX2+pj8rPz8fGjRsxa9Ysy222bt2Kl156yaoxly5dwpgxY3DLLbdg1apVUChcK+P9X0y27NOKvZcBANMTIuDnqZI4GiJqTqifO24fEAYAlplosg+u2MzUzKpk67nnnkNAQACSkpLwyCOPYMSIEUhISMAjjzxiuc2ePXvw9ttvWzVm4sSJqKioQHh4OJ5//nk899xzeO6557B169YOeIiO5+pky9V6jdmr7LJa/HrWVAPywM3dpQ2GiK7L/De64XgeiqvYBsIe6AxGnGycaXTFmS2rOsj7+fkhJSUFGzZsQFZWFiZPnozbb7+9yUzUhAkT4O/vb9WYRx99FHr9tQeIuru7t+EhOb6+4b5wU8hRVtOA7LI6RDfWcJF0Pt+XYWn30CPER+pwiOg6EqIDMCjKH8ezK/BFShaenBAvdUgu73xBFep1Rvi4KxEb5CV1OJ3O6uN63NzcMGPGjBa/n5ycjOTkZKvGLFmyxNownJpaqUCfcF8cz67A0exyJlsSq9Hqsa6xb8+DN8dIHA0RtcaDN3fHk2uPYXVKJh4bEwc3JY8/k9LVLR/kLri5iM8+O8V+W/bj2yM5qKrXIybIC6N7BksdDhG1wu39wxDio0ZxlRY/ncyXOhyX56r9tcyYbNkpFsnbB6NRYMW+DADAwhHdXPITGZEjclPKcV9jk9MVey+z/lViTLbILpmfkKfzNGjQsxOyVHZdLMal4hr4qJWYmegcjXOJXMXcpGi4KeQ4nlOJI409nqjzaep1SC+uBsBki+xMty6eCPBUoUFvxNl85ziKyBGZt47PSoyCt9rqEkciklCQtxpTBocDuNK6hTrfiexKCAFEBXqgi7da6nAkwWTLTslkMsvRPVxKlEZaUTV2XiiGTGbqGE9EjsfcBmLzqQLkV9ZJG4yLutLM1PX6a5kx2bJjrNuS1sr9GQCA8b1D0K2L621VJnIG/cL9MCwmEAajwOoDmVKH45JcvV4LYLJl15hsSaeyToevD+cAAB5guwcih/Zg4+zWFylZqNcZpA3GxQghmGyByZZdMz8xL5fUoKK2QdpgXMxXh7JR22BAr64+GBnXRepwiKgdJvTpigh/D5TX6vDDsVypw3EpOeV1KKlugEohQ79w5znD2FpMtuyYv6cbYho77XJ2q/MYjAKfNbZ7uP/m7pDJ2O6ByJEpFXIsGGFuA5HBNhCdyPze1SfMF+4q1z33mMmWneNSYufberYQOeV18PdUYdrgCKnDIaIOMGdoNDxUCpwrqML+S6VSh+MyuIRowmTLzt3UzbR7Y18aXxw6y2eN7R7mDI2Gh5vrfhIjciZ+nipMTzB9eDL/jZPt7U0rAXDlvcxVMdmyc2N6mY6HOZRZxrqtTnA2X4P9l0qhkMssyw5E5BzMbSB+PVuI7LJaaYNxAbkVdThXUAW5DLgl3rWPOmOyZeciAzzRO9QHRgHsvFAsdThOz/yJd1K/UIT7e0gbDBF1qB4hPhgVHwQhgM8b6zLJdrafKwJgmtUK8HKTOBppMdlyAON6hwAAtp0tkjgS51ZW04DvG3cqmT8BE5FzebCxlcu6Q9mo0eoljsa5bT9bCAAY17urxJFIj8mWAxjfx5Rs7ThfBL2B5yTaypcHs6DVGzEgws/l6wuInNXonsGICfJCVb0e3xzJkTocp1XboMfedFOtsfk9zJUx2XIAg6MCEOCpgqZej8OZ5VKH45S0eoOlY/z9I9nugchZyeUyLLyqDYTByDYQtrAvrRQNeiMiAzwQH+ItdTiSY7LlABRyGcb2Mn0yMK+BU8f67kguCjVadPVV465BYVKHQ0Q2NCsxCn4eKlwuqcGWUwVSh+OUtjW+V43vHcIPr2Cy5TDGNU7DbmOy1eEMRoGPdl0CACweFQu1ku0eiJyZl1qJhSO7AwDe35HGJqcdTAiB7eca67X6sF4LYLLlMEbFB0MplyGtqBqZpTVSh+NUNp/Kx+WSGvh7qjB3WLTU4RBRJ3hgZHd4qBQ4nafBroslUofjVE7naVCo0cLTTYGkmECpw7ELTLYchJ+HCkO7m560XErsOEIIvP9bOgBg4Yju8FIrJY6IiDpDgJeb5cPV+7+lSRyNczG/RyX3CHLpI3quxmTLgZhbQDDZ6jg7LxTjTL4Gnm4K3N+4rEBErmHxLTFQKWRIuVzGzUcdyFKvxV2IFky2HIi5buvApVJUsz9Mh3h/h2lWa+6waJdvukfkasL8PHD3ENMRPh/s4OxWRyiu0uJ443mI5o1dxGTLocQGeaF7F0/oDAJ7LrKbfHsdzizDwctlUClkeGhUjNThEJEEHhkdB5kM2Hq2COcLqqQOx+H9dt40qzUw0g8hvu4SR2M/mGw5EJlMZunEy27y7Weu1Zo+JBJhfjyah8gVxQV74/b+oQA4u9URtje+N5nLXsiEyZaDMa+B/3a+CEY242uzcwUabDtXBJkMeGR0rNThEJGEHh/TAwCw4UQ+D6huB63egN2Nqy7jeURPE0y2HMzQ7oHwVitRUt2AE7mVUofjsD5orNW6o38YYoPZ3ZjIlfWP8MOo+KDGnnvpUofjsA5eLkNNgwEhPmr0C/eVOhy7wmTLwbgp5bilZxCAK4d8knWySmux4XgeAOCxMXESR0NE9sA8u7X+UA6Kq7QSR+OYtl21hCiXs2v81ZhsOSBL3RZbQLTJf3enwyiAW3oGo3+En9ThEJEdGB4biCHR/mjQG/Hp3stSh+NwhBDYZu4az3qtazDZckBjegVDJjN16S2orJc6HIdSVFWP9YdyAACPc1aLiBrJZDLL7Nbq/ZnQ1OskjsixpBdXI7usDm5KOW7uESR1OHaHyZYDCvJWY3CUPwA2OLXWp3sy0KA3IiHan8dIEFET43uHoGdXb1Rp9Vi1P1PqcByKeQlxRGwXnsTRDCZbDmq8pZs867Zaq7JOh9UHTC+gj4/pwZPoiagJuVxmqeNcsfcy6nUGiSNyHOwaf31MthyUuW5rT1oJXxBaafWBTFRr9ejV1Yc1BUTUrMkDwxEZ4IGS6gasP5QtdTgOoaK2wXLcEbvGN4/JloPqE+aDMD931OuM2J9eKnU4dq+uwYBP95iKXh8bE8edMkTULKVCjkduMfXe+2jnJegMRokjsn87LxTDYBTo1dUHUYGeUodjl6xOtnJzc/HXv/4VCxcuxN/+9jeUlZV1yJiSkhK8/vrrmDNnDvbs2WNtWC7H1E3e9AliG5cSb2j9oWyU1jQgKtADdw0MkzocIrJjsxKjEOTthtyKOkubGGqZuXZ4HJcQW2RVspWVlYWEhAQcPXoUw4cPx9atW5GYmHjdhKs1Y7744gsMHjwYRUVFWLduHTIyMtr8gFyJeW18+9kiCMFu8i3RGYz4765LAICHb4mDUsEJXSJqmbtKgQeTTeelfrAjnad1XIfeYMSO8+au8Uy2WmLVu87LL7+MsLAwfP/993jsscewefNmaLVavP322+0ak5ycjPT0dLz55pttfyQuaGRcENxVcuRV1uMcD1Bt0Y/H8pBbUYcgbzVm3RQpdThE5ADmD+8GH7USF4uqsZUNpFt0JKsClXU6+HuqMCQ6QOpw7JZVydZPP/2E6dOnQ6FQAAA8PDwwefJkbNq0qV1joqOjoVar2xK/S3NXKXBzXGM3ebaAaJbeYMQHO03HbyxKjoG7SiFxRETkCHzdVbhvRDcAwHu/pXH1oAXmMpaxvUKgYC1si1qdbNXV1SE/Px/R0dFNrnfr1g2XLl3qsDGtpdVqodFomny5IvMa+TZ+8mrWN0dykFZUDX9PFeYNj77xACKiRg8mx8DTTYHjOZXYfKpA6nDs0varjuihlrU62aqvN3Uq9/Zuemivt7e35XsdMaa1li1bBj8/P8tXVFRUu+7PUZmf4EezK1BazfO8rlbXYMA/f70AAPj92B7wdVdJHBEROZIgbzUWjzLtTHx9yznuTPwfWaW1uFhUDYVchlt6Bksdjl1rdbLl7e0NuVx+TTF8aWkp/P39O2xMay1duhSVlZWWr+xs1+yHEubngb5hvhACliJFMvl072UUarSIDPCwLAcQEVlj8S2xCPJWI6O0Fl8ezJI6HLtibqo9tHsA/Dz4YfZ6Wp1sqVQq9OnTB6dOnWpy/eTJkxgwYECHjWkttVoNX1/fJl+uyrIrkXVbFqXVWnyww1Sr9afbekGtZK0WEVnPW63EUxPiAQD/3noRVTwz0cLSNb6xyTa1zKoC+Xnz5mHdunXIyzP1HTl79iw2b96M+fPnW27z9ddfY9GiRVaNofYxLyXuulCMBj2nuQHgne1pqNbqMSDCD5MHhksdDhE5sHuGRiE22AulNQ2WNjKurlqrR8ol06oV+2vdmFXJ1jPPPIPExEQMGjQIt912G0aMGIFZs2bhvvvus9zm1KlT+Oabb6weM2fOHMyZMwcA8N5772HOnDn473//297H5xIGRfqji5cbqrR6HMq4cZNZZ5dRUmM5A3Hp7b3ZLZ6I2kWlkOPZ23oDAD7efQmFmvbVHDuDPRdL0GAwonsXT8QGeUkdjt2z6mhutVqNTZs24fDhw8jKysJbb72F/v37N7nNzJkzMXDgQKvGBAcHY9q0aQBg+S8AxMXFWflwXJNcLsPY3iH4+nAOfj1biJE9gqQOSVJv/HweeqPAmF7BLv+7IKKOcVu/rripWwAOZ5bj7V8v4B8zBt54kBMz9x4b17srZDJ+oL0RmXCS5iEajQZ+fn6orKx0yfqtrWcK8dDKQ/DzUGH/0nHwdLMqj3YaR7PKcff7+yCTAZufHIXeoa73XCAi2zicWYYZH+yHXAb8/NQtiO/qI3VIkqiobcDwZdtQrzNi3cPDkRTbReqQJGFN3sFzS5zE2N4h6NbFE5V1Onx9OEfqcCQhhMCyzecAADMTIploEVGHuqlbICb1C4VRAK9tOSd1OJJZk5KFep0RfcN8MSwmUOpwHAKTLSehkMvw4M2ms7w+3XMZBhc8y2vb2SIcvFwGtVKOZyb2lDocInJCz07qBYVchq1ni3DgUqnU4XQ6rd6Az/ZlAAAeGhXDJcRWYrLlRGYlRsLPQ4WM0lqXO8tLbzDiH42fNBclxyDMz0PiiIjIGcUGe+PeYabTKJb9dNbljvHZcDwfxVVadPVV4y7u9G41JltOxNNNiXlJpheB5btda3vyV4dNx/IEeKrw6BhurCAi23lifDy8Go/x2XQyX+pwOo0QwvLecv/IGLgpmUK0Fn9TTmbhyO5QKWRIzSjHsewKqcPpFLUNesuxPH8YF89jeYjIpoJ91Hj4FtOHute3nHeZ/oZ700pxrqAKnm4Ky+wetQ6TLSfT1dcdkweZpnZdZXZr+e7LKK7SIjrQE/OH81geIrK9h0bFINhHjayyWqxJyZQ6nE7xceN7yuzEKPh58kOtNZhsOaGHkk0Hp24+VYCc8lqJo7GtkmotPtp55VgeTmsTUWfwUivx9ATTRpx3tqdB4+TH+FworMLOC8WQyYAHbu4udTgOh+9MTqhvuC9u7tEFBqPAZ3szpA7Hpv6z7SJqGgwYFOmHOweESR0OEbmQ2YmRiAv2QllNg+VDn7P6ZPdlAMBtfUPRrQs7xluLyZaTMs9urU3NdtqDUy8VV+OLlCwAwHO39+GxPETUqZQKOZ67vQ8AUzlDfmWdxBHZRnGVFt8dywUALL4lRuJoHBOTLSc1umcweoR4o1qrx7rUbKnDsYnXt5iO5RnXOwQj4lyzgzERSWtCnxAM7R4Ard6If/5yQepwbGLVgUw06I0YHOWPhOgAqcNxSEy2nJRcLsOiZNMnkBV7M6A3ONdumS2nCrDldAHkMuDPk3pLHQ4RuSiZTIald5hmt746nIM9F0skjqhj1esMWH3AtAFg8ahYNjFtIyZbTuzuIRHo4uWG3Io6bD5VIHU4Haaoqh5/+e4kAOCR0XHoFeqa55MRkX1IiA7A/OGmVgh/+vo4Kmudp3Tj2yO5KKtpQIS/B27r11XqcBwWky0n5q5SWFohLN99ySk6HQsh8Nw3J1FW04A+Yb6W3UBERFL6yx19EBPkhfzKevz1x1NSh9MhjEaB5XtM7R4eTI6BUsGUoa34m3Ny943oBjelHMdzKnEos1zqcNptbWo2tp8rgptCjn/dM5itHojILni6KfHP2YOgkMvww7E8bDyRJ3VI7bbjQhEuFdfAR63EPUOjpA7HofGdyskFeasxIyECAPDxLsducppZWoNXNp4BYOqpxeVDIrInQ6ID8LvG48Ke/+4UCjX1EkfUPh/vMrV7mJsUDW+1UuJoHBuTLRdgLpT/9WwhMkpqJI6mbQxGgWfWH0dtgwHDYwMtj4mIyJ78YXw8BkT4obJOhz99fcJhyzdO5VZi/6VSKOQy3D+yu9ThODwmWy6gR4gPxvYKhhDAp3svSx1Om3y4Mx2HM8vho1bizVmD2FOLiOySSiHH2/cMglopx64LxZadfI7mkz2m94o7B4Qh3N9D4mgcH5MtF/HQKFOT068O5aCitkHiaKxzOq8S/9pq6l/z4pR+iAzwlDgiIqKW9QjxwXO3m1rS/P2ns7hUXC1xRNbJr6zDhuOmmrPFje8d1D5MtlzEyLgu6BPmizqdAWsau647gnqdAU+vOwadQWBSv1BL/RkRkT1bOKI7bu7RBfU6I55ef9yheh1+vi8TeqNAUkwgBkT6SR2OU2Cy5SJkMhkeaqxz+nxfBhr0jvGH/+bP53GhsBpB3mq8On0AG+oRkUOQy2V4Y+Yg+LgrcTy7Au/95hhnJ9Zo9fgixbT0+RBntToMky0XMnlQOEJ81Ciq0lqmiO3Z/vRSfNJYY/b6zAEI9HKTOCIiotYL9/fA36b1BwD8Z/tFnMipkDagVvjqUDY09XrEBHlhfO8QqcNxGky2XIibUo6FjbtKXv/5HIrseFuypl6HJV8dhxDA3GHRGNebnYuJyPFMGRSOOweGwWAUeHrdMdQ1GKQOqUXZZbX497aLAExNTLkRqeMw2XIx94/sjvgQbxRqtHh09WFo9fb5h///fjyN3Io6dOviif+7s4/U4RARtYlMJsPfp/VHiI8a6cU1eG3LOalDalZtgx4PrzqM8lodBkT4YdZNkVKH5FSYbLkYL7US/12QCB93JY5kVeD//XhG6pCuseVUPr49kgu5DPjn7EHwYjM9InJg/p5ueGPWIADAZ/sysPtiscQRNSWEwLNfn8DZfA2CvN3w0X03wV2lkDosp8JkywXFBHnhP3OHQCYDvjyYhTUp9tMH5lRuJZ771nTI9GNj4nBTt0CJIyIiar/RPYNxX+NZtc+sP460oiqJI7rio12XsPFEPpRyGd6fdxP7atkAky0XNbZXCP50Wy8ApiW71IwyiSMCDmeWYe7HB1BRq8OQaH88OZ6HTBOR81h6R2/0DvVBcZUWsz86gNN5lVKHhB3niyxLmy9O6YdhMfyAawtMtlzYY6PjcOeAMOgMAo+tPoL8yjrJYtmXVoL7PjmIqno9hnUPxMoHh/GQaSJyKp5uSny5eDgGRPihrKYBc/97AEeyyiWLJ6OkBk98eRRCAHOGRmF+UrRksTg7vpu5MJlMhjdmDUTvUB+UVGvx6KrDqNd1fsH8trOFuP+zVNQ2GDAqPgifPzgMPu6qTo+DiMjWArzcsGZxEoZ2D4CmXo/5y1OwL72k0+Oo1urx8KpD0NTrkRDtj5em9mMfQxtisuXiPN2U+HhBIvw9VTieU4nnvzvVqQenbjyRh0dWHUaD3ojb+nXF8oWJ8HBjYSYROS9fdxU+f3AYRsUHobbBgAdWpGL7ucJO+/lGo8Af1x/DhcJqhPio8eH8m6BW8nXXlphsEaICPfHevQmQy4BvjuTg830ZnfJz1x/KxhNfHoXeKDBtcDjeuzeBf/BE5BLMH3Qn9OkKrd6Ih1cexqYT+Z3ys9/7LQ0/ny6Em0KOD++7CSG+7p3yc10Zky0CANzcIwh/ucPUz+qVTWexP73Upj/vs72X8ezXJ2BsbFr6z9mDoVTw6UhErsNdpcAH8xMwZVA49EaBP3x5BF8dyrbpz9x6phD/3HoBAPC3af2REB1g059HJnx3I4tFyTG4e0gEDEaB331xBDnltTb5Oe/9lob/t8HU3+uh5Bi8end/diomIpekUsjx9j2DMWdoFIwC+NPXJ7Byf4ZNflZaUTWeXncMQgALRnTD7KFRNvk5dC0mW2Qhk8mwbPoA9I/wRVlNAx5ZdbhDj5YQQuD1Lefwxs/nAQBPjo/H83f2YVEmEbk0hdz02vvgzTEAgL/+cBrv70jr0J+hqdfh4VWHUKXVY1hMIF64q2+H3j9dn9XJ1rFjx7BgwQKMHj0aDz30ENLSbvyEaM2YttwvdTx3lQIf3ZeILl5uOJ2nwayP9mF9ajZqtPp23W9maQ2WfnsS7+9IBwA8f0cfPH1rTyZaREQwfdh94a4+eGJcDwDA61vO4//9eLrdKwyaeh3WpGRi1gf7cam4BuF+7nh/XgJULNvoVFb9tk+dOoXk5GR4e3vjueeeQ01NDYYPH46cnJx2jWnL/ZLtRPh74P15CXBXyXEqV4NnvzmBoX/fime/Po7DmWWt3q2YV1GHj3ddwtR392D0GzuwNjUbMpmpTmDxLbE2fhRERI5FJpPhmYm98NztvQGYjvZJfu03zPhgHz7bexlFVfWtuh8hBPanl+KZdccw7O9b8fx3p3C+sAreaiU+ui8RQd5qWz4MaoZMWLHP/5577kFubi727NkDADAYDOjVqxcmT56Mt99+u81j2nK//0uj0cDPzw+VlZXw9fVt7UOi6yjS1OObI7lYfygbl0tqLNfjgr0wOzEK0xMiEezT9I+2uEqLzafyseF4HlIzrjTrU8hlGBnXBQ/eHIOxvUM67TEQETmizSfzsXJ/Jg5cLoX5XVomA4bHdMHkQeG4vX8oArzcmowpqKzH14ez8dXhHGSWXpkRiw/xxj1Do3D3kAh0YaLVYazJO6xKtoKCgvCnP/0Jf/7zny3XnnzySWzfvh0nT55s85i23O//YrJlO0IIpGaUY/2hbGw6kY+6xsanCrkM43qHYOZNkaiobcCG4/nYl14C41UvDEO7B1peGPhpiojIOoWaevx00vQB9khWheW6Ui5DcnwQJg8Mh4ebAusPZWPXhWLL66+3WonJg8IwOzEKg6P8WbJhA9bkHcrW3mlNTQ1KS0sRHh7e5Hp4eDgyM5s/yLg1Y9pyvwCg1Wqh1Wot/6/RaFr7UMhKMpkMw2ICMSwmEC9O7ouNJ/Kx/lA2jmZV4Nczhfj1TNNmfIOi/DF5YBjuHBiGMD8eaEpE1FZdfd3xwM0xeODmGGSX1WJTY+J1Ok+DHeeLseN8cZPbD4sJxOzEKNwxIBSebq1+iycba/W/hE6nAwCo1U1nJzw8PCzfa8uYttwvACxbtgwvvfRSa8OnDuLjrsLcYdGYOywaFwqr8NWhbPx0sgC+HircNTAMkweGI7qLp9RhEhE5nahATzw6Og6Pjo5DenE1Nh7Px6aTeahtMGDyoHDMToxCTJCX1GFSM1qdbPn4+EClUqGsrKzJ9dLSUnTp0qXNY9pyvwCwdOlSPPPMM5b/12g0iIpiz5DO1LOrD56/sy+ev5NbiImIOlNcsDeenBCPJyfESx0KtUKrdyMqFAoMHDgQqampTa6npKRgyJAhbR7TlvsFTDNhvr6+Tb6IiIiI7I1VrR8WLVqEr7/+GmfOmLp/79mzB9u2bcOiRYsst1m+fDluvfVWq8a05jZEREREjsiq6rlHH30Up06dwpAhQ9C9e3dkZmZi6dKlmDZtmuU2OTk5TWapWjOmNbchIiIickRWtX4wKy4uRk5ODrp3746AgKaHWObk5KCgoACJiYmtHmPNbVrC1g9ERETUWWzWZ8ueMdkiIiKizmJN3sHDkYiIiIhsiMkWERERkQ0x2SIiIiKyISZbRERERDbEZIuIiIjIhphsEREREdmQ0xwJbu5godFoJI6EiIiInJ0532hNBy2nSbaqqqoAgIdRExERUaepqqqCn5/fdW/jNE1NjUYj8vLy4OPjA5lM1uH3r9FoEBUVhezsbDZN/R/83TSPv5eW8XfTPP5eWsbfTfP4e2mZrX83QghUVVUhPDwccvn1q7KcZmZLLpcjMjLS5j/H19eXT+gW8HfTPP5eWsbfTfP4e2kZfzfN4++lZbb83dxoRsuMBfJERERENsRki4iIiMiGmGy1klqtxosvvgi1Wi11KHaHv5vm8ffSMv5umsffS8v4u2kefy8ts6ffjdMUyBMRERHZI85sEREREdkQky0iIiIiG2KyRURERGRDTLZuoKamBi+++CLGjRuHKVOmYP369Tcc8/PPP+PBBx/E+PHjsWjRIhw6dKgTIrWdX375BTNnzsSYMWPwxz/+ESUlJTYZ42jKy8vx7LPPYuzYsZg+fTp++umn695eCIHvvvsOCxYswIQJE/Doo4/i9OnTnRRt5/ruu+8wbdo0jB07Fs8//7xVx2i99tprGD58OL7++msbRiiN/Px8/OEPf8CYMWNwzz33YNeuXa0a9+2332LWrFmYOHEi3nzzTej1ehtH2vlWrlyJO++8E+PHj8err76K+vr6696+vLwcr7zyCu68807ceuutePLJJ3Hp0qVOirbzpKen49lnn8WIESOwatWqVo9ZvHgxRo8ejfvuuw9HjhyxcZSdT6vVYs2aNbjtttswadKkVo05evQofv/732PChAm49957sXHjRhtHeQWTresQQuCuu+7CDz/8gCeeeAITJ07EggUL8OGHH7Y45tVXX8W//vUvjBo1Cn/5y18QFBSE4cOH4+eff+7EyDvOjz/+iDvvvBMJCQlYsmQJDh06hDFjxkCr1XboGEej0+kwbtw47N27F8888wySkpIwdepUfPXVVy2OWbJkCVavXo1bb70VS5cuhVwuR0JCAlJSUjoxcttbsWIF5syZgzFjxuCpp57Cli1bMGnSJBiNxhuO3b59O1asWIHjx4+joKCgE6LtPBqNBiNHjrS8efbo0QPjx4/Hjh07rjvumWeewcMPP4xx48bh+eefR0VFBV555ZXOCbqTLFu2DL///e8xbdo0PPbYY/jss88wZ86c64657bbb8O2332Lx4sV49tlnkZGRgeHDh6OwsLCTora9H3/8EZMmTUJgYCDS09ORm5t7wzF5eXkYMWIEqqqq8Nxzz8HX1xfJyck4ceJEJ0TceYYOHYqNGzciLCysVRMamzZtwuLFi9GnTx8sXboUQ4cOxezZs/Hmm292QrQABLXop59+EgDExYsXLddefPFFERQUJHQ6XbNjqqqqrrk2efJkcfvtt9ssTlvq16+feOSRRyz/X1paKtzc3MTy5cs7dIyjWblypVAqlaKoqMhy7Xe/+53o0aNHi2Oae26MGDFC3HfffTaJUQoGg0GEhYWJ559/3nLt8uXLQiaTie+///66Y4uKikRUVJRISUkRXl5e4p133rF1uJ3q9ddfF/7+/qKurs5ybcaMGSI5ObnFMbt37xYAxNatW5tcv/o+HF11dbXw9PRs8u+9f/9+AUAcOHCg2TG5ubkCgNiyZYvlWnl5uQAgvvnmG5vH3Fk0Go0wGo1CCCEiIiLEsmXLbjjmmWeeEbGxsUKv11uuJScni5kzZ9osTilUVFQIIYR4++23RZcuXW54++rq6muuvfDCCyIsLKzDY2sOZ7auY9u2bejduzd69OhhuTZ16lSUlJTg+PHjzY7x9vZu9lpDQ4PN4rSVgoICnD59GpMnT7ZcCwwMRHJyMrZu3dphYxzRtm3bMHz4cAQHB1uuTZ06FWlpacjMzGx2jDM9N1py5swZ5OfnN/n37969OwYMGHDdf38hBBYuXIiHH34Yw4YN64xQO922bdswYcIEuLu7W65NnToV+/btQ21tbbNjVq1ahZ49e2L8+PFNrl99H45u//79qK2tbfKcGT58OLp27dricyY4OBjR0dHYt2+f5dqePXvg5uaGQYMG2TzmztKWs363bduGO+64AwqFwnJtypQpTvX6C7T+mBwzLy+va6515usvk63ryMzMRHh4eJNr5v9v6Q31f509e9ZSv+JozI+xud9BS4+/LWMcUUc8Nw4cOIBt27Y55HOjJW3993/rrbeg0WiwdOlSm8YnpZaeM0ajETk5Oc2OOXPmDBITE7FixQrceuutmDp1Kv797387VYJufl6EhYU1uR4WFtbic0alUmHnzp3YsGEDYmNj0a9fPzz88MPYunUr4uLibB6zPWvpeVZRUWFV7aSzq6ysxHvvvddpr79OcxB1a+Tk5GDmzJnXvc2dd96JF154AYCpLud/O896eHhYvncjRUVFmDJlCsaOHYvHH3+8jVFLx/wYm/sdtPT42zLGEbX3uZGRkYEZM2Zgzpw5N6xNcSTX+/evq6trdkxqaipee+01pKamNvk07mza8pypr6/Hpk2bUF1djT/96U8oKSnBX/7yF+zevdtpNhDodDrI5XKoVKom16/3miGEwFNPPQWj0Yg333wTXl5e+PDDD/HYY49h165dCAwM7IzQ7VJ7X5tcQUNDA2bOnAl3d/dOq9lyqWQrODgY//rXv254G7PAwECcOXOmyfdLS0sBAF26dLnu/RQXF2P8+PGIjo7GN998A7nc8SYRzS9YZWVlTa6Xlpa2+PjbMsYRBQYGNvsYgRs/N7KysjBu3DgkJSXh888/t1mMUrj63z8kJMRyvbS0FFFRUc2O+eKLLwCgSdJZV1eHt956C9u2bcN3331nw4g7T1ueM4GBgZDJZFi/fr3lDdTLywvTpk1DXl7eNTMYjigwMBBGoxGVlZXw9/e3XL/ea8b27dvxww8/4OzZs+jduzcAYOzYsYiMjMSHH36Iv/zlL50Rul1q6XmmVCqtXnpzRuZEKyMjAzt27GjynLMll0q21Go1hg8f3urbJyQkYP369aivr7fUSKSkpEAul2PgwIEtjispKcH48eMRFBSEDRs2WD5VOJr4+Hj4+PggNTUVycnJAEyfKFNTUzFv3rwOG+OIEhIS8NFHH0EIYampSElJgYeHB3r16tXiuOzsbIwdOxYDBw7EunXroFQ615/gwIEDoVQqkZqaankTbGhowPHjxzFlypRmxzz11FO45557mlwbN24cZsyY4XTPmf379ze5lpKSgtDQUISGhjY7JjExEWlpaU1mKszLbeXl5U6RbCUkJAAwzXDeeuutAEzJwaVLlzBkyJBmx5iTiYiICMs1Nzc3BAUFXZNouJqEhASkpqY2uZaSkoL+/fs73euNtXQ6HWbNmoUzZ85gx44dTZ4/NtcpZfgOqrCwUPj4+IhXXnlFCCFETU2NuOmmm8TUqVMtt8nLyxNJSUli586dQgjTzruBAweKsWPHipqaGinC7lCPP/64iIuLE8XFxUIIIZYvXy4UCoU4e/as5TZLly4VDz/8sFVjHF1aWppQqVTigw8+EEIIUVZWJnr27CkWLVpkuc3Zs2dFUlKSOH78uBBCiJycHBEXFyemTZsmGhoaJIm7M8yePVsMHjzYsvvyH//4h/D09BS5ubmW2zz22GPi2WefbfE+nHE3YkpKipDJZOK7774TQgiRnZ19zc7NPXv2iKSkJJGdnS2EEOLixYtCrVaLTZs2CSGE0Ov14oEHHhCRkZEt7oh2RLfccosYO3as0Gq1QgghnnrqKREUFCQ0Go3lNvfcc4947bXXhBBCZGZmCrVabXltFkKILVu2CLlcfsNdr46qpd2IP/zwg0hKShK1tbVCCCF+/PFHoVQqxa5du4QQQpw5c0b4+vo63d+T2fV2I06YMEF8/PHHQgghdDqdmDZtmoiLi7P8fXUmJls38NNPP4mgoCDRrVs34evrK0aMGCEKCwst3798+bIAYHkB/cMf/iAAiIEDB4qkpCTL17Rp0yR6BO1TXV0t7rjjDuHp6Sni4uKEt7e3WLlyZZPbzJgxQ9x8881WjXEGa9euFb6+viI2NlZ4enqKW2+9VVRWVlq+n5qaKgCI3bt3CyGEmDt3rgAgEhISmjw37r//fqkegk2UlJSIUaNGCR8fHxETEyMCAgLEDz/80OQ2o0ePbvKh5X85Y7IlhBDvvfee8PT0FPHx8cLd3V3MmjVL1NfXW76/YcOGa9rNrF+/XgQGBoo+ffqI0NBQ0adPH5GamipF+DaTmZkpBg8eLAICAkRkZKQICwuzfIA169WrV5OWMmvXrhUhISEiOjpa9O7dW3h4eIgXXnihs0O3qYyMDMvrhJubm4iKihJJSUniqaeestzm448/FgCatJZ56aWXhFqtFj179hRubm7ikUceEQaDQYqHYDN//OMfRVJSkujWrZtQKpWW31N6errlNn5+fuLFF18UQgjxySefCAAiLi6uyetvUlJSs215OppMCCE6bx7NMTU0NODs2bPw9va+ZqeLVqvF0aNH0atXLwQEBODy5cvNNtVTq9UtTok7gqysLJSWlqJXr17w9PRs8r2LFy9Cp9Ohb9++rR7jLOrq6nD+/HkEBASgW7duTb5XU1ODkydPol+/fvDx8cHFixctNTpX8/b2Rv/+/Tsr5E5z+fJlaDQa9O7d+5qC3TNnzkCpVKJnz57Njk1NTUV0dDS6du3aGaF2qqqqKqSlpaFr167XLANWVFTg3LlzGDJkSJPfmVarxblz5+Dv74/o6Gir2wE4igsXLkCr1aJPnz7XLHkdP34cvr6+iImJsVzT6/XIyMhAQ0MDunfv7nSvM/X19Th27Ng11/39/S3L9MXFxUhPT8ewYcOa1AaXl5fj8uXLiIyMbFI/6SzOnz+P8vLya64PGjTIUrpz6NAhhIaGIjIyEkVFRS2eMDB06FCbb85hskVERERkQ463RY6IiIjIgTDZIiIiIrIhJltERERENsRki4iIiMiGmGwRERER2RCTLSIiIiIbYrJFREREZENMtoiIiIhsiMkWERERkQ259hHgROTShBA4cuQIcnJyMHDgwCZHwRARdRQmW0TkknJycjBt2jTk5uYiISEB58+fx4IFC/DXv/5V6tCIyMkw2SIil3TvvffCx8cHaWlp8PLygtFoxKZNm6QOi4icEA+iJiKXc/nyZcTGxmLHjh0YPXq01OEQkZNjgTwRuZysrCwAQM+ePSWOhIhcAZMtInI5/v7+AIDS0lJpAyEil8Bki4hcTr9+/RAeHo6VK1c2uV5cXCxRRETkzFggT0QuR6lU4qOPPsLMmTNRWFiI5ORknDlzBllZWfjmm2+kDo+InAwL5InIZZ09exYrV65EQUEBEhIS8PDDD0OtVksdFhE5GSZbRERERDbEmi0iIiIiG2KyRURERGRDTLaIiIiIbIjJFhEREZENMdkiIiIisiEmW0REREQ2xGSLiIiIyIaYbBERERHZEJMtIiIiIhtiskVERERkQ0y2iIiIiGyIyRYRERGRDf1/VD/W+YFfrVQAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "psplot.figure(figsize=(7,4))\n", "psplot.sympy_function(bulk_free_energy_density.subs(A, 1), (-0.2, 1.2))\n", "psplot.xlabel(\"c\")\n", "psplot.title(\"Bulk free energy\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To minimize the total free energy we use the Cahn Hilliard equation\n", "\n", "$$\\partial_t c = \\nabla \\cdot \\left( M \\nabla \\frac{\\delta F}{\\delta c} \\right)$$\n", "\n", "where the functional derivative $\\frac{\\delta F}{\\delta c}$ in this case is the chemical potential $\\mu$.\n", "A functional derivative is computed as \n", "$\\frac{\\delta F}{\\delta c} = \\frac{\\partial F}{\\partial c} - \\nabla \\cdot \\frac{\\partial F}{\\partial \\nabla c}$. \n", "That means we treat $\\nabla c$ like a normal variable when calculating derivatives.\n", "\n", "We don't have to worry about that in detail, since pystencils offers a function to do just that:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle {c}_{(0,0)}^{2} A \\left(2 {c}_{(0,0)} - 2\\right) + 2 {c}_{(0,0)} A \\left(1 - {c}_{(0,0)}\\right)^{2} - {\\partial_{0} (κ {\\partial_{0} {c}_{(0,0)}}) } - {\\partial_{1} (κ {\\partial_{1} {c}_{(0,0)}}) }$" ], "text/plain": [ "c_C**2*A*(2*c_C - 2) + 2*c_C*A*(1 - c_C)**2 - Diff(κ*Diff(c_C, 0, -1), 0, -1) - Diff(κ*Diff(c_C, 1, -1), 1, -1)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ps.fd.functional_derivative(free_energy_density, c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case we could quite simply do this derivative by hand but for more complex phase field models this step is quite tedious.\n", "\n", "If we discretize this term using finite differences, we have a computation rule how to compute the chemical potential $\\mu$ from the free energy." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle 4 {c}_{(0,0)}^{3} A - 6 {c}_{(0,0)}^{2} A + 2 {c}_{(0,0)} A - κ \\left(- 2 {c}_{(0,0)} + {c}_{(1,0)} + {c}_{(-1,0)}\\right) - κ \\left(- 2 {c}_{(0,0)} + {c}_{(0,1)} + {c}_{(0,-1)}\\right)$" ], "text/plain": [ "4*c_C**3*A - 6*c_C**2*A + 2*c_C*A - κ*(-2*c_C + c_E + c_W) - κ*(-2*c_C + c_N + c_S)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "discretize = ps.fd.Discretization2ndOrder(dx=1, dt=0.01)\n", "\n", "μ_update_eq = ps.fd.functional_derivative(free_energy_density, c)\n", "μ_update_eq = ps.fd.expand_diff_linear(μ_update_eq, constants=[κ]) # pull constant κ in front of the derivatives\n", "μ_update_eq_discretized = discretize(μ_update_eq)\n", "μ_update_eq_discretized" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "pystencils computed the finite difference approximation for us. This was only possible since all symbols occuring inside derivatives are pystencils field variables, so that neighboring values can be accessed.\n", "Next we bake this formula into a kernel that writes the chemical potential to a field. Therefor we first insert the $\\kappa$ and $A$ parameters, build an assignment out of it and compile the kernel" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "μ_kernel = ps.create_kernel([ps.Assignment(μ_field.center, \n", " μ_update_eq_discretized.subs(A, 1).subs(κ, 0.5))]\n", " ).compile()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we formulate the Cahn-Hilliard equation itself, which is just a diffusion equation:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle - div(M \\nabla \\mu) + \\partial_t c_{C}$" ], "text/plain": [ "-Diffusion(mu_C, M) + Transient(c_C)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M = sp.Symbol(\"M\")\n", "cahn_hilliard = ps.fd.transient(c) - ps.fd.diffusion(μ, M)\n", "cahn_hilliard" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It can be given right away to the `discretize` function, that by default uses a simple explicit Euler scheme for temporal, and second order finite differences for spatial discretization. \n", "It returns the update rule for the concentration field." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle {c}_{(0,0)} - 0.04 {μ}_{(0,0)} M + 0.01 {μ}_{(1,0)} M + 0.01 {μ}_{(0,1)} M + 0.01 {μ}_{(0,-1)} M + 0.01 {μ}_{(-1,0)} M$" ], "text/plain": [ "c_C - 0.04*mu_C*M + 0.01*mu_E*M + 0.01*mu_N*M + 0.01*mu_S*M + 0.01*mu_W*M" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c_update = discretize(cahn_hilliard)\n", "c_update" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, we build a kernel from this update rule:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "c_kernel = ps.create_kernel([ps.Assignment(c_field.center, \n", " c_update.subs(M, 1))]\n", " ).compile()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we run the simulation, the domain has to be initialized. To access a numpy array inside a data handling we have to iterate over the data handling. This somewhat complicated way is necessary to be able to switch to distributed memory parallel simulations without having to alter the code. Basically this loops says \"iterate over the portion of the domain that belongs to my process\", which in our serial case here is just the full domain. \n", "\n", "As suggested in the book, we initialize everything with $c=0.4$ and add some random noise on top of it." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def init(value=0.4, noise=0.02):\n", " for b in dh.iterate():\n", " b['c'].fill(value)\n", " np.add(b['c'], noise*np.random.rand(*b['c'].shape), out=b['c'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The time loop of the simulation is now rather straightforward. We call the kernels to update the chemical potential and the concentration in alternating fashion. In between we have to do synchronization steps for the fields that take care of the periodic boundary condition, and in the parallel case of the communciation between processes. " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def timeloop(steps=100):\n", " c_sync = dh.synchronization_function(['c'])\n", " μ_sync = dh.synchronization_function(['mu'])\n", " for t in range(steps):\n", " c_sync()\n", " dh.run_kernel(μ_kernel)\n", " μ_sync()\n", " dh.run_kernel(c_kernel)\n", " return dh.gather_array('c')\n", "init()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can run the simulation and see how the phases separate" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "if 'is_test_run' in globals():\n", " timeloop(10)\n", " result = None\n", "else:\n", " ani = psplot.scalar_field_animation(timeloop, rescale=True, frames=600)\n", "\n", " from pystencils.jupyter import display_as_html_video\n", " result = display_as_html_video(ani)\n", "result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now there are a lot of places to change and play with this model. Here a few ideas:\n", "\n", "- try different initial conditions and/or parameters $\\kappa, A$\n", "- the model can be generalized to 3D, by altering the DataHandling and the plot commands\n", "- modify the free energy formulation, make one minima lower than the other, maybe even add another phase" ] } ], "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.12.3" } }, "nbformat": 4, "nbformat_minor": 4 }