{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from pystencils.session import *" ] }, { "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": [ " ⎛ 2 2⎞\n", " 2 2 κ⋅⎝D(c[0,0]) + D(c[0,0]) ⎠\n", "c_C ⋅A⋅(1 - c_C) + ───────────────────────────\n", " 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": "iVBORw0KGgoAAAANSUhEUgAAAbUAAAEWCAYAAADhIgmdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl81NW5+PHPkz0hC0lIQsgKJEEiIEtYBQVRi1rFWkWsC1WstS23tz/b3tre1ttr23u73NbWW28tFituFZerUqVaxRVlC/sOIWSHkJB9387vj0y4aQxkEibz/c7M8369eDGZ75mZ55tk8sz3nOecI8YYlFJKKW/gZ3UASimllKtoUlNKKeU1NKkppZTyGprUlFJKeQ1NakoppbyGJjWllFJeQ5OaUucgIkZEMhy3nxKRnzr5uAQR+UhE6kXk18MbpVKqtwCrA1BquIhIAZAAdALtwKfA/caY4mF+6fuASiDS6ERQpdxKr9SUt7veGBMOJALlwH+74TXTgIPnSmgi4pEfJj01buVbNKkpn2CMaQFeBrJ77hORD0Tk3l5ff1lENg30XCISISLvi8ijIiJ9jj0FrAD+RUQaRORKEfmxiLwsIs+KSB3wZRHxE5EHReS4iJwRkRdFJKbX88wRkU9FpEZE9ojIwvPEM0ZEXhGRChE5ISLf7HXsx47nftrRHXpARHIG8di+cYeKyFoRqRaRQyLyLyJS4mj/XRF5pU9s/y0ivx3oe6qUq2hSUz5BRMKAW4EtF/g8scBG4BNjzDf7Xo0ZY74MPAf80hgTbox513FoKd1JdaTj+DeBG4HLgTFANfCY4zWSgDeBnwIxwHeAV0Qkrp94/IC/AnuAJGAx8C0R+VyvZjcALzheez3w+0E8tm/c/wakA+OAq4A7erV9FlgiIiMdzx9A9/f8mXN8O5VyOU1qytu9JiI1QB3df4R/dQHPNQb4EHjJGPPDQT52szHmNWNMlzGmGfgq8K/GmBJjTCvwY+BmRyK4A9hgjNngaP8OkAtc28/zzgTijDEPG2PajDH5wBPA8l5tNjmeq5PuBHPJIB7bN+5lwH8YY6qNMSXAoz0NjTEngY+AWxx3LQEqjTE7Bvm9UmrItI9cebsbjTHviog/3VcdH4pItjHm1BCe6zqgAXh8CI/tW5ySBrwqIl297uuku7AlDbhFRK7vdSwQeL+f500DxjgSdw9/4ONeX/c+1yYgxJE8nXls37jH9Lmv7/G1wNfoTo53oFdpys30Sk35BGNMpzHmf+lOHPMddzcCYb2ajR7gaZ4A3gI2iMiIwYbQ5+ti4BpjzMhe/0KMMaWOY8/0OTbCGPPzfp63GDjRp22EMaa/q7qhPLZv3CeB5F5fp/Q5/howRUQmAZ+nu8tSKbfRpKZ8gnRbCkQDhxx37wZuEpEwx3y0lU481SrgCPCGiIReQEiPAz8TkTRHfHGO+KB7bOp6EfmciPiLSIiILBSR5H6eZxtQJyLfcxRx+IvIJBGZ6UQMQ3nsi8D3RSTaMfa3qvfBXgU5zwPbjDFFTsShlMtoUlPe7q8i0kD3mNrPgBXGmAOOY48AbXSX+q/FiasKR2HIfXRf5bwuIiFDjOt3dBdt/F1E6ukuYJnteI1iurtKfwBUOF7ru/TzfnWMk10PTAVO0D0/7k9AlBPnMpTHPgyUONq/S3cCa+3TZi0wGe16VBYQnRuqlBoqEfkasNwYc3mv+1KBw8BoY0ydZcEpn6RXakopp4lIoohc6phnNwH4NvBqr+N+wAPAC5rQlBW0+lEpNRhBwB+BsUAN3fPf/gfAUTxTDhTSXc6vlNtp96NSSimvod2PSimlvIbtuh9HjRpl0tPTrQ5DKaWUjezYsaPSGPOZpeL6sl1SS09PJzc31+owlFJK2YiIFDrTTrsflVJKeQ1NakoppbyGJjWllFJeQ5OaUkopr6FJTSmllNfQpKaUUspraFJTSinlNWw3T00pX1TX0s7fD5RTdKbxnG2iwoJYMmk0SSMvZBs3pbybJjWlLNLa0ckHRyp4fXcp7x46TVtHFwAi/bc3Bn7yxkFmpcdw47Qkrp08mpFhQW6MWCn7s92Cxjk5OUZXFFHeqqvLsK2gitd3l7Jh3ylqm9uJHRHE9ZeMYenUMUxNGYmcI6sVnWni9d2lvLa7lOMVjQT6CwsnxHPj1CQWT4wnJNDfzWejlPuIyA5jTM6A7TSpKeUeuQVVfPflvZyobCQsyJ+rsxNYOi2J+RmjCPR3fnjbGMOBsjpe21XK+j1lnK5vJSIkgH+9diK3zkw5Z1JUypNpUlPKJto7u3h04zEeez+PpOhQvn3VBK6+OIGwoAvv/e/sMmzJP8Nj7+fx6fEzXJWdwM9vmkxseLALIlfKPjSpKWUD+RUNfGvdbvaW1HLLjGT+7YaLCQ92/VB2V5fhyU9O8Mu3jhAZGsivbpnCognxLn8dpazibFJzqs9DRJaIyBERyRORB/s5Hiwi6xzHt4pIeq9jU0Rks4gcEJF9IhIymBNRyhMZY3huayHXPbqJoqom/nD7dH51yyXDktAA/PyEexeM4/VVlxI7Ioi7/7ydh17fT3Nb57C8nlJ2NWBSExF/4DHgGiAbuE1Esvs0WwlUG2MygEeAXzgeGwA8C9xvjLkYWAi0uyx6pWyosqGVe9fm8q+v7icnPZq3/vkyrpmc6JbXnpgYyeurLmXl/LE8vbmQ63+/if2ltW55baXswJkrtVlAnjEm3xjTBrwALO3TZimw1nH7ZWCxdI9WXw3sNcbsATDGnDHG6EdH5bX2l9ay5Lcf8XFeJQ99Ppu1d89idJR7OydCAv350eezeXblbOpb2vnC/3zCS7nFbo1BKas4k9SSgN7viBLHff22McZ0ALVALJAFGBF5W0R2isi/9PcCInKfiOSKSG5FRcVgz0EpWzhyqp471mwlOMCfv66azz3zx+LnZ10l4vzMUbz9rcuYPTaWf3llL6/vLrUsFqXcxZmk1t+7sm91ybnaBADzgdsd/39BRBZ/pqExq40xOcaYnLi4AXfrVsp2jlc0cPufthAc4MfzX5nNhNERVocEwMiwIJ64K4fZY2N44MU9/G3fSatDUmpYOZPUSoCUXl8nA2XnauMYR4sCqhz3f2iMqTTGNAEbgOkXGrRSdlJ0ponbn9gKwHP3ziEtdoTFEf2j0CB/1qyYySXJUXzzhV28d7jc6pCUGjbOJLXtQKaIjBWRIGA5sL5Pm/XACsftm4H3TPdcgbeBKSIS5kh2lwMHXRO6UtYrq2nmtie20NLRyTMrZ5MRH251SP0aERzAU/fM4qLRkdz/7E42Hau0OiSlhsWASc0xRraK7gR1CHjRGHNARB4WkRsczdYAsSKSBzwAPOh4bDXwG7oT425gpzHmTdefhlLud7quhS89sYW65naeuWc2ExMjrQ7pvCJDAnn6nlmMGzWCe5/ezrYTVVaHpJTL6eRrpYbgTEMry1dvobSmmWdWzmJGWozVITmtor6V5as3c6q2hWfvnc201GirQ1JqQC6dfO1p9pfWaveKGja1Te3cuWYbRVVNrFkx06MSGkBcRDDP3TuH2PBgVjy5TeexqWFzoKyWj45W4M6LJ69Mao+9n8e/rd9vdRjKC3V1Gf553S7yTjew+q4c5o6PtTqkIRkdFcLzX5lNeHAA967NpaqxzeqQlBd6bmsR31q3262LbHtlUkuNCaO4upmuLnt1rSrP98TH+XxwpIIffX4il2d59vST5OgwVt+VQ1VjG999aY9bP00r31Bc1URKTJhbX9Mrk1pKTBhtHV2crm+1OhTlRXYWVfOrt49wzaTR3DEnzepwXGJSUhQ/uPYiNh4+zZpNJ6wOR3mZoqomUjWpXbieb2JRVZPFkShvUdvUzj89v4vRUSH8/ItTvGrPshXz0rk6O4FfvHWYPcU1VoejvERHZxel1c2kxoS69XU1qSk1AGMM33tlL+V1Lfz3bdOICg20OiSXEhF+efMU4iNCWPWXndS16Jrj6sKdrG2ho8volZorjBkZip9oUlOu8eyWQt46cIp/WTLBa8vfR4YF8ehtUymraeH7/7tPx9fUBSt2/P3VMTUXCArwIzEqlBJNauoCHSir5SdvHGLRhDjunT/O6nCG1Yy0GL5z9QTe3HuS57cVWR2O8nDF1d1/f/VKzUVSY8L0Sk1dkIbWDv7p+V1Ejwjk18umWrrivrt89bJxXJYVx8N/Pcihk3VWh6M8WFFVEwF+QmKUjqm5hCY1dSGMMfzotf0UnGnkd8unETMiyOqQ3MLPT/jNskuIDA1k1fM7aWrrsDok5aGKqppJjg7F380fBr02qaXEhHK6vlW3s1dD8srOUl7dVco/L85izjjPnGA9VKPCg/ndrVPJr2zkodcPWB2O8lBFFsxRA69Oat3fzJJqvVpTg3O6voV//+sBZo+NYdUVGVaHY4l5GaNYtSiDl3eU8NFR3bhXDZ4VE6/Bi5OalvWrofrPDYdpbe/iP2+a7PauEztZdUUG40aN4KHX99PSrj0eynn1Le1UNba5vUgENKkp9Q82Hz/Dq7tKuf/ycYyLs+feaO4SHODPw0snUXCmiT9+mG91OMqDFFc1A+6vfAQvTmoxI4IYEeSvSU05ra2jix+9vp+UmFC+vsg3ux37mp85iusvGcNjH+RRUNlodTjKQ/T83dWk5kIiQkpM2NkJgEoN5E+b8sk73cDDN0wiJNDf6nBs44fXTSTI34+H1h/QSdnKKVZNvAYvTmqgZf3KeSXVTTy68RhLLh7NoovirQ7HVhIiQ/j21Vl8dLSCt/afsjoc5QGKqpqICg20ZEk5r09qxVXN+ulSDejH6w/iJ8JD12dbHYot3TknjezESP79rwdpaNW5a+r8iqvdvzp/D+9OarFhNLd3UtmgGyCqc3vnYDnvHirnW1dmMmake1c/8BQB/n787AuTKK9v4XfvHrU6HGVzVmw508Ork1qKVkCqATS1dfDj9QfISgjn7kvHWh2OrU1LjWb5zFSe/KRAl9BS59TVZSiparZkPA28PKn1fFLQYhF1Lr9/L4/SmmZ+euNkAv29+u3gEt9bMoGo0EB++Np+3Vle9au8voW2zi69UhsOSY6uJL1SU/3JO13PEx/nc/OMZGaNjbE6HI8wMiyI719zETsKq3l5Z4nV4SgbKjrTU/loTVe+U0lNRJaIyBERyRORB/s5Hiwi6xzHt4pIuuP+dBFpFpHdjn+Puzb88wsJ9Gd0ZIgmNfUZ3QsWHyAsKIDvX3OR1eF4lC9OT2ZmejT/ueEQNU06Xq3+kZVz1MCJpCYi/sBjwDVANnCbiPQtEVsJVBtjMoBHgF/0OnbcGDPV8e9+F8XtNC3rV/15/8hpNuef4TtXZxEbHmx1OB7Fz0/4yY2TqGlu5w8fHLc6HGUzxVVN+AmWFV05c6U2C8gzxuQbY9qAF4ClfdosBdY6br8MLBYRWyyapxOwVV9dXYZfvnWE9Ngwls9KtTocj3TR6EhumpbMU58WcKq2xepwlI0UVTUxZmSoZWPUzrxqElDc6+sSx339tjHGdAC1QM9+HWNFZJeIfCgiCy4w3kFLjQnjVF2LLsiqzvrr3jIOn6rngasnaHHIBfjWlZl0GcPvNh6zOhRlI1aW84NzSa2/K66+ZU/nanMSSDXGTAMeAJ4XkcjPvIDIfSKSKyK5FRWu3eYiNTYUY6C0ptmlz6s8U1tHF7/++1GyEyP5/OREq8PxaCkxYdw+O40Xc4vJr2iwOhxlE0VVzbZPaiVASq+vk4Gyc7URkQAgCqgyxrQaY84AGGN2AMeBrL4vYIxZbYzJMcbkxMXFDf4szkNX61e9rdteRFFVE99dMgE/H95WxlVWXZFBcIAfv35HJ2Sr7nmflQ2tls1RA+eS2nYgU0TGikgQsBxY36fNemCF4/bNwHvGGCMicY5CE0RkHJAJuHUPi7ObhWpS83lNbR08+l4es8bGsDDLtR+efNWo8GDunT+WN/eeZH9prdXhKIuVVFu35UyPAZOaY4xsFfA2cAh40RhzQEQeFpEbHM3WALEikkd3N2NP2f9lwF4R2UN3Acn9xpgqV5/E+cSFBxMS6KdXaoo/f1JARX0r31syAZvUMXmFey8bx8iwQH759hGrQ1EW65mjZmVSC3CmkTFmA7Chz30P9brdAtzSz+NeAV65wBgviIhoWb+itqmdP354nCsnxjMjTSdau1JkSCDfWJjBzzYcYvPxM8wdHzvwg5RXsnqOGnj5iiI9UqLDKKrSQhFf9ocPj1Pf2sF3PjfB6lC80p1z0xgdGcIv3z6su2L4sKKqJsKDAxgZ5v4tZ3r4RlJzzFXTN5tvKq9r4c+fnODGqUlcNPozxbfKBUIC/fnWlZnsKqrhnYPlVoejLFJc1URKTJil3fs+kdRSY8JoaO2guqnd6lCUBR7deIzOLsP/u/IzhbfKhW6ekcy4USP4r78foVMXO/ZJ3XPUrN2+yWeSGmhZvy8qqGxk3fZivjQ7ldRY6/r5fUGAvx/fvnoCR8sbeG1XqdXhKDczxlg+8Rp8JanFalLzVb955yiB/n6suiLD6lB8wjWTRjM5KYpH3j1Ka4eu4uNLKupbae2wbsuZHj6R1FKidV81X3ToZB3r95Rxz/x04iNCrA7HJ/j5Cd/93ARKqptZt7144Acor9Fz0WDlxGvwkaQWGuRPXETw2TkUyjf8zwfHCQ8O4L4F460OxacsyBzFzPRoHv/gOG0dXVaHo9zEDuX84CNJDbq/0cXVmtR8xYnKRt7cW8Ydc9KIsrC82BeJCF9flEFZbQuv79axNV9RXNWMCCRFa6GIW+gEbN/y+AfHCfT3Y+X8sVaH4pMWZsWRnRjJHz44rpWQPqKoqonEyBCCA/wtjcNnklpKTBhlNc20d2p3iLcrq2nmf3eVcOvMFOIidANQK4gI31iUQX5lI2/tP2V1OMoNeuaoWc1nklpqTBhdpvsPnvJuT3ycjzFw32XjrA7Fpy2ZNJpxo0bw2Pt5uvCBD7BDOT/4UFJLcfTzahekdzvT0MpfthWxdGoSydHWv8F8mb+fcP/C8Rw8WccHR127T6Kyl5b2Tk7VteiVmjvpXDXf8OdPCmjt6OJrC/UqzQ5unJrEmKgQ/uf9PKtDUcPIDlvO9PCZpJYQEUKQv25B483qWtpZu7mAJRePJiM+wupwFBAU4Md9l41je0E12064ddcp5UbFNpmjBj6U1Pz8hOSYUJ2A7cWe3VJIfUsH31ikq4fYyfJZqYwKD+IxvVrzWnaZowY+lNRAy/q9WXNbJ2s+PsHlWXFMSoqyOhzVS0igP/fMH8uHRyt0d2wvVVTVRGigP6PCg6wOxQeTmq4q4pXWbS/iTGObXqXZ1B1z0ogICdCrNS/VU/lohx3lfS6p1bV0UKtb0HiVto4uVn+Uz8z0aGaN1V2t7SgyJJAVc9N568Ap8k7XWx2OcjG7zFEDH0tqPd90XS7Lu7y2u5Sy2ha+rldptnb3pekEB/jxhw/yrQ5FuZAxhmKbzFEDH0tquq+a9+nsMjz+wXGyEyNZmBVndTjqPGLDg7ltViqv7S7Vgi0vUtXYRmNbp+Wbg/bwqaSWoknN67y1/xT5lY18Y1GGLfrz1fl9ZcE4/KR71RflHc5WPtpkE16fSmrhwQHEjgjSpOYljDGs/jif9NgwlkwabXU4ygljRoZy49QkXswtpqapzepwlAvYqZwffCypASTHhGnXh5fYWVTNnuIa7pk/Fn8/vUrzFCsXjKWlvYvnthZZHYpygZ6/p3ZZls6ppCYiS0TkiIjkiciD/RwPFpF1juNbRSS9z/FUEWkQke+4Juyh07lq3mPNphNEhgTwxenJVoeiBuGi0ZHMzxjF05sLdBNRL1BU1UR8RDAhgdZuOdNjwKQmIv7AY8A1QDZwm4hk92m2Eqg2xmQAjwC/6HP8EeBvFx7uhUuNCaW0upkO3YLGoxVXNfHW/lN8aXYaI4IDrA5HDdLKBWMpr2vlzX1lVoeiLpBdVufv4cyV2iwgzxiTb4xpA14AlvZpsxRY67j9MrBYHKP2InIjkA8ccE3IFyYtdgQdXebsApzKMz31aQF+IqyYl2Z1KGoILs+MY3zcCNZsOqHb0ni4E5WNpMWOsDqMs5xJaklAca+vSxz39dvGGNMB1AKxIjIC+B7w7+d7ARG5T0RyRSS3omJ4t6jISuhe6PZouU4A9VT1Le2s217MdVMSSYyyRxmxGhw/P2Hl/HHsL63ThY49WG1zO+V1rWQlhFsdylnOJLX+RuD7frQ6V5t/Bx4xxjSc7wWMMauNMTnGmJy4uOGda5QR3/3NP3b6vCEpG1u3vZiG1g5Wzh9rdSjqAtw0PYnosED+tOmE1aGoITrmuDjI9LCkVgKk9Po6GejbEX62jYgEAFFAFTAb+KWIFADfAn4gIqsuMOYLEh4cQNLIUL1S81AdnV089WkBs9JjmJI80upw1AUICfTnjjlpvHuonILKRqvDUUNwtLz74iDTRls9OZPUtgOZIjJWRIKA5cD6Pm3WAysct28G3jPdFhhj0o0x6cBvgf8wxvzeRbEPWWZC+NkfhvIsfz9YTkl1M/foVZpXuHNOGgF+wp8/0as1T3S0vJ6wIH+SRtpnGGDApOYYI1sFvA0cAl40xhwQkYdF5AZHszV0j6HlAQ8Anyn7t5OshAiOVzTQ2aUD1J5mzaYTpMaEcVV2gtWhKBeIjwzhhkuSeGlHiS407oHyTjeQGR+On43miTo1T80Ys8EYk2WMGW+M+ZnjvoeMMesdt1uMMbcYYzKMMbOMMZ9ZA8cY82NjzH+5NvyhyYwPp62jS+ereZhdRdXsKKzm7kvTdbK1F1k5fyxNbZ38ZbtOxvY0R8vryUywT9cj+OCKIqAVkJ5qzaYTRIQEcEtOysCNlcfIHhPJvPGxrP20gHadP+oxapvaOV1vr8pH8NGkdrYCUpOaxyitaeZv+09x26xUwnWytddZOX8sJ2tb2LDvpNWhKCcddeyLZ6ciEfDRpDbibAWkFot4irWfFgCwYl66pXGo4bFoQjzjRulkbE9y1Ibl/OCjSQ0gKyFcux89RENrB3/ZWsQ1k0bbqspKuY6fn3D3/LHsLaklt7Da6nCUE46VNzDCZpWP4NNJLYL8ykZdA9IDvJRbTH1rB/cuGGd1KGoYfXF6EiPDAlnzsZb3e4Jjp+vJSIiw3T6GPpvUMhMitALSA3R1GZ76tIAZadFMTdHJ1t4sLCiAL81K5e8HT+n2UB7gaHkDWfH26noEX05qjh+GjqvZ24dHKyg808SXdSzNJ9wxJw0R4dmthVaHos6jpqmNivpW242ngQ8nNa2A9AxrNxcQHxGsO1v7iDEjQ7k6O4F124tpae+0Ohx1DmeXx7LZHDXw4aQ2IjiA5OhQjurCxrZ1orKRD45UcPvsNAL9ffZX1eesmJdOTVM763frXmt21VNkl6VJzV6yEiL0Ss3GntlcSKC/cNtsnWztS2aPjWFCQgRPfVqg5f02day8nvDgAMZEhVgdymf4dFLLTAgnv0IrIO2osbWDl3KLuXZyIvER9nvjqOEjIqyYl87Bk3Xs0PJ+Wzp2uoGM+HDbVT6Cjye1rPgI2jq7KNRKK9t5dVcp9a0d3DU33epQlAVunDaGyJAA1m7WghE7OlreYLvlsXr4dFLrqdzRLkh7Mcbw9OYCJidFMT1Vy/h9UVhQAMtyUvjbvpOU17VYHY7qpbqxjcqGVtstj9XDp5Nahpb129Lm/DMcLW/grrlptuzeUO5xx5w0Oo3h+a26er+d2HV5rB4+ndTCggJIidFdsO1m7acFRIcFcv0lY6wORVkofdQIFmbF8fy2Ito6dNzbLnoqxu1Y+Qg+ntSge1ztmF6p2UZpTTPvHCxn+axUQgL9rQ5HWWzFvHQq6lv5235dvd8ujpXXExEcQKINKx9BkxqZCRHkVzZoBaRNPLuluzDg9tmpFkei7OCyzDjGjhpxdpcGZb1j5Q1kJNiz8hE0qZEZH057p6HgjFZAWq2lvZMXthVxVXYCydFhVoejbMDPT7hzTho7i2rYV1JrdTiK7oWMM2245mMPn09qPf3CWgFpvb/uKaO6qZ0VWsaverk5J5mwIH/Wbi6wOhSfV9XYRmVDm23H00CTmmMCoVZAWs0Yw9rNBWTGhzN3fKzV4SgbiQwJ5KbpSazfU0ZVY5vV4fi0/6t81KRmW6FB/qREh53dmlxZY2dRDftL67hrXrpt++qVde6am05bRxcvbNfyfisdO7vmo3Y/2lpWQrh2P1rs6c0FRAQHcNO0JKtDUTaUlRDBvPGxPLelSIu6LHS0vIGI4ABGR9qz8hE0qQHdl9InKhtp1zeLJSrqW9mw7yQ35yQzIjjA6nCUTd01N53SmmY2Hj5tdSg+69jpejJtXPkITiY1EVkiIkdEJE9EHuzneLCIrHMc3yoi6Y77Z4nIbse/PSLyBdeG7xo9FZCFZxqtDsUnvbCtiPZOw51z0qwORdnYlRPjGRMVwjO6HqRljpU32HZ5rB4DJjUR8QceA64BsoHbRCS7T7OVQLUxJgN4BPiF4/79QI4xZiqwBPijiNjuo3hPJY8Wi7hfR2cXz28rYkHmKMbF2befXlkvwN+PL81OZVNeJccr9L3qbmcaWjnT2Gbb5bF6OHOlNgvIM8bkG2PagBeApX3aLAXWOm6/DCwWETHGNBljOhz3hwC23BxpfFxPBaSOq7nbu4dOc7K2Ra/SlFNunZlKoL+cnaSv3KfnQ7+dy/nBuaSWBBT3+rrEcV+/bRxJrBaIBRCR2SJyANgH3N8ryZ0lIveJSK6I5FZUVAz+LC5QaJA/qTFhulyWBZ7ZUsCYqBCuuCje6lCUB4iLCOaaSYm8vKOEprbP/ClRw+jYafvudt2bM0mtvxHBvldc52xjjNlqjLkYmAl8X0Q+UzZjjFltjMkxxuTExcU5EZLrZcZH6JWam+WdbuCTvDPcPieNAH+tWVLOuWtuGvUtHby+u8zqUHzK0fJ6IkICSIgMtjqU83LmL0kJkNLr62Sg72/T2TaOMbMooKp3A2PMIaARmDTUYIdTVkK4VkC62bNbCgn0F26dmTJwY6UcZqRFMzExkqc3F2KMLUc0vNKx8gayEiJsXfkIziW17UCmiIwVkSBgObC7ZyOqAAAgAElEQVS+T5v1wArH7ZuB94wxxvGYAAARSQMmAAUuidzFMhPC6egyFFRqBaQ7NLV18MqOEq6dnMiocHt/8lP2ItK9HuShk3XsLKq2Ohyfcex0g63XfOwxYFJzjIGtAt4GDgEvGmMOiMjDInKDo9kaIFZE8oAHgJ6y//nAHhHZDbwKfN0YU+nqk3CFnjJVrYB0j9d2lVHf2sFdc7VARA3ejdPGEBEcwNNa3u8WlQ2tVDW22Xp5rB5OldcbYzYAG/rc91Cv2y3ALf087hngmQuM0S0y4sPxc1RAXkei1eF4NWMMT28uYGJiJNNTo60OR3mgsKAAvjgjmee2FvLD67KJi9Cr/eF01AOWx+qho/MOIYGOCkhdA3LY7Sis5vCpeu6am2b7/nllX3fOTaO90/BibvHAjdUFOeYh5fygSe0fZCZEaPejGzy9uZCIkACWTh1jdSjKg42PC2d+xiie21Ko60EOs6Pl9USGBBDvAVfEmtR6yYwPp6CykbYOfYMMl4r6Vv62/yQ3z0gmLMh2i8soD3PHnDTKalt0Pchhdux0A5keUPkImtT+QVZCRHcFpK4BOWzWbdd1HpXr9KwHqSuMDB9jDMfK6z1iPA00qf2Dnv7iQyfrLI7EO3V0dvHcVl3nUblOz3qQHx+rJF/XgxwWp+tbqW5q94jxNNCk9g+yEsIJCfRjd3GN1aF4pZ51Hu/QqzTlQj3rQT6jV2vDYldR99/DS1JGWhyJczSp9RLg78eU5JFnf4jKtZ7dUsiYqBAW6zqPyoV0Pcjhtau4miB/Py4eE2l1KE7RpNbHtNSRHCyro7Wj0+pQvMrxigY25VXqOo9qWPSsB/naLl0P0tV2FdWQPSaS4AB/q0Nxiv516WNaSjRtnV0cKNNxNVd6ZnMhQf5+LMvRdR6V6/3fepAFuh6kC3V0drG3pIZpqZ7R9Qia1D5juuOHt7NQ15RzlYbWDl7eUcJ1UxJ15Qc1LESEFXPTOHyqnu0F+t51lcOn6mlp72KaB638o0mtj/jIEJJGhrJLi0Vc5tWdJTToOo9qmC2dmkRkSABrNxdYHYrX2OVYMHq6Xql5tmmpI9mtxSIuYYxh7eZCpiRHMdVDqqeUZwoN8ufWmSm8vf8Up2pbrA7HK+wqqiEuIpikkaFWh+I0TWr9mJYaTWlNM+V1+sa4UJuPnyHvdAN3zU33iNUIlGe7Y04ancbw/LYiq0PxCruKa5iWMtKj3rua1PrRMyiqpf0Xbu3mAqLDAvn8FN35QA2/tNgRLJoQz/Nbi3S5uwtU3djGicpGjxpPA01q/bp4TCRB/n7sKtYB5wtRWtPMOwfLWT4rlZBAzygHVp7vrrlpVDZ0rzGqhq5nEQpPqnwETWr9Cg7wJ3tMpF6pXaDnHCs83D471eJIlC+5LDOO9Ngw3UD0Au0qqsZPYEpylNWhDIomtXOYljqSvSU1uqXFELW0d/LC9mKunJhAcnSY1eEoH+LnJ9w5N50dhdXsL621OhyPtau4hotGR3rcbhqa1M5hWmo0Le1dHD6lm4YOxZt7T1LV2MaKeelWh6J80M0zkgkN9OfpzQVWh+KRuroMu4s8a9J1D01q5zAtpadYRMfVhuLpzQWMjxvBvPGxVoeifFBUaCBfmJ7E67vLqG5sszocj5NX0UB9a4fHFYmAJrVzSo4OJS4iWMfVhmB3cQ17SmpZMU/L+JV17pqbRmtHFy/mFlsdisfxxEnXPTSpnYOIMC1lpK4sMgRPf1pAeHAAN01PtjoU5cMuGh3J7LExPLOlkM4uXQ9yMHYV1RAVGsjYUSOsDmXQNKmdx7TUaE5UNmr3xSBUNrTyxt6TfHF6EuHBnjXArLzPinnplFQ38/7h01aH4lF2OcbTPLGnxamkJiJLROSIiOSJyIP9HA8WkXWO41tFJN1x/1UiskNE9jn+v8K14Q+vnkFS3TTUeeu2F9PW2cWdc9OtDkUprspOYHRkiK4HOQj1Le0cPV3PtBTPG08DJ5KaiPgDjwHXANnAbSKS3afZSqDaGJMBPAL8wnF/JXC9MWYysAJ4xlWBu8OU5Cj8RItFnNXR2cWzWwq5NCOWjPhwq8NRikB/P740O5WPj1VyvKLB6nA8wt6SWozxvEnXPZy5UpsF5Blj8o0xbcALwNI+bZYCax23XwYWi4gYY3YZY3p27TsAhIiIx+w9EhYUwEWjI3VczUnvHirnZG0Ld+lVmrKR5bNSCPQXntHJ2E7p+RB/iYcuQO5MUksCepcPlTju67eNMaYDqAX61nJ/EdhljGkdWqjW6Fmxv0sHmgf05KYCkkaGsviieKtDUeqs+IgQrpucyEu5xdS1tFsdju3tKqohIz6cqNBAq0MZEmeSWn8jhX3/wp+3jYhcTHeX5Ff7fQGR+0QkV0RyKyoqnAjJfaalRlPf2kGedl2c176SWrYVVHH3pekE+Gv9kbKXlfPH0djWyYvbtbz/fIwxZ1fm91TO/PUpAVJ6fZ0MlJ2rjYgEAFFAlePrZOBV4C5jzPH+XsAYs9oYk2OMyYmLixvcGQyz/1uxX8fVzmfNpnxGBPmzbGbKwI2VcrPJyVHMSo/hz58U6NJ351F4pomqxjaPnHTdw5mkth3IFJGxIhIELAfW92mznu5CEICbgfeMMUZERgJvAt83xnziqqDdadyoEUSFBuok7PM4VdvCG3tPsmxmCpEhntllobzfPfPHUlrTzN8Pllsdim317EwyPc2Lr9QcY2SrgLeBQ8CLxpgDIvKwiNzgaLYGiBWRPOABoKfsfxWQAfxIRHY7/nnUgIuIMC11pCa183h6cwGdxnD3vLFWh6LUOV2VnUBqTBhrNp2wOhTb2lVUw4ggfzLjI6wOZcicmh1rjNkAbOhz30O9brcAt/TzuJ8CP73AGC03LSWaD48epb6lnQi9EvkHzW2dPL+tiKuzE0iN1dX4lX35+wlfnpfOw28cZHdxDVM9eNxouOwqquGSlJH4+3nepOseOqLvhGmpIzGme/6G+kev7CyhpqmdlfPHWR2KUgNaNjOFiOAAvVrrR3NbJ4dO1nns/LQemtSccImu2N+vri7Dk5+cYHJSFDPTPXdgWfmO8OAAbp2ZwoZ9JymrabY6HFvZX1ZLR5fx2JVEemhSc0JUaCAZ8eE6rtbHh0cryK9oZOX8sR65RpzyTSvmpWOM0aWz+uj50D5Vr9R8Q8+K/cboJOweazadICEymGsnJ1odilJOS4kJY8mk0fxlaxGNrR1Wh2MbOwtrSI0JY1S4xyz61C9Nak6alhpNVWMbhWearA7FFg6fqmNTXiV3zU0nKEB/jZRnWTl/LHUtHbyys8TqUGzBGMPOomqPH08DTWpOOzsJu1jH1QCe3HSCkEA/bp+danUoSg3a9NRoLkkZyZ8/KdAl8ICTtS2crm/16JVEemhSc1JWQgThwQFsO6FJrbKhldd2l/HF6cmMDAuyOhylBk1EWDl/LCcqG3lP91pje0EVANPTPLtIBDSpOc3fT5ifMYr3D5/2+XG1Z7cU0tbRxT3zdbK18lzXTBpNYlSIlvcDGw+dJnZEEBePibI6lAumSW0QFk+M51RdCwfK6qwOxTIt7Z08u6WQRRPiGB+ne6YpzxXo78eKeelszj/DgTLfnYPa3tnFB0dOs+iieI+edN1Dk9ogXHFRPCLwjg+vHbd+TxmVDW062Vp5hdtmphIa6M+TmwqsDsUy2wuqqGvp4MqJCVaH4hKa1AYhNjyY6anRbDzsm0mtq8vwp4/zuWh0BJdm9N0uTynPExUWyLKcZNbvKeVkrW9Oxt546DRB/n4syBxldSguoUltkK6cmMD+0jqffAP8/WA5R8sb+NrC8TrZWnmNexeMo8vA6o/yrQ7F7YwxvHuonHkZsYwIdmopYNvTpDZIV07s3mRg4yHfqpgyxvD794+RFhvGdTrZWnmRlJgwbpyaxF+2FVHZ0Gp1OG51vKKBwjNNLPaSrkfQpDZoGfHhpMWG8e4h3+qC/PBoBftL6/j6wvG6s7XyOl9fNJ7Wji6fq4R852D3h/OeD+veQP86DZKIsPiiBD49fsZnltgxxvD79/IYExXCF6YlWx2OUi43Pi6caycn8szmQmqb2q0Ox202Hirn4jGRJEaFWh2Ky2hSG4Irs+Np6+ji42OVVofiFltPVJFbWM1XLx+vS2Ipr7VqUQYNrR089WmB1aG4xZmGVnYUVXtN1WMP/Qs1BDPTY4gICWCjj3RBPvZ+HqPCg7l1ZorVoSg1bCYmRnLlxHj+/OkJn+iFef9IBcagSU11T9pcNCGe9w6fptPL143bXVzDx8cquXfBWEIC/a0OR6lh9Y1FGdQ0tfPc1kKrQxl27x4sJyEymElJkVaH4lKa1IZo8cR4zjS2sbvYu/dY+/17eUSFBnLHnDSrQ1Fq2E1LjWZ+xihWf3SClvZOq8MZNi3tnXx0rILFExO8bnqOJrUhWpgVT4CfeHUV5OFTdbx7qJy7L00n3EvmsCg1kG8syqCyoZUXc4utDmXYbMk/Q1NbJ1d5WdcjaFIbsqiwQGamx3j1uNpj7x9nRJA/X56XbnUoSrnNnHExzEiL5o8f5tPW0WV1OMNi46HThAb6M3e8960MpEntAlyZncDR8gaKvHDj0PyKBt7YW8adc9N1exnlU0SEVVdkUFrTzGu7Sq0Ox+WMMWw8VM6CzFFeOU6uSe0C9ExY9MYuyD98cJwgfz9W6vYyygctzIpjUlIkf/jwuNcVgx08WUdZbYvXVT32cCqpicgSETkiInki8mA/x4NFZJ3j+FYRSXfcHysi74tIg4j83rWhWy8tdgQZ8eFet8BxSXUTr+4q5bZZqcRFBFsdjlJuJyJ8Y2EGJyobeXPfSavDcamNh04jAosu8p5VRHobMKmJiD/wGHANkA3cJiLZfZqtBKqNMRnAI8AvHPe3AD8CvuOyiG3myokJbM2voq7Fe1Yh+OOH+YjAfZfp9jLKd33u4tFkxIfz2Ht5dHnR1dq7h8qZmjLSaz+wOnOlNgvIM8bkG2PagBeApX3aLAXWOm6/DCwWETHGNBpjNtGd3LzSlRPj6egyfHikwupQXKK4qol1ucV8cXoyY0Z6z9I5Sg2Wn5+walEGR8rrveZqrbyuhb0ltV7b9QjOJbUkoHdta4njvn7bGGM6gFrA6bIaEblPRHJFJLeiwrOSw7TUaGJGBHnNuNqv3j6Cn8A/X5lpdShKWe76S8Zw0egIfvHWYVo7PH/eWs/uIr6e1Pqbmdf3WtyZNudkjFltjMkxxuTExcU5+zBb8PcTFk2I5/3Dp2nv9Ozy393FNazfU8ZXFozzqgVOlRoqfz/hh9dlU1LdzFovWBNy46FykqNDyUoItzqUYeNMUisBei/6lwyUnauNiAQAUUCVKwL0BFdlx1PX0kFuQbXVoQyZMYafvnGQUeHBfPXy8VaHo5RtzM8cxcIJcfz3e3lUNbZZHc6QNbd1simvkiu9cBWR3pxJatuBTBEZKyJBwHJgfZ8264EVjts3A+8ZY7xnZHUACzLjCPL38+iJ2G8fOEVuYTUPXJWlq4co1ccPrp1IY2sHj248ZnUoQ7Ypr5LWji6v7noEJ5KaY4xsFfA2cAh40RhzQEQeFpEbHM3WALEikgc8AJwt+xeRAuA3wJdFpKSfykmPNyI4gHkZsby576RHdkG2dXTx878dJishnGU5ul+aUn1lJUSwfFYqz24pJL+iwepwhuT13aVEhAQwa2yM1aEMK6fmqRljNhhjsowx440xP3Pc95AxZr3jdosx5hZjTIYxZpYxJr/XY9ONMTHGmHBjTLIx5uDwnIq17pyTxsnaFt7c63lVUs9sKaTgTBPfv3ai7mqt1Dn8vyuzCA7w4+d/O2x1KINWXNXE3/af4rZZqV6/J6J3n50bLZoQT0Z8OH/8KB9P6nmtbWrn0Y3HWJA5ioVZnlWko5Q7xUUE87WF4/n7wXK25p+xOpxBWbPpBALcfWm61aEMO01qLuLnJ9y3YByHTtaxKc9zdsT+7/eOUdfSzg+unejVg8dKucLK+eNIjArhZxsOecyE7JqmNtZtL+aGqWN8oqpZk5oLLZ02hriIYFZ/lD9wYxsoPNPI2s0F3DIjmYmJ3rVRoFLDITTIn+9cPYG9JbWs39O3CNyent1SSHN7J19Z4BsrBGlSc6HggO5tWj4+VsnBsjqrwxnQL986QoCfH9++eoLVoSjlMb4wLYlJSZH86u0jtt9ItKW9k6c+LeSyrDif+eCqSc3F7pidRliQP098bO+rtR2FVby57yRfvXwcCZEhVoejlMfw8xN+cO1ESmuaefKTE1aHc16v7SqlsqGVr/rQOq6a1FwsKiyQ5TNT+eueMspqmq0Op1/GGH765iHiI4J10WKlhmDe+FFcOTGe/3n/OJUNrVaH06+uLsPqj/O5eEwk87xwM9Bz0aQ2DO6Zn44Bntxkz09xz28rYldRDd+5egJhQTrRWqmhePCaibS0d/LwX+05S2nj4dPkVzRy32XjfKoITJPaMEiODuPzUxL5y7YiapvttSXN8YoGfvLGQRZkjuLmGTrRWqmhyogP558XZ7J+T5ktd8he/dFxkkaGct3kRKtDcStNasPkvsvG0djWyV+2FVkdylltHV1864XdhAT681+3XIKfn+98elNqOHxt4XhmpEXzo9f2U1zVZHU4Z+0sqmZ7QTUr54/1uQUVfOts3ejiMVHMzxjFnz85QVuHPZbO+u27R9lXWsvPb5qsxSFKuUCAvx+/vXUqBvj2i3votMnctSc+yicqNJBbZ6YM3NjLaFIbRvddNo7yulZe321918S2E1X84cPjLMtJZskk3+qOUGo4pcSE8e83XMy2gioe//C41eFQUNnIWwdOccecVEb44OLkmtSG0YLMUUxMjOSJj61dOquupZ3/t243qTFh/Nv1F1sWh1Le6qbpSVw3JZFH3jnK3pIaS2P506Z8Av38WDEv3dI4rKJJbRiJCPddNpaj5Q18cMS6Hb0fem0/p+pa+O2tU33yk5tSw01E+I8bJxMXEcy3XthNU1uHJXGcaWjlpdwSbpqeRHyEbw4xaFIbZp+fMobEqBD++JE13RKv7y7ltd1lfPOKTKalRlsSg1K+ICoskF8vu4QTZxr56ZuHLInh6c2FtHZ0ce+CsZa8vh1oUhtmgf5+fGXBOLbkV/Hc1kK3vnZJdRM/fG0/01NH8o1Fupu1UsNt3vhRfGXBOJ7fWsQ7B927afD+0lr++NFxrs5OICM+wq2vbSea1Nxgxbx0Fk6I499eP+C2LSs6uwwPvLiHri7Db2+d5nNlvUpZ5dtXZ5GdGMn3XtnL6foWt7xmRX0rX3k6l5iwIP7jpslueU270r90buDvJzx62zRSY8P42nM7h30+S1eX4SdvHGTbiSp+fMPFpMaGDevrKaX+T3CAP79bPpXG1g6+9uxO6lqGdwGG1o5O7n92BzVN7TyxIodR4cHD+np2p0nNTSJDAvnTXTl0dHbxladzaWwdnoHkjs4uvvvyXp76tIB7Lh2rq4YoZYHMhAh+s2wqe4pruG31lmFbH9IYww9f3c+Owmr+65ZLuHhM1LC8jifRpOZG4+LC+f2XpnO0vJ4HXtzt8k0GW9o7+fpzO3llZwkPXJXFjz6vG38qZZXrpiTyxIocjlc0sOzxzZQOwwLnT35SwEs7Svjm4kyum6LzT0GTmttdlhXHD6/L5u0D5fx24zGXPW9Dawf3PLWdvx8s58fXZ/PNxZma0JSy2KIJ8TyzcjYV9a3c8odPOV7R4LLn/uhoBT978yBLLh7NtxZnuux5PZ0mNQvcfWk6y3KSeXTjMd7ce/KCn6+6sY3bn9jC1hNV/GbZJXz5Ut8t51XKbmamx/CX++bQ2tHFssc3s7+09oKfM7+igVXP7yQrIYJfL9N1XHvTpGYBEeEnN05iRlo0335p9wX9kp+qbWHZHzdz6FQ9j98xg5um6xiaUnYzKSmKl+6fS0igP7et3sK2E1VDfq7a5nbufTqXAH8/nrgrRxdU6EOTmkWCA/x5/I4ZxIQFcd/TuXx8rGLQY2xHy+u5+fFPKatp5qm7Z3JVdsIwRauUulDj4sJ56f65xEUGc+earby1/+Sgl887cqqerz+3g6IzTfzh9umkxGhlc1/izDdVRJYAvwP8gT8ZY37e53gw8DQwAzgD3GqMKXAc+z6wEugEvmmMeft8r5WTk2Nyc3MHfyYean9pLXc9uY2qxjbGRIXwxRnJ3DwjmbTYEf22L65qYsO+k7yx9yT7SmuJDgvkqbtncUnKSDdHrpQaijMNraz48zb2l9Z173c2JZHPT0lkclJUv+PgNU1trN9Txku5JewrrSXQX/jZFyazLMe3VuAXkR3GmJwB2w2U1ETEHzgKXAWUANuB24wxB3u1+TowxRhzv4gsB75gjLlVRLKBvwCzgDHAu0CWMabzXK/na0kNuqsW3z1Uzku5Jd1XbAZmpcdwc04y101OpKa5nQ17T/LGvpPsKe5eLPWS5Cium5LI0qlJuo2MUh6mua2TN/ed5M29ZXx8rJKOLkNKTCjXTR7D56ckctHoCD7Oq+Tl3BLeOVhOW2cX2YmR3JKTzNKpScSMCLL6FNzOlUltLvBjY8znHF9/H8AY85+92rztaLNZRAKAU0Ac8GDvtr3bnev1fDGp9XaqtoX/3VXCy7kl5Fc2EhzgR6tjP7ZJSZFnf+m120Ep71DT1MbfD5bzxt6TfJJXSWeXISTQj5b2LmJGBLF06hhunpHs83PQnE1qzowwJgHFvb4uAWafq40xpkNEaoFYx/1b+jw2qZ9g7wPuA0hNTXUiJO81OiqEry/M4GuXj2dnUQ1/3VNGXEQw101OJH1U/12SSinPNTIsiGU5KSzLSaG6sY23D5xiT0kNl2fFc8VF8QQFaOnDYDiT1PqrFe17eXeuNs48FmPMamA1dF+pORGT1xMRZqRFMyNNV9ZXyldEjwhi+axUls/y7Q/3F8KZjwAlQO8RyWSg7FxtHN2PUUCVk49VSimlXMKZpLYdyBSRsSISBCwH1vdpsx5Y4bh9M/Ce6R6sWw8sF5FgERkLZALbXBO6Ukop9Y8G7H50jJGtAt6mu6T/SWPMARF5GMg1xqwH1gDPiEge3Vdoyx2PPSAiLwIHgQ7gG+erfFRKKaUuhFPz1NzJ16sflVJKfZaz1Y9aVqOUUspraFJTSinlNTSpKaWU8hqa1JRSSnkN2xWKiEgFUOiCpxoFVLrgeexAz8V+vOU8QM/FjrzlPMB155JmjIkbqJHtkpqriEiuM5UynkDPxX685TxAz8WOvOU8wP3not2PSimlvIYmNaWUUl7Dm5PaaqsDcCE9F/vxlvMAPRc78pbzADefi9eOqSmllPI93nylppRSysdoUlNKKeU1vCapiUiMiLwjIscc/39md00RmSoim0XkgIjsFZFbrYj1XERkiYgcEZE8EXmwn+PBIrLOcXyriKS7P8qBOXEeD4jIQcfPYKOIpFkRpzMGOpde7W4WESMiti3DduZcRGSZ42dzQESed3eMznDi9ytVRN4XkV2O37FrrYjTGSLypIicFpH95zguIvKo41z3ish0d8foDCfO43ZH/HtF5FMRuWTYgjHGeMU/4JfAg47bDwK/6KdNFpDpuD0GOAmMtDp2Rzz+wHFgHBAE7AGy+7T5OvC44/ZyYJ3VcQ/xPBYBYY7bX7PjeTh7Lo52EcBHwBYgx+q4L+DnkgnsAqIdX8dbHfcQz2M18DXH7WygwOq4z3M+lwHTgf3nOH4t8DdAgDnAVqtjHuJ5zOv1e3XNcJ6H11ypAUuBtY7ba4Eb+zYwxhw1xhxz3C4DTgMDzlB3k1lAnjEm3xjTBrxA9zn11vscXwYWi4i4MUZnDHgexpj3jTFNji+30L0juh058zMB+AndH6pa3BncIDlzLl8BHjPGVAMYY067OUZnOHMeBoh03I4CytwY36AYYz6iew/Kc1kKPG26bQFGikiie6Jz3kDnYYz5tOf3imF+z3tTUkswxpwEcPwff77GIjKL7k96x90QmzOSgOJeX5c47uu3jTGmA6gFYt0SnfOcOY/eVtL9SdSOBjwXEZkGpBhj3nBnYEPgzM8lC8gSkU9EZIuILHFbdM5z5jx+DNwhIiXABuCf3BPasBjs+8kTDOt7fsCdr+1ERN4FRvdz6F8H+TyJwDPACmNMlytic4H+rrj6zrdwpo3VnI5RRO4AcoDLhzWioTvvuYiIH/AI8GV3BXQBnPm5BNDdBbmQ7k/SH4vIJGNMzTDHNhjOnMdtwFPGmF+LyFzgGcd52OW9Phie8J53mogsojupzR+u1/CopGaMufJcx0SkXEQSjTEnHUmr364TEYkE3gR+6Lict4sSIKXX18l8ttukp02JiATQ3bVyvq4LKzhzHojIlXR/GLncGNPqptgGa6BziQAmAR84eoFHA+tF5AZjjN22b3f292uLMaYdOCEiR+hOctvdE6JTnDmPlcASAGPMZhEJoXtRXTt2pw7EqfeTJxCRKcCfgGuMMWeG63W8qftxPbDCcXsF8HrfBiISBLxKdx/1S26MzRnbgUwRGeuIcznd59Rb73O8GXjPOEZebWTA83B02f0RuMGm4zY9znsuxphaY8woY0y6MSad7rECOyY0cO736zW6i3gQkVF0d0fmuzXKgTlzHkXAYgARmQiEABVujdJ11gN3Oaog5wC1PcMsnkREUoH/Be40xhwd1hezumrGVf/oHlvaCBxz/B/juD8H+JPj9h1AO7C717+pVsfe6xyuBY7SPc73r477Hqb7DyV0vzlfAvKAbcA4q2Me4nm8C5T3+hmstzrmoZ5Ln7YfYNPqRyd/LgL8BjgI7AOWWx3zEM8jG/iE7srI3cDVVsd8nnP5C91V2O10X5WtBO4H7u/1M3nMca777Pr75cR5/Amo7vWezx2uWHSZLKWUUl7Dm7oflVJK+ThNakoppbyGJjWllFJeQ5OaUsC/weYAAADLSURBVEopr6FJTSmllNfQpKaUUspraFJTSinlNTSpKWVDInKXY++pPSLyjNXxKOUpdPK1UjYjIhfTvaTQpcaYShGJMcbYbY1PpWxJr9SUsp8rgJeNMZUAmtCUcp4mNaXsR/Dg7UWUspImNaXsZyOwTERiAUQkxuJ4lPIYOqamlA2JyArgu0AnsMsY82VrI1LKM2hSU0op5TW0+1EppZTX0KSmlFLKa2hSU0op5TU0qSmllPIamtSUUkp5DU1qSimlvIYmNaWUUl7j/wMV+6/WZzjgnwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(7,4))\n", "plt.sympy_function(bulk_free_energy_density.subs(A, 1), (-0.2, 1.2))\n", "plt.xlabel(\"c\")\n", "plt.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": [ " 2 2 \n", "c_C ⋅A⋅(2⋅c_C - 2) + 2⋅c_C⋅A⋅(1 - c_C) - D(κ*Diff(c_C, 0, -1)) - D(κ*Diff(c_C\n", "\n", " \n", ", 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": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAAcCAYAAABIx99TAAAABHNCSVQICAgIfAhkiAAADIhJREFUeJztnXuwVVUdxz8gmchVsshJGxPNqKthmIylTczRRDLG6SVRjeF20h7iCDnNiGLjrTEwCUVT/jCzm9lMzoBYRibdEsIwRaMJHyWZNxWLxxBahiB4++O39pzNvvux9uvsx/l9Zs7ce/Zznd93/X57r73W+m1QFEVRFEVRFEVRGs9s4M/Ay+bzEDC91BIpUVwBrEe02gbcC7y31BIpiqIoiqIoSo34OPAx4F3ABODbwGvAiWUWSgnlfuACpNEzEVgB/At4c5mFUhRFURRFUZQ6swP4ctmFUKzoAfYB55RdEKVruAPYCowpuyCK0iWcDAwBXyy7ICnQeKEo5ZM6hnzB7DgEXJhzoarEAcDngD1I70IUPcALiE3uKbhctnwYWA78E9ht/q5CeriayhGIBh+K2a4qer0F8aEVwN+AXcBLwIOIY44sr2gdoypapGEy8DpwWdkFaTDqI0oQK5BrWk/ZBUmAxovy0Dii+EkcQ44CdgL/obkNoInAf4G9wL+xazB8F7HH68BzxRXNmquQ8mwDfggsAG5F5stcV2K5iuYuYAPSeI2iKnp9xZTjReAnwELgdsTHhoBlwIjSStcZqqJFGlYhWo0uuyANRn1ECeIURP8ryy5IAjRelIfGEcVPohgyAhgAngEW0dwG0IHAccjTmoXAdqIn1h+P9BLdB/wJscu4gssYxQxThl8DhwSsf0Nni9MxFiGt+eNitstbL8cco5Vi3zOQ4Xr+p09vQxoDQ8CnM5St6lTNd5IwAWm03Vp2QRpOt/tIGA7p405TeAr4B/EPvKqAxoty0TgSjEN3xxHrGDIHceApQB/xDaAW8kT+BdpDsO4HPpGltCUwANwWsf63SKKEXmR87xAwtQPlCmIk8HfgFeCtCfdtUV+9FgNbkBvqOPLWy6GYAHKlOe73Qta3qK9eLlXynV5z/rt9yz8AbEaC5CTP8mvN9h+JOW6L+utUVeJ8BJprf4d84k6L+trnasQG00o4t8aL5tAN19owHLLHkRb1tc+wGBI0FrIXceAbgd9ZHHQJ8ADS6h4ArjffTwZOy1bejjMSOChk3eeB04GlSEtyo1n+/g6UK4jTgGOAXyLD96YDlyON11Mj9quzXjcB5yE6PBmzbdX0iuI183dvwLo66+VSNS3c827wLDsfWIOMGZ+M9FK5nIkk3PhDxDGboFOVifIRUPvHUXf7/N78LeOhicaL5tD0a22R1N0+sTFkFPAo8FfaY1f7CO8BWkB7TKU/00kPcHT6shbOtUjygPHIXKCFSK/X2QHbHoI86dkOHGaWTUN++11FFzSEr5nz34y8z2jI91nD8J6hOuu1FHkH0BlIV7b7CZrUVpReDvn3AI1CGgRBTzfrrJdLFX1nsTn/dKQ7/AbzfSnDh42OQS6WGwmnCTpVmSgfgebb3yFb3GmCfcYiv+GREs6t8aIZNP1aG4dD+jjSBPvExpBvIU8uvD0IfQQ3gE4y265H5tLUjX6k63o3kqpygPDu9esRG8z2LDvSLHu6uCJGstCcfy+wCelu7wFOAH5l1q32bF93vfwNPPfTF7BtUXo55N8AchMDrPQtr7teLlX0nQfM+Y9H5s/tJnyI7wSz7aqQ9U3RqcqE+Qh0h/0d0sedJtlnF/Lut06j8aIZNP1aG4dDujjSJPuExpBTkJtpf+awPoIbQD/Gbpyry2HI3I13Wm6fF8vIloryBKTb9HGGT57ajvQaHWpxnEHCb+KDPnfGHO86s90+4H2+daOB5816tzGrehWrV3+Ksl5q9n2K4S90bYJeeWgxSL5+AzJkdDcyfGUIeDhi21OJ7q2qi05JSOKDg+Svj5coH4H62N/WpoPkG3eaZJ/NhA+BdBlE40UVqFIMgeZfa/0Mkl8caZJ9AmPIKGTY25PAG33r+ghuAG1FXhxqm1N9EfCjgOUXA88CrwKPIcPSkhC3/4mmnGMTHtdlNfL7g8YNuk+Gplgc5zfAXxJ84lJYX2HOvSlk/W1m/RzzXfVKptdcpO57P/fQDhb+dUknAc42x3oCGcrnpwl6rSa7Fnn7zbG0g/462sNHw16oO8ms/1nI+rx0smUKcC+S4nUIODfFMfL0wbz18RLnI9AMP/GSd9zJwz5VqXM7kNdyRNEN8SJrXc2KTX2oSgyB7rjW+skzjlTBPnnVucAY8ibsW4pLkEQBQ+w/KTCKg5Ec7P4fPBN5QnwRknzhJuS9PO+wPK7t/o+x/xAcW87DziZzUxw7K58y514fst5NXz4P1SsvvRyzfyvl/i5zzXE2AocHrG+CXlX1nXPNee9GuvPPoa1FUIB3h+s9GLAuL52ScDZwDW3/T3ozWrQP5kWcj0Az/MQGh3RxJy/7VKHOjUR6jJ9JeO6sVC1eZK2reWBbH8qOIdAd11pbHJLHkarYJ486FxpDRiM9BkGfP5oTrjXfZ5rth5DuJBtmIC0v/8unHga+71u2CZnfYoPt/lfTzgBhy6FIir89wA8Its0AYocsT3TTMg6pUDsJHpd5H1K2z6J65aWXQ/YG0OW0A0rYe3DqrleVfcedzOl9D8RDZtmsgO1HIE/AtgWsy0untKS5GS3SB/PCxkeg/n5ii0O6uFNE/SyrzrmpqJcnPHdWqhYvstbVvImqD2XGEOiOa20SHJLHkarYx0vaOpcqhvQRPATOzabxmYB9JrD/mP8bkcmDXg5ExuLN8C2/BclgFkeS/T+K3IwleSvzEuT3fSdim8m0ny6UwZ3m/Nf4lk9FWro7kd49UL0gu14O2RpA3zD7P0rwfAYvddaryr7jJgg51rPsTLPsWYIfJiwz64NevJtVpywkvRkt2gfzIImPQL39xBaH9HEn7/pZVp27wJz7kgTnzoMqxYusdbUIoupDWTEEuudamwSHdHGkbPv4SVvnhsWQUSlO7jIP+DnwUyQn/hPIzfYk4CjgCM+245Gxe17GIYbb4lu+BQkwcSTZ/0UkXeWR2HWhT0S60Z5HMuOF8TiShKAXMfgui2PnyWXIy9jmI2MkH0HSEX7SlOsipBEEqheUq9f5tLMsrkUmZfoZpD0Zsa56HUy1tTgJeAm5eXEZQOYrtYCvIkHcy3LkCfA0ZCK0l6w6dZIifTAPkvoI1NdPOmXTsutnXvY5C6kXYXNriqJK8SJrXe00ZdR36J5rbV1iSCfrbaIYkqUBtBIJAPOQzCdnIZmdNjI8IByETHwKYsj3fUTAsihs9ndvrmxbzbcgtpkDvBKx3atIN957kAlYUdlhimAr0gC6Cmn0fBCZ4LUS6Vr0vohN9SpXr2PM3wMIn/eyhnZQrqteVdbi7cg48DUM/w3zka7z+cDt7D9RcjkSrGchv89LVp36kG77KE5n/5T2WSnCB/MgqY9Aff2kU+Rln6xksc9YZJL2L5AHK52iivGCgLIkrasufRQbe8qo79A919pOURX72JAohtg0gPoIftcKSOt6rcUxvC9B9C7bx/CsHIczvKUYdkzb/d0u0KBxuUHYZAlz6U2wbRHsQHqCbNIjql7Z9OonXbpriPajMOqoV5V9ZzPh8xvWRazbgwT6BcgTYf+E0Cw63Yw8WYviOYtj21CkD+ZBH8l9BOrpJ0noJ33cgWz2yUoe9pmF3FgtzrlscVQtXmStq36Kjj1lxBDonmttUvpJH0fKtE8SEsUQ27R2WdmAvETMyx4kY4M/Re5UJLjEkWT/iUjXWN7GbiqqV71QvYrnBuRmIGpYXxxBOm0nPu3r/zKc00u3a6p+Ek2QfbKS1T6jkVc+LMfuBqwqFBEvstZVP0XHnqrX9zRoDImmCPskoZIxZCLSAvRn4piJGOdC5GnwEiQ13tGWx7Xd/w4kG5Vih+pVL1SvzjAFGTIyJuX+YToloQcZe+2+b+Tr5v8k6Va7VVP1k2jC7FNmnetFnuaPtzxXlSgiXmStq3lgWx+qXt/ToDEkmqLsk6XOVSKGrCM4g8vFyIS03Ugr0T+ExkF+8PiQ48btPxp4GZkfo9ijetUL1asehOlkS4vg9yn1e7ZxUE3DSOMnDtH2jNsf6mPTIPu0iK5zDt1jn04TpEfWmJ6VFvExqMl66rU2miLs06LmdW4a8DT7p8yz4ZtI1om0CRsuAVal3LebUb3qhepVD9LqlATVNJw09s9qT6iPTdU+1aIsPbLSZD31WhuN2ieES0neVbseyUSSli8B786wfzejetUL1asepNEpCappNEntn9WeUC+bqn2qRRl6ZKXpeuq1Nhq1j6IoiqIoiqIoiqIoiqIoiqIoiqIoiqIoiqIoiqIoiqIoiqIoiqIomfk/e6vkWt5I/bUAAAAASUVORK5CYII=\n", "text/latex": [ "$\\displaystyle 4 {{c}_{(0,0)}}^{3} A - 6 {{c}_{(0,0)}}^{2} A + 2 {{c}_{(0,0)}} A - κ \\left({{c}_{(-1,0)}} - 2 {{c}_{(0,0)}} + {{c}_{(1,0)}}\\right) - κ \\left({{c}_{(0,-1)}} - 2 {{c}_{(0,0)}} + {{c}_{(0,1)}}\\right)$" ], "text/plain": [ " 3 2 \n", "4⋅c_C ⋅A - 6⋅c_C ⋅A + 2⋅c_C⋅A - κ⋅(c_W - 2⋅c_C + c_E) - κ⋅(c_S - 2⋅c_C + c_N)" ] }, "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": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAALIAAAAUCAYAAAApz2ebAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHwUlEQVRoBd2aiZFUNxBAWcoBrHEGSwYcGUAGUI7AkAEUEVA4A8jAhgwMEWDIADIwbAb4PVmtlTTS/Dl29rCqtLr6VLdarT979OPHjxvXvRwdHZ2gwwN0eXNIXeBzDP1bu/BAtq+74IkD30c0n/ehsSvvdXjI9YR19+SXDPf60mSE8Y1DVBS7Q/2L+p2qgg0f5t663s9vO4aGGzmlw5qbLS9PrPXREg9gnmVYZVeHhEP7Jc8HrW3a4yW+69bhqw4n62Aucg1ZXtXy0E92oDWgNLa+iPFPMD5IQfjPEH7IqdXYOsOoaJx9ywcIPJ4RQQ6j9BvkeE2rUxu9pwU4D+D9DPAY/PcV8HP6yqxuL6v5dV353YbO6TqgDdZ+A0Zd724Ae1CQHInvoFO5ZdSPefdHGX8+qAAj4ghwsNMDvwdUHXmvaDSTEbo65kq07+GB0ZmENaK+7dfrcYZLEbyej36moU4bRUf5Uc9Ff+gYBZ+FLLu04Jum7EvDoPB9xJ959+bJaO2QczdhesjyEOJfUWDfaDST0QigcZeKB8rIagSZRuQcaf4EJuDprhR5Whb5Qk9e385Rf2+BFzK/zII+T6mzqOseX/itcWhHXucQe9kiO4m3Sbne1hC8m+Gmjpzpue5jLvK9FZLQecekcI9ChhWgswmdPRz/bHbHHrwNCF/h675eekGOO9QnVKN8FA/vNFgE0Hm35+bIKHNCfU19RX1GNdeMB1+Sm7njvCacD6pSFta+sG5qUBc3r85f67VZ39RCJx0VH3TSCydZRzuccxqVkVdjTqOx61Qfk01hzr1b93ZQrumboCF2oAHy6cCfIH+P6g32mbF2196WTYJLAsx46hw1aKT1jf9gvL1zZJhpfI1S8kb6PvCa/JjxK/lREnzNe2FtJbcFXvqLuRgwOlSCo9X5FcCHStE7z6c8lv40P+5wlElaRedu3Txymhuz5iH4UuPYd861fj7GrHmgP8V42xZc92DnHBncsHWzh8pBCZsv0gdWu6hr+YqU8Ye6sabe0tc+7m34kgHxZO+IzEkywkn8OcrUJ9G+3z69Dv0WquAf7VOMKAV2sPYtQZ39MQIW+DwtvR7uDOOs58ZHdPVrg0XcVLL8OnWSk8kaPkMNm2lUzvrUNEcEVvhkWep9GuGp8+xWGcGf21zWS1u/ZL9iL2v6sa9lDZxP1FGU1SnfQcdULYo2aL4GyVMazBsUHlL9kvRUOOal8YJx8g03RcBtajmN4Hk6PI5N9GHcRJZ6PcOXU9utGdnLGkIyTM5VTm6eE27xmyUwzVcNxhIs9Lu+hmjW5TWrwA6jMvNL0dg9l0+vU9wYzV7W/MHRKZiay+UaxYg/sqkyu3ejtaUvOjqxxIfyZbrla4ZwwveyMjf0mQFc2GNoZ+ikKJ307ZG3HUPMTWl+kGAchloRgLWpsVgzSq1sFHPF8UI+5nZ1ZPGSc9MqS0kN6Lsx8i8HNfiNWuBCl+IAzLn5zeHpcSu8xiHEow6v1qDB+kaOHPB9m3mv7GcPNxqD694M5WM+bF5oZ16j9GnFZyb8mmDYw0DfvUj8btLZuRDaFd5arpJMTIf0KMaVnqdT47XglXJaT+b+7HNd/ARaoyxesV5LIBh56uI15HWl3LeQo05Z5H/KXK9PjV/6wI2+YJhyRNpRYLvOTM+SbmT5OrQ0vMXfWuYRzLnPVfL8PSFuELD441N6+NP3VvCHEh/3Op0p5sxnXC4FuJT7MtGkGgXgv477kP4tYS9HrohG7htTv9JJzoBA5TNVVkJj/SEgYxWti47XOJH4zCX4GpC+SmjUdaU4RgWkIe5RfQD2/5sh/MxQFYmmG07rq135PcCjQ1ojyadxxowrfuxlOEaNZ1+YJfo9znmOjZKjYoD6Xd2p/nbgWDnNp/3unOxKG7IP6WR7S1/89xW8c01xLdZvNitbDjIRo24yoOgIogE8fWGo+8BFXwfS0O+y4XplhKtpHTMWv3FuaVCcW/rw3j9AxZOndJvbAnmU2XkfEBsXdQFYuT1wpgbh2HRXC3zkoY5FzwwVeH7KqvevJ6LOjew9wCHG6Hma+d7u6SOvubPOGzoEiHqMbGcA8VYqxX2hun/hK+7PCLfgNB2Y61g7V4hpGBXxKrAqvELEXJMn53mdffjpLK8bqacwyksxqq3kX3lN3uZhAuq4dd4mXj3WAXXeGl78Jn9dt0fASkNea3PjLFuBBV49077lNftpbsaPdVOljXL4EQ1w5V/0H8HM5sDTru6nNOJA6nwrejOnH5SHX0+TNfUsdqbfyMRYPus+Q8q/+NDODtwLdhnjrGx5rF2GDMETWTToovNnuOGDKWjNWulTh4d3htPPg7+zIwetTCOccLj/wHgo0yOYNjl94G/SguN+Nh8RajzWGie/7o5s1F6JBrXCV62PvEb+JvpsKqPG2xU3eICvUy0euIDftYWHjpj0pG2cblOa4BmVS9QVj6L87kNzgI7yogDXspBXmRb4mIjc6srqgaxxNfu/H5vnf2hkDknzAbyld8GV0B95TS1eUD9SfbRtpW8oAZ3ktIy/UU+p/1DfQM9+Kf8HR9bAXmHN46FoeIU6GMX83Btk5cG0JCa41+bALulyiPVr78huCkY20vmo7D+nHWLPLpwm+plCGdWu/K1z4ZuTGf4LWdnvAy60a7EAAAAASUVORK5CYII=\n", "text/latex": [ "$\\displaystyle - div(M \\nabla \\mu) + \\partial_t c_{C}$" ], "text/plain": [ "-Diffusion(μ_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": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvQAAAAYCAYAAACVxWnGAAAABHNCSVQICAgIfAhkiAAACedJREFUeJztnXusHFUZwH+XFi0PbdEG1Ki9cmMRhXopSFGxLpCAMSFiotQQHxtTSKhYE9LE8lDbJioJsdRHE1FQQBMlEWI0ilYeKghcLaXyTptqJVisYOVpqYr1j+8MO3fuzOycOWdmzux8v2Szu3Nmvp0789vvnjNzzllQFEVRFEVRFEXpEK8HvgPsAvYBO4ENwGEVx/oQ8HXgduAZYD/w/RKf2QV+hRyf/cDKnPWujq13dQ37laQpl0B9cqUNjvn0yzXexxgch+UlP78ttMEN0PzTFOqH1nWyaIsbI8EEsBs5iD8GLgNuNe8fAV5dYaytpuxZ4GG6Jbkt/wD+gxyjqzLWWQL8D/ivWe/8enbtJZp0CdQnV0J3zKdfrvHeADyFuNaFCn3oboDmnyZRP7Suk0Ub3AiGPnIAeiW3/6XZ/tOJ5evN8m9WGOsU4M3AGLL/XZLchgnk2PwO+Bvw+5R1DgA2I0nmbrP+Oyw/p097XQL1yYW6HHPBp18u8caAm4EdwOWMfoVe80+xWF3NP+rH8Fjqhmc3esD1wGPI7ZHHkZNyloedbpI+5SU/0mz7Z+SgxnkF8BzwPHBIDbF6FJN8pVnvkyllc5FW3q0F9rdNLEP+5q8BNyHHMXmMzzfrfAJ4Avg38HLLz+kzGi5BMZ+66FIWdTlWFp9+ucb7DOLGUmAN+RX6UXBM8499rB7d+X+mftjF6qFuxMl1I7kySL+m24BTkSsr683744F3ed39dnGqed6EyBHnWaRVdTBwUs2x8lhsnu/JKBvLKGszJ5jne4AtyHGciJXPB74I3AX8xry/H2m41oW61G5Cd8y3E2XjHY3cXv8q8NsCnzMKjoXuBrQz/4D6URdt9EPdYGaF/kvIFZUbgHGklXcRcI55v9Hr7reLo8zztozy7eZ5Yc2x8lgMvAA8mFJ2vHne4vgZoRF9KTYD95rXi2LllyEt9k8xOAab69m1l1CX2k3ojvl2oky82cD3gEeBiwt+zig4Frob0M78A+pHXbTRD3WD6RX644DPmhXOQS73x3kO+It5fRjSh2eCevkRcGHNnxkx1zw/nVEeLZ9Xc6ws5iBXyO5DBk8kiYRItlqbOrc2ZHkwhnj8L2SwTfQFjr4US5BbclciX5imEqa6FD5tdsy3E2XifR45Tn1gb4HPaJNjbXYD2pd/oJwfodVT1A+3WFlo7jDEK/QXmverkX45eVwM/AIZ6BRnBdJX6gXk4L1nSJwkw7ZfB1zKQJI8djKY1id6fNeU3ZZSdo3lviYZM8/7HeP4irUIuUqWdZvpBOSW1/bE8qxzW5SlwE+R6an2I1NQ2VLWg4Vm2VbgReBPyOwaixC3NyKjyC8168dbxHnsRF3y5ZJrjvBBEUdDcCzv4TJIzKdfafFORM79V5Dbw0UIyTHNP0Io+QfK+VFFPaUtuaNLfvjKHa2vu8yOvT4D+CciQB4HA+cCZyaWL0P6Sq4A7kA6798EvBW57TqMItvfhxysjzK8+88GZrb6JoEPANciX4I4W4fEi1qSWY2JVybWqytWFnl9yuYirdLbmf5Fyjq3NhwC/BFJKDeU2N7FgzTJtwLHAuchrdpzgT2mbDHS/+yBIfukLgmuLrnmCF8UcbRux3Yg/wSKsiunzLcTNvGirjbbgM8VjA/hOKb5J7z8A/Z+VFVPCTF3dN0PX7ljZOouc5A/9t60wgQfNkHHEsungG8nlm0Hvlwgps32X0AGUpShT/mR38vNtldmlEdTM51WQ6wew6/SfcusM5lS9n5TdkVieda5LUuZVq6LB1eYz/x4bNl6ZDDOHmQaqOhve5NZd8py/yL6jIZLMNwnXy655ogqyHO0acfK4tMv23jzKH6XYUMsRiiOaf4JL/+AvR9V1VPihJw7+oyGHz2aqeu0su4SdbmJVjy8wE6fjLSE4q2dlyEtiE2JdTdRbGYcm+2nkHk3DyoQ1yfRnYvTSZ9+6d1IX9G7a46VRdRq/WtK2UfMc7JFm3Zu68TVg/go8YgtiN/zgAsY/G1N9U+EbrrkmiOaoK2O+XbCJt4+5JcL0x7RBaM7zPt4d5wQHNP8E2b+AXs/qqin2KB++I+VRQh1nSByR3SA9yKX7l8HnJ2yswuBWeb1ODNv98435bsTy3cDr0mJl8Rm+13AgWZf62QHcnLGkVHGcdYit2uuY/pg4gngLcj+usay4UDgGPP65ETZ2cjtHpApj+KMk38rv2pcPDgAaaE/j/y6XMTPgQ8i02fFf6ihaP/EKuiiS645ogna6lhZJ3w4the5Kpf2+InZ5lrz/nrzPhTHNP/Y+1F1/oFyfozjv55iQ9f8aPv/J1eCyB3xPvSrkYT7Q2TS+geRlsEk8tPdrzXrzSG7r2eytTOWsiyPIttHMybUfYUepG/UncjE/6chB38J8stm24BLEuvfAixAbpHsdIx1FoMf9ooEeSeDAS5PAqvM62OQHxvYDfwAuBEZUPF2ZCT134EjkJkoLmfQOk47t2uQW0R5nAL8esg6NpTx4GjgUORWVny+2z3Iz00nafIKCDTrEhT3yadL4J4jItZQvZdtdqyME74dK0pojmn+sfOjyvwD5fyosp5ShK75MQp1HR80mjviFfqfIf2VViMH73TkwN2PdPSPeBKZ7ofEsheZ2RI5nJktljRstn+VeX6iQFzf7EBaR+uA9yH9sx5HRF3LYMBCFbEmkYZWnCPNA2RK0Ujy6BbUOlPeR+SZQr5E7zVlC5g+ECjt3H4DaeTl4WtAo4sHWVNTZbEY+WI9ZLmPvmjSJSjuky+XXHNEkjq8bLNjPv2qIl6cUBzT/BNe/oFyflRRT7FB/Sgfq6m6jgutzR2rSB9ZO4UMToizDbvBSEW2X056XyllwEakVXii5XZZ57YsZQeWqAfh4NMl1xxRBXmOqmP1EJJjmn/Co4wfVdVT4mjuaJ6q6jqdqbsci7RE5ieWL0Pmr1+O3ELYgPwY1YKCcYtufx0yoErJ5i7kBxZsuyVlnVsbDkVa2JPIl2KVef3GgturB2Hh0yXXHOGLoo6qY/UQkmOaf8KjjB9V1VM0d4SFz9zR2brLncio2yQrkP5T+5DbB0sT5X3kQI1nxB22/UHAM8BJ1nvcHWYhgyvKXmnPOrdF6ZE+Rd01sXX6qAdtoAqXXHOED3oMd1Qdq4cQHdP8Ew4uflRRT+mhuSMUfOeOHh2tu5yB3EqYNWzFBGuRwbazh62YwQXMnBZImc7bGD5vax5lz60N6kE7aMIlVzd8oY7VQxsdUzfqw8UPraeMNpo7PLIS+1vlf0AGKpTlPOAoh+2VYpQ5tzaoB93B1iVXN3yhjrWHuh1TN9qD1lOUPDR3KIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKEoL+T8EqW4OLF8VyAAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle {{c}_{(0,0)}} + 0.01 {{μ}_{(-1,0)}} M + 0.01 {{μ}_{(0,-1)}} M - 0.04 {{μ}_{(0,0)}} M + 0.01 {{μ}_{(0,1)}} M + 0.01 {{μ}_{(1,0)}} M$" ], "text/plain": [ "c_C + 0.01⋅μ_W⋅M + 0.01⋅μ_S⋅M - 0.04⋅μ_C⋅M + 0.01⋅μ_N⋅M + 0.01⋅μ_E⋅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 = ps.plot.scalar_field_animation(timeloop, rescale=True, frames=600)\n", " result = ps.jupyter.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", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" } }, "nbformat": 4, "nbformat_minor": 2 }