test_entropic_model.ipynb 102 KB
Newer Older
1
2
{
 "cells": [
3
4
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
5
   "execution_count": 1,
6
   "metadata": {},
Markus Holzer's avatar
Markus Holzer committed
7
8
9
10
11
12
13
14
15
16
17
18
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<module 'scipy.ndimage' from '/home/markus/miniconda3/envs/pystencils/lib/python3.8/site-packages/scipy/ndimage/__init__.py'>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
19
20
21
22
23
   "source": [
    "import pytest\n",
    "pytest.importorskip('scipy.ndimage')"
   ]
  },
24
25
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
26
   "execution_count": 2,
27
28
29
30
31
32
33
34
   "metadata": {},
   "outputs": [],
   "source": [
    "from lbmpy.session import *\n",
    "from lbmpy.phasefield.analytical import *\n",
    "from lbmpy.phasefield.eos import *\n",
    "from lbmpy.phasefield.high_density_ratio_model import *\n",
    "from pystencils.fd.derivative import replace_generic_laplacian\n",
35
36
    "from pystencils.fd.spatial import discretize_spatial, fd_stencils_standard, fd_stencils_isotropic\n",
    "from lbmpy.phasefield.fd_stencils import fd_stencils_isotropic_high_density_code\n",
37
    "from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method\n",
Markus Holzer's avatar
Markus Holzer committed
38
    "from lbmpy.forcemodels import *\n",
39
40
41
42
    "from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments\n",
    "from scipy.ndimage.filters import gaussian_filter"
   ]
  },
Markus Holzer's avatar
Markus Holzer committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAToAAAArCAYAAADov0TWAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJYklEQVR4Ae2d25HVOBCGzdQEwEIGkAGXCBgygN0IdskAat94oyADIIIpyACIgEsGkAGzkwH7fx63S9bx5fhYln1wd5VGlixL3b+O2t2SrLn269evIic9f/78ltp7p/huzna9LUfAEdguAicLiP5Obd5ZoF1v0hFwBDaKQFZFJyvukXB2JbfRH5uL7QgshUBWRSchHyq8R9jKheXSyRFwBByBWRHIpuik2F5KEsKPSqLrs0rmlTsCjoAjUCGQRdFV1ttPxSi5n1Xbt7wXHAFHwBHIgUAWRSdBnknJvaoEMovuRpuAKve6Ld/yUJoK3xV8rs9A8dgRcAR6EZhd0UkhnYmDDwEXpuhuB3nlpcqyWBGWjYtYGmuQep0cAUfAERhEYHZFJw4eS4GVCxAVN5dVjLKK6WFUNr5f6D6KEuvQ6tkp4xmOgCPgCIQInIaJ1NdSSk9VJ65m6I6ay9pYjFAZLDT22NWkPMqwuZjV2pC+K/ExzPBrR8ARcAS6EJhN0VVKqmhRUuTxOUZs0WHNPYsYvad0m+V2W2XfRGU96Qg4Ao5AKwKzKTq19lLK6Elrq1dbTMyysyKx4iOfBYdzK0CsOtvKhUX82hFwBByBBgInjVSihJQRbujXnuqw0hquK2k9V6+kVtfsu6vdVuWVrqzyQldYSSdHwBFwBLoRuJb6o34pIyyur4r/6GpW91hZRRnigparsIpRav8o4JJSB/cfKKAwmY+j3J8KbFVxt1VAODkCjsB+CCRTdFI+Zm2hoKBvCg+Uj/VWkq7ZPoI7G5b5qHyUV/g8z7Ja+0MB5YcShF4obfvxrnJW+Fc8oqj9hJYV9o2ztE0Ekim6bcLXLrUUHVboHcXX2kvkyxUPvCgeK9jLJdzqAyMo5QsFXja8YJyOAAH1FdM3ZkRgANj1EXB/OIuHyn16eJP+ZBsC6ohVndAifnDz3yj+T/EXxSi9BimPbT1MN9xVcGXXQGd9CfURHg6LeSi7TwqfFeIXmLJ+L5oi98nvBcV+0giwMwUsnTkoyQktKXlUXVhtTA10fXViiztdq+TJcUopX3Lm1l8hv91zBZTdF4XV7SmdqX8PlnurFh2DnpCU1Lm8aQmmMKa0kZJHc1u7BgSKEGIjdi5KKV8unhdvBwUiJsCOue1L4sWZamcgaf9OlXuTFl17v0zLVUegLNZ6Qku5RUc8drmlKGYGzW/v/kzr5VU8TV+ySEd/bYkmyb1Vi26OHwiT+WbJ2cEFuBZrIKyAnTe/+OWt+1aBmBVy47vQ9WsFk0e3m6R7KHZcYVbHuxRo86GEKbUJz/8q2LFfHBJBHzQUQMUnVjbzWDcVXijgAmERlXwrXrWs4rUQj/DMYKcvL5RmXpW4tY+UT/9sSu4+mV3R6dcwlQQwP75w/ssUxs4JLVPbGvu8eGMTNkqh0DXfHhvd1wWDgUHeWKBQeuwpMlkVnfgDb+YVayWrPORE8dWfEQbl+Lyw7BPFKIhHissVccWrllW8liQ+bVGJzyd7t1mprOGzGbmHZD6pcPRoGgJjTmiZ1tL4p/nRQ0/0Y3gVBJQb1g2f6mEthLTaU2TEq1mS8B0q2MbLRvdQ7ig1ytmLx2QMn1utrMasxZXsJHes86DM5uTep69Pq0IsUZdvfQNsIK7fpFZO9eT9v4lVw2q3c6+a7vHWt4FurBKXLqXut5n935TfsHDCB+Nrld37hJb4WdJ6fm4ecXcu1U482Gn7vQKWHVYdW0+QHbxQEDUprxw8iqkrpMFTZPRMavmoD+KTQbNQcUk/KB0qANw2ysRf0WD5lXORupdMVtUFRpPHkeroI3gv1FaoqOPyueVO3b+xPKSH5B6UGUV3qYom/49V1dOpcNo4z5EnntoUWaF83BWOj5r0lYWeL18OimMFUCiv7YSWHbFVblYe1SCDuW+RwVZaKccAwsKp3T+loXsK/E5iGjxFZgb5SnlU71Df8blgqPgKPYM1aBYhsiSTVXUnGUcw1UO8lHZeWFH53HLP/ftFvCG5B2U+jUDy5DgExp7QMq72iaU1+Mo3oaoJ5w/jWk1Jm5WAIoiJes7DzEpphFk5rz/3NSbeeAERYrlRkoXumwI8Bllh2Yh+sH6yvDreotz7ynxSo+QXoxAQwAyasSe0jGojQeFyYKseG9iNKiUDAwfrFpfVyqz9FBksGlzVHZIMyBJSbP0wJREqirXLGsrCNZZ1r5KvHtii3L0yu0VX/TLGRBpQWAJ8tN95QovuX1AnZRXiTuBWDvpLjVy2ta88lCBzcQx8TokxIv1J95nbQk7KMbXBJ2KkkQVXga0cS8jFfEzDXRIfWG/kEwqlkRnFDb8lKc2CC7KELu/aZb1iXn/FP7IgJzy3kspsTu59ZXZF1/qTac8UqPzQUA4MmEJpLLqhE1pQiOUJLTyTg9QePMJr6bpW6bBpBs2Fwt+6F8/fvaiee6qYQcX3r2xQRbmgSLAq2N4QT/Ire36iXQUsMXix+UXSDeWne1hvb5WPHJAp5dCdXbWsV2zXf60vzfKub0QXW5R7UOZFTi/Rj4+BZhPeDBwGXbbTM9Q+Ls7kxQjVMRsdA49ThM8tn9rDomN1eXWLZvvgKL5Z3byhmEG9Ny0lt9pNMsYOkbtN5pO9EUtUUEyg5Mqd6LpmbxduEZYDrlFpKSVqqq+aS90krJmOgccp+OWWj0UXfmdHQxoPbGzGOocwCM7Lq3F/lpL74P5NIPeOzEu4rrgcDTdDgmHN8calU/vmvcZ1cUdptTVk/nc8mS/7GHicgsYC8vESXcTdnoATY4Vpg+vUoTieZtin6kXknti/U+XekTm7RaeegYnv1nlBT6F8mGvB4nNyBJIgoN8TW4CYS0VZYCExiI6F4JU5ReKxLuvm5O7r6+xzdGIGq43TdxvfgVb5+PV++KNAcHIEHIF0CGRXdF2sS9GxgsYCwVFOFnfJ5fmOgCOwPAJLuK47Uku5sXSOy2orsTtlPMMRcAQcgUMRWIWiE/O4s3xgHm7mPFQmf84RcAQcgQYCi7uuUm7sD8JltW8uGwx6whFwBByBqQgsatFJubGlhE2QruSm9qQ/7wg4Ap0ILKbopNxYYeWYn3rZXNdYdszVOTkCjoAjkAyBRRSdlBmLD/cVx4sPKD8+B3NyBBwBRyAZAtnn6CqLjU2QbV8nnOl+Y39dMkm9IkfAEdgsAqcLSI6Swz1lfi6mo/oWMWbe046AI7BOBP4HJQMqc3w9Y1oAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle - \\frac{A \\omega}{2} + A + B \\omega + eq \\omega - fq \\omega + fq$"
      ],
      "text/plain": [
       "  A⋅ω                             \n",
       "- ─── + A + B⋅ω + eq⋅ω - fq⋅ω + fq\n",
       "   2                              "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A, B, C, D, F = sp.symbols(\"A, B, C, D, F\")\n",
    "\n",
    "tau = sp.Symbol(\"tau\")\n",
    "omega = sp.Symbol(\"omega\")\n",
    "eq = sp.Symbol(\"eq\")\n",
    "fq = sp.Symbol(\"fq\")\n",
    "\n",
    "test = fq - omega * (fq - (eq + ((1/omega) - sp.Rational(1, 2)) * A + B))\n",
    "test.expand()"
   ]
  },
77
78
79
80
81
82
83
84
85
86
87
88
89
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementation of high density difference model\n",
    "\n",
    "According to *\"Ternary free-energy entropic lattice Boltzmann model with high density ratio\" by Wöhrwag, Semprebon, Mazloomi, Karlin and Kusumaatmaja*\n",
    "\n",
    "Up front we define all necessary parameters in one place:"
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
90
   "execution_count": 4,
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 0.037\n",
    "b = 0.2\n",
    "reduced_temperature = 0.61\n",
    "gas_constant = 1\n",
    "κ = (0.01, 1, 1)\n",
    "λ = (0.6, 1, 1)\n",
    "χ = 5\n",
    "φ_relaxation_rate = 1\n",
    "ρ_relaxation_rate = 1\n",
    "external_force = (0, 0)\n",
    "clipping = False\n",
    "\n",
    "domain_size = (100, 100)\n",
    "stencil = get_stencil('D2Q9', ordering='uk')\n",
    "\n",
    "fd_discretization = fd_stencils_isotropic_high_density_code\n",
    "target = 'cpu'\n",
    "threads = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Free Energy\n",
    "\n",
    "The free energy of this model contains a term that is derived from an equation of state. The equation of state and its parametrization determines the density of the liquid and gas phase. \n",
    "\n",
    "Here we use the Carnahan-Starling EOS, with the parametrization from the paper"
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
127
   "execution_count": 5,
128
129
130
131
   "metadata": {},
   "outputs": [
    {
     "data": {
132
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAA4CAYAAAAcn+fJAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAYP0lEQVR4Ae2d4bXcNBOGN/fcAkJSAaSDcFNBoAMgFSR0ACf/+JdDOgAqSKCDkAoS0gF8FRBuB/d7H1+PI8uyLdvaXe965hyt7ZFGGr2SrPFI9t65ubnZOTkCjsDpIvDTTz89k/b/6Pjn6dZiO5qrnb5Qbb9S+KjA+SOFF+J/0NHpiAioDe6q+Oe1CrQN9FT869tT/z1VBNSGvyh8P6T/5VCkxzkCjsC6EdAA/0Yafqnjr+vW1LULEPhZ5xiwVZvp+IOu3yp8FqTx0+Mg8LPao5k0df6L1PhL4cFx1PFSCyJA275R+Lovz4u+COc7Ao7AuhHQwK6eZnVsbuDr1ti1qxF4oeOrAI37Osdb5HR8BJ5pPOG9M8J4/UK8h8bw42kioDb8R5r/riMPIElygygJizMPgUDqJiPeXQVzVe84V+AmxeS/WpJ+nRsmOis0ddmD8r8rT55gnfaIQOm2VX4fCIHKePmYeJ0mIlC6bVQ8DxfvJ6pxFsn3gOXqcFEd8cp+r2NyPvEls9U12XSF1Lh2M/1X0rh2cQ1iDWfRHHnJYADQsVreCV1jAPxYF3ylI0++P4ofTgB19O6t+HRMi7NO+qUl0JFymPRZ/9WhQ9fiV0sNdV62/k9C8uOJoLO3Rjz0NN1JxzW4pdLm4JtTFxVRhmr9r3Tsdf8uLUl559S7t5hc+dLpUEh50m8wGFlOvIYXknil+2mYfda5dOBJ9YkCffvoS57SYTXtPaF9io67RDtwj2B50+5RWW1bOtGB2qYolqUxyMlPOA2O+zoP5pPfFL6N87yIGX59WgioA7C+/U5HjI6XOscYYZ2UG/4oLZBnsrkXFlCXWW1c0znGEoYNN5K/dB66oU3so06YrOjEGCV/KMQTGBtOMVKYMOKA0fdUwQiDBhwscDNDD57AG9I1ZVkajqRjkIBbnDYX35y6NDoUOKGd9zaJCofceierkitfMp3ywiOHAWw3vOQYUDz80v00icMQU3owXh8rfK3zXjf+UB6l4lT+atp7YvvsbdxJD+5L3A/CB7RSkGfnc6i2kUJ7wzK7sjMSCp+scR9kzX3zG+QCXnXqHqIYkRO6VoPydhGdAUOiIp3jMeGaSWHQezBXXnJ9N2+eMM3rYvpgcKAnBlS8aZSlg46VXgkGP0rTqYd41WSnY1V3HSmDm3pMGExsWG0w0jlpWYbDADK+eYaeW9o6z1x8s+qivEsRdejgUiLzifXuFJkrv4d0GNdVf1Le9FEmtBTtpZ+mChrjSU/GK8YtDw1/KnwYk0nFS44HDpaXJxvJkqEv5fbzTvG58rnpVMCU9tnLuJOu3F/QI35A69R/jKG8Vt82dR32guUYPsQvxCh33FeqqCzGHOPsO4XWeLmoUvjPqSLAzT91A30n/ldq9I4FHFV0srzyZJKhAxJiYuD/nSgXY4MbbmXExEIj19QlRXiDQuOLpcJcAwHMWnVQXqn6TMYnpWhpnnS1if596bzr/JbWO1e+dLpcOPbRT7PKVtsxDv4L2hA563ssMc+luxIkzKHcdujLO1c+N93R2ocKqm24T/Egh+eOyRNDc869y/A6hbYxXY91XILRHJ25d9IfW3TRuvKLU0OAG8fHhNK2f4j4IZoj/0Q3h5ZVHRSA4cN6+3XAC08n37CVl3lwmnzE48ntRcO4PcFwwg2K1ycsh7R4yxpSPE/in3E0ps5tqSxMOwcfy3KfR/Qawnlp2UvrnStfOl1uvYv309yC1c8YG4xZG6OIggP811wcgXLboU+1XPncdMdsHwwf7gEsqT4k6BwPXuo+K/beKRezPkWWyvfle+r8v1WBzgPI5anXaqv6a6CGk34fDPf6IubIS4ZliNBgaGWv+I7FXSeoPBqKxzPTIvEqd72YvHrMzehFKp0JKY403Ki4STWk6z8IYmDY8AROPF6jcFlMl11SWm4alZte55Wxp+NkfCUzqS5oIhnKBjcmROrGk2k4WYrVoUfikL44qezJ9Q6VyJUvnS7UYexcZR+inw61K+U/lx7/1rrSnouXZsbqnYrPbYeULLxc+dx0dZ6T2kd5Txp3Sj/UNiy7MwZay++SCb3RqLl3moJZSpk58pIpiWVKrbXwmIvw1hKae+nlWrRzPSYjYMZO05iJHIYmt0ny6jRM1riPxybrlhpKjzFUTfStiNsL9HttHbIugyU3XNWN9yaSw3AhdEgy3ypgsDGoSYOuybTi75QW3bg5MiExQMIlqEn4SHZyXVQ+ujEgq5utjujCjTjeayVWiyhrX0+sU+vdUkwXufKl08V6TLqu+0KRfjrWroqnrxHWQLnt0KdrrnxuumQ5A+0zadxltM3Y2EvqtyfmIsyk01T5oljuCZNS2dr9kzHfjMWLUrl7PqtEAK/LEgrleVurb6lsqAw2U+O94Y2aFolXrdEbU9cYMBhCSS+U4um8X+mYNJbExzt0rYBniDSkx8Cy5TBdfiLx2UT4UoEn0lcKbGxNpv0k1Tpr8JHc1LpQDkt84ZMnBhkGEobREHGjo57HoqbeMxXIlS+dbkjdIv207j9z23VIv2PG5bZDn4658kPpku0zZdx52ySbp8F8Y1ja/dOMxgoc9xAl+8hJMM3CTSlrjWwu+VSabHkNFDwuSSMllbHxJIcMe136XOCWNDxiFGH0sJEx9kZhPMS8SlZp0ZFlBzMwMFAwOriR/qZz9g3ZIKhkwh/FYbQRz2vbPCVm4xPmE50P1YXvYLSW/XSNAQfdvT30/hLfq5/0J54368byCQvAu/ZhKF/FlexXvfpH5eSmC+sy6Vz1LtlPl7TrqN61rimDuWobxVv/D/MaentoKb658rnpQr2r8xnt0zfuvG0+oZszlkldFMs99N9PNZp2luyPbhBNA3E1qdWxWL5Cn9SkZ7yk8YBQrrzSMUnjtejNi/xiUnoMlHs6Jt/8Ev9NHf9lLFtfWx3CaAycPj1Yfvo8TKwyMHLwFrGBjkmk2qAt3kPS6fiBY0B4aEiHQYYsUSk9jFfponST6qL0lEEe8SZa+FCs1y330++1Tu2G9olbnyl/4vtw7aQPGcgqwLI6htHG62uDXa586XShklPOpUexfqq8lrbrqOoqI2Xw7MRnbPAQ0fHEDmWq9Ktq71hX6dfbPorLHndK623TBrc1lg+Fpcop2n/bVZp0ZffPlmF0OSmLjSZWI9J5ntfVt6f4p+JfHxkSWxaK1bDGTi4tBYlz5K+U/pHqiqclJIwKbsDw8QI13g6dc3N+oGPjGdJ5hZuONpmSb6sz6hqqdFe6llGga9qAPFp8BOo4jLZOe4iHbhhChgki1YZJ8XnTrCNDgppy8CHppLooPbjw1B6X/UR89DWMdHkUyq13n3K58qXT9emT5Avn0v107e2axEHM3HZYKj+pnIz2mTLuvG3arWf3Q9oE2gKWtzW9/WU+gVr3YDeIbkEZ+21980YDFRc7kyreh2MSxkhq0zDegdSEG+uaI8+AsUHTyAuD/+Dr2Bg9ROoaQwkDqjGQ4IuYfMI9SL8m0pCOJ7lOeeIxYKGOEaV8eMolpJbZkKHzh3kyCFJLaFaGpc3Bh/yn1qVTR+kOboQczw4Gkw1onRan3Hr3FZwrXzpdnz4dfo136X66tF07eh6IkdsOferkyuem22W2z5Rx523Tbr14jtgCliECZhC25pOLMIWf9yLwTAOUAWWEEcLkywR2NFL5GBgfdcTYqEjnTJR8gZMvNFcET+FGofKMBPwseUsfHSmnNSkrfzw43PQoj+94NEE8NmVf62hUxdkFR8XzWj/UMrJuWU1ZYR51VHVAhv0/sU643Hn1PvS6YKzhbm9I8WCIbKOneLn4ZNel1g+cbEDuah64UXbHAyZ+TO/EMOMtjlt8LR2y6o3eCrP71YRysvSJKn6/vm5wtniVW7yfgoXyX9qupuJBjxPa4VDtnds+WePO22Z8jlCHOxcse8d9NKjoY9fqG6355DJK5JdpBFj3ZH/JGglLHw/WIx3ZRM3xsa6biZVGV8AgSNVhVF5yDSkfvGN0Joi9Nkzk/JcaexcwMojDCImp0YcIpWdpCL3JD2Liwlr/XLxWJyVShDx8jIEOSQaPD0YgG6jJxwgjyTw+FU/XPA2hu5UNH73ZiN1KK94oPpKZUhczrPnekBmAeBptU7NORwksmJwIKaxGM8hIkFPvEv1qtJxa16x0woP+CBnOtD99H6MYwwraRz+18pa06612x/kdxZe+VmO55D4yWk5u+6CLQs49xNtmfI44aSzVD3LGfTiy6IfxvX535+bmJkx08ucChomNAdB836Z0pRiEypNXa4+9ZFa6ap7fnhFQn8EIu9KRATmbJM+SJUZUZ1DPztQFZyNQql3nKqDy8XBO3lQ9t7xTkvO2Kdda+8Ly0P1X5bFaglfMHpIqkC5LQaWMbS8LXgoMBSz3cJlisKgceaXB2LFd6nd1zjXlhJMCy1hMOlRWhw7xlDP741uSJX9uPosmtI5WztgKAhjrbPJeSryhxht8Yd9fmqfLz0egVLvO1QBP4b68hXN1Wouct025ltgXlgfrv5rDsR2Yx+O3fHeXJXBSAVhbL3SsbvR1gXzkjiWIUaNIaUbl6zxxR5tBtNM5hgmucJ6UbZJhyYhJIlUujRlv9hUrj1RGZYAp9VE+s5+npadaKwLqPwxE+lBy2W+i3hj9fGtodn+eWJ4n70GgcLv2lDLMlg5uGCcg8rZJgDKTtU8sD9x/2WNr351roXHRuppxoYqwX4S9DGaQ7HSOtcd1uEcjmfsEecphczNGkJHdBJ4bg6PSYIixQbUJYuPBYp9Joydpc0lyTGQYZOSNlwn3NLxZJFn2sKT22szKz4VOAoGrWkvrt7OVVt9hHxH9ORwPs/NzwUUIFGvXRVq4cAoBb5sUKvN454IlD5HJB8nFBpEy5u0ebs4x8RTMpM9T8RDlylMGhhahIuXdnBtPx76nb5bWGu9SkH70VHIYPrYMV+ofkMGF4LQdBN6rqpVBXajKbCJvPQwUytezmYZA6XadVrqnHkLA22YInWlxJ4+l5nKcEGynSa0gFVkyYxmqtTGpxtgKHFtzzJJXBXiqbu39Ec+ejhtPlHh/1OU3B/HwDr1oGNNPWNLDeOHYkPKdZWA1GfjJphBQf8GAX+wdMtDIT4Glaoz95BOPpfXj/hCgHZR7sXbdn6bby9nbplybnzqW0h/HBttr2HuZpMskN5OpjHM8HLxOnaQl8pLFkMLQYVksZZBVZSoOEPDqDE4YdX54q7i5IcPyWGXU6dgyxBTn5AisAgH1TdbC+YsUvLE+Ka+iVVwJR8ARWCEC2ADM8b20yCBSrmbsYET00ZDRNFleN312h2MMPVJgGQ033hBhNBF6SXkSzz6oyuOjI/njDXJDqBc1j1gLAuqvvQ8Ea9HR9XAEHAFH4JgI2Pw+pMPFUGShuPsL82nJq1L8JcVLBSy9Vwq8zWZLZ62ixMfTM/jkXMvyTaFw+QsjCwMJw8jJEXAEHAFHwBFwBM4cgUtN+nhweH13yJMTw2Bf1f0YRwTX5v35N+DFp4vkpbu9OsfXaFN/1ImRY3uZ4rLt+jedxMtpGFLQFExuJaJf6cX+ppRhVeGj+NAQM2mMvkHXnuLP64uaVnM/OgKOgCPgCDgCR0AAg+ha5c76yCCyCqidMhyM12uQTJFXWpbKdjqyTBYS3hwMDkK8oRrP0VD5yKBn/IEm+FBc1i13wq/0TRk81APdZn9ZVvJ3JqjhSR0BR8ARcAQcAUdgAIESS2Zs5DSPSliUeYjGNnrmyrOnh+UxM7TCsjrndTr0uu5EfmLghUn9K/wT8flvl15j6lMWfuYIOAKOgCPgCDgCp47AZYEK8KdqqU3LeJ1SxkZcZK48hg0foosNnKs6w9jwMv7QshyeoJac8scTRZjlNat18YMj4Ag4AieNgO6FPFByj+QeyjkvsvCZh8Wec+Xj5AisDoGLpRppcPCGy0cdWQKqSOd4cfg8Nh+Oqwiewo1C/C2fLHllwj4f/qW6IeVFmZTFq/exoQQfivkVU+mJZ5CbJ2tX8zDQyM8HfYWU/zgCjsBGEeBB94HuhezVfKlzPnrLflMnR+AsESjhIQIYvCl8HI4nCDZRc3wcGhU6Z78RS1Cp1+Rz5H+VPG+MNR9hVF4YNHz5t+XlEQ/CoMEYYhCniCcfiO8N/XB7Wv0prW0Yr1l+cAQcAUdgkwi8iGp9X9d4i5wcgbNE4M7NzTZfVqoNqysdj7I0pnLxbs3eVH2WvdEr5Qg4AqtFQPesv6UcD77+3avVtpIrtgSBUh6iJTocSxYPUfxW2iF1wXuVXM47pBJeliPgCDgCQwjIAMKDzosm/AeUG0NDYHncSSOweA/RKdZeg9r2D/Utp+29WtKBDeJ+c9k70l7AKSOgMcKyuFNhBKbgqrTsH3qswPYE215QWCPPzhE4PgKbNIgEu72Bltp7dPxWcQ0cAUdgV0++vPG5eRIWHRzE40WVuQYjy/XZxo3S4s3mxRaWzDq6bL6BHICzQGCrS2Zs7OZpx5eszqIbeyXODQGNTfbY3dcR70SS6omZt0K/POZYVtn22RFeKHmggNGQ/Q2zTPm3Sodn295+5RyatQdSeeGh5k+vnyl0PNXikf//FMKXY+x+yQOl6aFTJ0fgPBDY7Kbq82g+r4UjcH4I1JMxBkBnsq/j+Lsd3nZiYsZbkfrbHrH3T9KHz4jwbZ5qP2KtHzweuEaNolx5pWND8z0FDBXypTzKNSNFl9OpLh+jp5NPXWZjbOr6mUrA+Ps8lX566S7hCKwLga16iNbVCq6NI+AIhAgw6Yaf12ji6omYL8zvdM6Sz9GWb1Q+BgLLVs3LGTrn8yJco//XCr00UX70/w17CxqOQE/wTv3FEDg/l554viA+p9IYSBXHfxyBM0LADaIzakyviiNwJgh8p0k4NUGvrXoYDKmlI17W+EF1wFjqeF6CSiyVD7Kadyr9+L7bfwp8j62lq66pW6p+8wpzKUdg5QhsdVP1ypvF1XMEtomAJmH2Do0uNa0EHT7dwdJdTKY/8UO0VH4o7ylx6PvdFAFP6wicIwLuITrHVvU6OQKniwDLTKt/+1OG290MiNnzk6Q58pKpluiUIV+M5u0y9hD1enAUh8GFFwrPD+nxApmxpsuGwBvcO5urmxR+4ghsAAH3EG2gkb2KjsAJIcBGaTYQr53M2GktM0VKDxlNU+XJ67UMmpcKvP5O+Ks2eqJiq/1V7Avib4j4X0bSsleIzd4pAm8MJidHYNMIuEG06eb3yjsCq0OAiT+1DLU6RTMUwpOzhBp5GTWtz4ToGk8Pnp3O5nPFsez4jY7hPqz34rGnKbWMB95uEAkEp20j4Etm225/r70jsDYE8JwMeV1m6ytjAGPrrQLHXOr7s+cho828P/Z2VqqspfLkiVHEH17zkcVwKYzPEuAVCskMnlTdkU3xQ3k/dwTOHgE3iM6+ib2CjoAjAAIyGq516HzbaA465KWAaMqQMF5opLSKmSKvtG8kfE/HPt2tvJ3S4AHi+nWrQBlO9fWHiM/l3ozQRFnOcgRWi4Avma22aVwxR2CTCOA5aSb4lSPAkpV5XkJVzUM0tjk8V559VSlMqnJkBIVGTvUqv3jXoUI6f6Lwj/gpI428hzxWUVZ+6QicJwJuEJ1nu3qtHIFTRYAJO2VkrLE+/G0IxkpMeHL4kGJslMTpcuX5VhB/CRITXp/Y6ILHfqGGJPtQFwSMpRRhWKUMpVRa5zkCZ4uAG0Rn27ReMUfgJBHA2/EoU3PbdFx5SjJliiWTocFr6h91ZBNzRTrH28I3fZ7ecqqlOjYz3yi03vLSdZa88vlFaVubp3Vtf8zaGDniUTbGZINHzcPw4m2z0JMkVkOVAddc+YkjsFEELjdab6+2I+AIrBOBV1KLCbyXNLFbvO2L+V08PBxvdDz0t3QwJvgzV4w4NlFzDP8Qdac49huhX8tzo2soR56lLsowowiDhyWu+D/FDA++N2QGE56lvo3hiqoIufCNtJrtB0dgWwj4n7tuq729to7A6hHQZM53ccYm8dXX49AK1gbTlY4YWVmktHiUMCRTS3JZeXgiR+BcEPAls3NpSa+HI3A+CPT92ej51HA/NcHTE+8pGiuJ1/PB28kR2DwCbhBtvgs4AI7AuhCQt4JlL76tg/fCKQMBYWX7h95lJK+S1PiC86GXGXNV9HSOwEERcIPooHB7YY6AI5CJAJuFbc9Mpsimk9nbblM8RODre4c23W288iECbhCFaPi5I+AIrAIBeS14ZZ03o2xz8Cr0WrESbNhu/b3HkK41ruDrr9sPAeVxm0Lg/wpVaQzKjlJqAAAAAElFTkSuQmCC\n",
133
      "text/latex": [
134
       "$\\displaystyle - 0.037 \\rho^{2} + \\frac{0.042578305 \\rho \\left(- 0.000125 \\rho^{3} + 0.0025 \\rho^{2} + 0.05 \\rho + 1\\right)}{\\left(1 - 0.05 \\rho\\right)^{3}}$"
135
136
137
138
139
      ],
      "text/plain": [
       "                           ⎛            3           2             ⎞\n",
       "         2   0.042578305⋅ρ⋅⎝- 0.000125⋅ρ  + 0.0025⋅ρ  + 0.05⋅ρ + 1⎠\n",
       "- 0.037⋅ρ  + ──────────────────────────────────────────────────────\n",
140
141
       "                                             3                     \n",
       "                                 (1 - 0.05⋅ρ)                      "
142
143
      ]
     },
Markus Holzer's avatar
Markus Holzer committed
144
     "execution_count": 5,
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ρ, φ, c_l1, c_l2 = sp.symbols(\"rho, phi, c_l1, c_l2\")\n",
    "\n",
    "critical_temperature = carnahan_starling_critical_temperature(a, b, gas_constant)\n",
    "temperature = reduced_temperature * critical_temperature\n",
    "\n",
    "eos = carnahan_starling_eos(ρ, gas_constant, temperature, a, b)\n",
    "eos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we use a function that determines the gas and liquid density, using the Maxwell construction rule. This function uses an iterative procedure, that terminates when the equal-area condition is fulfilled with a certain tolerance. "
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
168
   "execution_count": 6,
169
170
171
172
   "metadata": {},
   "outputs": [
    {
     "data": {
173
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAAVCAYAAAB2SvVRAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQFUlEQVR4Ae2d65UdtRKFD7MIwJgIgAywnYGdAY8IDBnA4hf884IMgAh4ZABEADgDuBFgJgPf/WlUfdVqdWur58y1PW6t1UevUj12VbdqdHpm3nj+/PmpLF999dW7uv4ux472gcCBwIHAgcCBwIHAgcC5EGjlGhclcxF8pv775dir2pYtCzs0dgcQ9tjEOl2fwGPP+tdljfA5cH9dnP0S2HkT9+W5Y/glgOn/psJtx+4m4m3EOS9a/oiuBi17KjnHVN6MliY+UPtt1d/EWNQa+zq3/1H9nq6vNWafooysL2hD/A8aexodavW/Lfp31X6ssctijOZvGiN5iLWRSNyr6OBHkvJ5Hr+v+hl9jcdapthokfutxlUtyqXG32JUNfw+zRTIpQ9mv+axqdIY819MA6cT/Z9atNBoPHwRSxb4VHSrPhMvW8/MEwx+0nVPay8ZaxQL9yy7hznsbdxrXSSDtZ+qDl/UJLv74jnkt5Yg8Qhf9nzk4HQSP+e+SKo4smudtWYTzzzfi48Z2y2emnPsGYoP8exiLgWtGMYQ8XOeHclmUzY8HbsTz3N9SCY4dn1n0NnYObpneYt7WOND95/ou34XjePL0XhzeE5QFHrGWP18H5IPE/G0fBsCM/0C88zLskc8oIvnLr6iP9sDRfMrsnTxg/538E8JiTrJuapbm/WfonuiuZ9ZkGn/VP1IVzcpEY21XnQozA1BIpA2btX0uUiCTuqjJ/xIClLipBqw/6OaTbLUh6SCwjzj6I8dl6qnoj5y4fcoBtUmeMPGSCIeaIx2KSOWPFQjbRpai47YEM44qU2y94vqD3UlHGOhapw00TKuPgkJpzkTrdpdfIKnaLuYw1/0XT0z3feiBU+SNfTYKl3csy0O5sixcF9RiNgpk8oVsl3Dlt/WOJs+smIz+8i9L06O7BW9F3hm2SPxUbNe4+naY8fHgN3dGMYI8bP8k2nde9K1u8ZxuO/6zqXLCljYDSi7iI+81r7/pL+DvevLkXhzeZ6kI7TYurr/Zbst+YM+y6ynqol51rH73M6yu3sL0kT7jS722h91XcYJCRtwmZUnzUTwiRr1xshJABsl9NMmnhZUH4PrAeEHrYkEAG5smmUCwIPvrmimUxy1n+r6Q+O1Pox/qPFewfY6IQBMbEendOoBE40t7NUYgcRcJA+sI+MjAYmxsOkLzcUYa6DlZqnLYw38pmuiVdvBJ3g6PrP0lI4kcAlHtZ2v9BzcbcwlG5sc3CGdStZ16p+zId4jfluIzusdH7k42ffFgOyZ3mt4anw0Pia+azxFYNsDM/Hpxodo8JmDOSydGIbO8s+A7CG7UeA6RXpZvnPpsi4udl3VJXd2nB8LMp7WczPTOn63fIkO4tmNt6yrzVP01vPdlS8dLd9mPadK65qYZwLXHmtvmYRe7d2J90Ue/EiKpCOTgogmG1HrJ8zfNf5Qa+5AtFGs9eLDCQInGTMdNM4pTOl86MoEJUSjo6NP0Jc1pxt/NWwhiSCQU8KhNja3Sp2powvBwJWKeEztGMs1Jz+lfdX0VVfrXXxYYGEuuhE9rxQ536eLORJd3CfthBexNPPBNHmehuW3DVGuj1ycRu4LV/ak/k3g2eE5Yo8bH8N2TwCsN1z/uLJH7F7X6hbMdOJj5P5zsXd96cYbXrB4ytaR5/uI/KFI6GBu2yPCob1Fctn3yUHuXGQwWpt8KPCMRlWCHsC3CvPOek4oLqXL2sZ9QtksqMXvnzx3P9cjVfoaZkN2kqv58rQi8dcYWd2TUpjG+F7sLeoYV5uAo9SnUATXB5rnNCXsgw6+JW0XHxblYmGOfrpcPYP3uWoLc4RJRwv3SrGPtW6W3Fbz1+26fluTY/lIi7s4yc6IG/e+cGWXut8Enk2eo/YMxMceu0sMWu2uf/KiruxRu1vK3LKxZnxkG0fuvy72mafly4F4g63FU3T2831QfjbNrrYwh4llj3Tcs7eQU3z0pj74CR1Bs1LcILPxqnO36k/dwfUkEn9rDT/ZfqyLBIMseHq5U3MkLBo6tWS+zYRKnGakjujTMa06zDPHOyRP02T+UJ8MulXQ5VTTB6HG4ccLOfHCYUzNas1zQ6TjKLVnm6T6P3NpnoTlX7Xhhd3l1z3qpnc3NvGBSOtjc6K7Vlr4sXZVzzVGa+PitYm75ndhjjyt3cRd8xw5lsncmpq7xyXD9dtChtbaPhKthZPokNPy6+y+GJEdimvN2fHc4qm54fs8dKXW+kV8aMzGPHhpzWYMZ1ld/7iyRXctu0Pvl6F2sNvSU+s3Y07z1v0nOtvvou36sqWz1i3iLegGeHb3v+BZ11vya9qtvvhsYs7aAXtmorTO2VvIQR5d6AMw/ppxuOrEA+6yMRdDWw4fWR987kv5z3XxogtZIwkJm3UUNm8CoC7v54HgQ5c2L8rAi42ei5dnAGeziAZ+yNlKNkgyuJoFHrpwMnaQBP3RIhQNN0IkKvBDv1nSpH7Y1cNnBHOxTUFm6ZmIvY9duJuYo8Eq7uKBz3iwxwmep/EOKtNvLc7DPiqZrODk3hdDsm8CT5Ona08JTbRb8TFktxjtimEUaPhnRPZ17A77X3S9G7uMn3UPC2fnuTmC/QK3hi8XNBpoxVuLLo2t8HSf7y2+Q/JbDKSThfnK2vc13twrsVVXdw/MfMlB3r3QB2C0jnsz3WYVP4FtEm1M8mvG4QyUj405lvyoxvcFzWMm1J+SCrUB5JJxlWkj0jjvn8T4SW3myMKcn55/Eh1Z+PTyrPpT0TgOeKh6cbIURJrj5S6SIW6cH3SRDJXJVSLNY+jJyQj84M07LYlW9Qg+WtotM5+Jv6Vnl2smEL+9uG9iDnvx7uHOr6rVMeSqPkQnOfhn1W9DzJbEMx9V0y2c7Pui4tXqlrJvAk+H5y57jPho2Rtjk93iszeG4dXyT8hYq0P2LrvXmL6I8Wtih8pOfJwk51z3X2DfgmvTl9Kh9zzq8hSP3c/3nfJbOlmYtxZqbBUj6Teyt5CDpISELJIHa122kpTIPOPdjXot/dH1UzJRMONNahx2nzEZiJ7v6OLXZz/TRXZIUPyui9LicTVz9ck8f4yFNc2iORIWvh5ZO8JjHacePVnQpSJe/OSD7unXea9Gkz0cC/Pw41QImY/UD7llIsaSlrwZPqIZxRy+U1nTcyLY39jEXXIdzJG+irt4gKWTbO63Iq/Msly/1fJ2+2gNJ42794Ut+ybwdHkO2FNjuxYftt01w6K/GcPQSe9WHNuyr2F3oeZL2exil/Gz7mHh5D43bexr1FZ8WZOtxVtNl/odns7zveY7JL9eTD9jueu52bFnJk60zT2wIML+OxfFwKwpBjzkKHeuqtlnjLVATITu+oIu5M0E5c6UQECvi4wufRWjGkMjy036aIx3MNis10roP5vXGgL9rmoSg61Cdt4Mdq3lpIdTm7r8kQem0x31SahmXwtpLfZwWoKOD9UPXKLW8KIkfAraln0xFhiN6LkQ2BqQ/GHctcbFHJFN3MUD+++oXo3Hlr7XGOv6bY23dAw/hj9K0hhb2NHDCb66eveFJVt8zo7nKE/HnhK43F6LD8tueEjucAzndc04xo6sW/g2d1MVY5O/d9pd8nxh7b3YobDWjsScdf+BZQYjcC6xibEJ+5jUuqYvY76om/FWzE/NNZ6FjqHrtKZopOd70Y+mLT8WlLVkj2BeLj2t2QOR5vbsLelg5E2tZ2MN58CvLPEVQjlGO05IVr+yyAvc9Wt0IXcRNDGRaxIA3uwNp95Xv5UwJL1F97RaD4g49z3VcULBWAoE1ZN8tcGK8QWPzDMlQqLjN1hCnzz1vyrzYRNd0GiM0xISkxLntaCE6aSf2mtYlrxYY+kJ4UAZwl02WpgjX7RbuIPNA9FwfFgW4oITMcbBdJb8lYRuO+vh+m2NreujtF4ybZwqgdhf3hdMO7Lx47nxPIePWvZg0yn7Zeu+dOyG1VAMZ9k9/7iyYdcqq3a3iF/g2DB2ha5WfIj+ia6R+28Ie/deM+JtMs3guaZj8Cif72lsRH4wadQW5pI1e24a9uzZW3i+PyMhwVgUaxUe5GSjdbmnAb4fWmymFaG7niOjejOBFXL4yQ+HnTIQ/PGgd0K2agzh1AHaKN9pfAZinoAu8QpCatFyw/MArtfwoKnfSeCmo7QSHsbBpN4EGI91YQt2cbFhLgJO9NgVulr4IETFxdzS84ql/WnjLptHMEeBwG+Bu3iBU2A1KavxfxlXPSWZTKrffKBNCzcaWjvitzVOro/QtYuTaIhT575AH0f2EJ5rRpbj0tHmOWhPiFmNj0zg2A2pHcMQS9euf0RmyR61W/S74xjdb6AMYVfKly0j8eE+NxFhYQ+hdHB8CSmlF2+JyOQ58nxPfF35QdyqRzCP9aY9e/aWu5Lx94U++En/QQgsawlnM36mmgdeKmqzUX6kK72AxSBjup7rmn1Nor67ntMANo4p+VF7IUc0JE71hkTAcVSNHVH487az78XU521fSr05wRMe2JDWRa0x+AJuWdCLUo9fjV59BfNLdKjFA/xYV/NDl9l7JepDz5EhR8cpUVHt4sNaC3PxJ/ly9RRpKvHVGMHTKhbu0nEUc2T1cG/pw5pYl+Ylmz6/Yj2L1dbijTHXb9e9L1yc3PviJLvd+GiZv8CzIurFR0Weui2etj0Fw/DzZTE2NQfstmIYxuJp+WdAtm23eGLvdeN4wkcN13dbdBZ26K5rsV+UyhTtVnxY9x88JMeKd9FZvqz0otuMNyZcnqKzn+/wzQVcKKvyr6bT55bPCrKp2cLctkdc9uwt97Tu6RtffvklWSGbIu8tLIrGUY5EAcP/0fVA1xONlwkAyvJrOyQVvGgzFfWt9SwQLXKgp7DpteRE0hJ06N766ZgAixMPeJHI8PIodkxF/fTrRtPAvMEpEEBNRX34spk9VptAWhSNcxLDTROFNfxF15ae4P+FrjLRWrOpiw8CJcfC3NVTdCRsFOyCN74nWSJp4oafivpd3EUzhDnMM99N3EMJ0ZKMogf6UvDT7xqP/3+EfMrWPwm8olj5FC/Lb6LbfV/ktdjRKrPYFK11X8BItFZ8hFDR9/C042OAp21PtgmcNuPDtVt03RjOMu04HpBt2y2e54hjy3eS5dKNYLfYLwbiw7r/sp+68Z6xtO61zNOJNzs+Mk/r+T4g3/IZ/CjCoHef2/aIl70HZtnw/vSN58+fowgdfnNllmRAeJQDgduIgGKdUyseiJe30b7DptcDgSOOXw8/32YrFcMkd/xw+95FNpTMbHaycZsBOGw7EBACvDN0JCNHKLzqCBxx/Kp78NCfbzLSyWBKSPRg5th98+9zHJgdCNwWBBTvHOHy9eNRDgReWQSOOH5lXXconhFQDHM6Qu6RvvqPExKmeeeB75COciBw2xH4RDdA86/w3nbDD/tuFQJHHN8qd76WxpBzTN/OpHdIAoacrfDfZ4+HdYBy1AcCBwIHAgcCBwIHAmdFQHkGv/nKv2iZ/uzFfwEZmXhdZ0r+4QAAAABJRU5ErkJggg==\n",
174
      "text/latex": [
175
       "$\\displaystyle \\left( 0.0695273860315274, \\  8.02904149705209, \\  115.480272671423\\right)$"
176
177
178
179
180
      ],
      "text/plain": [
       "(0.0695273860315274, 8.02904149705209, 115.480272671423)"
      ]
     },
Markus Holzer's avatar
Markus Holzer committed
181
     "execution_count": 6,
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ρ_g, ρ_l = maxwell_construction(eos, tolerance=1e-3)\n",
    "(ρ_g, ρ_l, ρ_l / ρ_g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this information we call a function that assembles the full free energy density:"
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
200
   "execution_count": 7,
201
202
203
204
   "metadata": {},
   "outputs": [
    {
     "data": {
205
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABO4AAABMCAYAAADa4/hFAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2d69HcNtKFRyoFIMsReJ2BJUVgOQNfIpA2A2/p3/5z7WZgOwJfMrAdgbzKYPeLwPKbgb7z8AVokEMOG7zM9XQVBiDYaDQO7k2Q8+D9+/c7kxEwAkbACJw3Av/85z9/kYb/kP/2vDW1dueAgNrJ3+T+dw66WAcjYARuDwGPQbdX5y6xETACRuDaEdDc9lhl/EnuC4Xvjlneh8fMzHkZASNgBIxAPQKaGJggfpJvo109fDeXQu3kaxX6k7UKLnkYAV/JsVgxXQgCqq+9NkAdyv3tQopgNS8bAcYNxiKTETACRsAIGIGrQEDzGsa6f8j9hzXVMQv1wCfujgm38zICRsAI1CGgSQGj3f/kM0mYLgQB1RdGE+ruqcJVT+TE/69UzD/kfyz3L8WFTs+J73PxP5c/2F4UX61XkklZxuhOPB/kmwpP6i8ejEd/T2lY+HBNOX9NcR1P8ZN6R2UmvozPM2X0Tm7wNGsNb0fhM7iQ7n9KDbDNBv+8wBxsk+KfxHisWEo7WeekFd+3hYwnCr9UXKd/6Lq2bUzKzHlKdtYzR/2guIxPjtvN4DvYVwt5B/lQIMq7Nl9N3ol3sr1IRwx3jA/fkcZkBIyAETACRuAaENC8xnqbdSvr9KOQDXdHgdmZGAEjYATqEdBk8Eqp/i7/aX1qpzg2AqonDCPfy2EIwiDExvYDxXcME4obJfH+Rze/kf8zTPKRSdxnCh803iXe3+R32kuKn62X0mPsoCxD+b9QPEavrO+k/kkfFjvZcLdTmAUQxkFePciywngqDbwRmRiFvhX/Z/IbUpjyYWAA49ZwqHCY917Sef1K//9KI4xjYEPdgStt605+QwqHMc5p+r5kROscPrD/NzLk06Z+k8OQ2LQt+aF6TOnhnZSZeKlL2hdttalj+Vx/Ir9ddCsc4ksyJ8tdw1fDKz1PkrfyrW4vSddP5d9RRpMRMAJGwAgYgWtAQPMa60dOl39xjPI8PEYmzsMIGAEjYATqENAkwAaSCeEok0GdduYeQoCNKZO3HAapH4Z4DsUpHYZaXmVsDFfwKsxml+vyVBG3hoj2sseHDLnZepGR0mPUwojcOkWT36+Ky4a2qP7w8eotxrpM2WD2Okfofo3eIZmSjc6twZC8lA+n78C5f6qwhhdR50ZvVTYMxw/kPqacch3jCddys9uG0kbrHMPxE/E3RjuAUpiTbr/LlW02Wo+IiMqEl7rldF1uZ8Q9lusbo0N8khMqd5QPZaK8a/NV5j2nvVC/9CWTETACRsAIGIGrQUDzMetHHgCyJticbLjbHGJnYASMgBGYhcAvSsXpmP7GcpYwJ7oIBDDS7r22p7g3ci/UFjA0HKIvxbPFK2nkP0SdE25iiOpPGTEgtUYk6d2GhzIKxEVlckLwvwNYYtB5rHgM5plqeHOaW/OjdY6Rdmgso97Kth2tR3AOyVSdwsfpvk7fUDzG6PLkZYiPjEXRckf5tpC5Rd5N4Wt+hDG4MzZNjV81Ys1rBIyAETACRuAcEMB4x3p48znu0TmU1joYASNgBIzAXwho8OfJTed0yl93HbpiBDAUdYwLqazZ4MH99jReiYPazJgRo2SbFZbsvTwVxwmab3oCQ/orLUay9pt4yFAc+kPl6av7mMBvhUzy5unomKGwXHjV8Aa0vEqWyTpPdUPh3w0g8EeKeyaf05uhtiG+XE+TMiWXE5acFBur86RCmA/+yXKLh34T5dtC5hZ5o+ccYgz7Um5ofJsjz2mMgBEwAkbACJwcAa0tfpbjbZG9tzTWVu7R2gItzwgYASNgBBYjwODvDc5iGC9HQGGIOKQ03ysbI04OYfTYnKQrp9IwfuU/eNgt0V9pMTA0Cx6FV2n3YzIVzymkIeJE1k73OfHVUA1vToOvdJSHfDAUgRWvqGbjqy6PR8q3eaVTOX4ohy6c4m3LuEQTycnGs0NieACB0QyeofaLXhC67ZHSDbaNSpkYBfmDH+r4KzmMhXzXjn/qLvtMiE9pouUO8UmP3doyo/Jq8oZ3AYEzY9Qq/XuBHk5qBIyAETACRmBtBFgP/6K5l5N3m633Hq2tteUZASNgBIzAfAQ04PORfDZ8/dNM84U65SUgkI0adweUPWQIwOgw67TagfzGbmFkw5VUrb/aOoYUDDPP5TAm8b2zRTRHZkrTGNimMp/i1X1w4ZXb5jt68ikffyTQOWE4lc9K92kvP0qHpk3Jp4y8Jtz5E44FedXUeT591s+ONgB12rZ0jLSNqMws+5nklsbmP3X9Ug45UJQvWu4oH3lHedfmq8kb3rnEH6XQF0xGwAgYASNgBK4KAa0jeGOAdSxrwLEHxIvL/HCxBAswAkbACBiBNRFgw8+x62azvaZgy7p4BPLppKGCYHR4N3RjzTi1S4w/L1ikzJDb0Z9Fjty/5Vjk/CD3H4U/nyG3TTJTJn9IQJ9r/zihFbgfGOVNun8uvzHapaQYIzHkHd1ooTwx0LXjiMI8BabejmXgBYJc5y+5kA4tDgp/oqisX+cJte5F2sakTMmhX0CcEO2f9vpR8d/DE+VrJMV+crmnuKN8yInyrs1Xk/dYeRmbGDtMRsAIGAEjYASuEQEOXLAGzOuO1cv4aHWJFniRCKRGxvvZUF5c8SQ6L6rv7/jXCCxEwG1tHEBhw0aW/rfZ05rx3H3nxAgcMrrlUzb5e2BDqsJzjPEao1THyJKUWaS/2n42VvP6Iv+EurgsEZniwYjFa5STfS7Ay7+ctie6Ei55Ph1dyEku936TG+VJskqPf4J9W0YEw9Qdhte/yQ3VY1BMwxauc+XF67IfKRWvkTDOYVziT09wGGtHdRH/YNuolDkkn5OQvErMaVUMrNAUH4bPaLmjfOQb5V2bryZveOcSuNa077n5OJ0RMAJGwAgYgVMgkB9os66IPAiu1tGGu2rIrjZB598BtSBmM8Oilu/AmIzAmgi4rY2jifGcDe7bcRbfuUYEkhGCog1tbnPckFHh2HAMGllq9Bcvhpud/H47x3jCiSxcfn1RwWmaI1Np8p/A8O2tgzTFq/voTD1xiqsk4qF+We9j9au0d/KethErBCSTf6Xm+3JjcnObmp0besuRfkhWjmvbLPziLU8j7hSXX7lu+HRd1TamZHJfDh3Je4wwYvKaC/cP8sEgviwzl5HoTDkOY3CIj4RR3rX5avLOBZzpH+vBwkz1nMwIGAEjYASMwHwE0vyM8Y51ziaGu4fz1dsmpQr9WI6NPY4n77i8ENomU0sFgVfCOW8wuGYxzWK2WUQTcU1Em5JzOztNpd5UW6uEmD6Yn9hUJjX7FSBA3ecTWmVx2PRCh9oGJ3E2nSsZN5UH+t3JDVFUfx4K8VrsmvpWyVTeGCA/lt+etFOYOW8P/yAvcni9s48Nf4aAEac1YOn6GMQpsiF8m7YkfUYNiZXKRet8TCxrDIxmGbeqehwR2peJjkNY5OS5bqJ8pIuWO8q3hcwt8s6Y1frgzxhlMgJGwAgYASNwrQjwOZXBteQaBX60hpCVZfg0zsqABsVhHc6vigSTXDSb29npqu/W2loIaW1c2WyyueGkjOk2EWDCzyeQSgQ4NTVkFCp5MD7sGZ1KhhXCGIOgsQ14VH+MNKWxphGqnywfg0MthWWmvvZcfv+1Vox5ne+gVfDuGd1TWvr12Km32jLW8H83UD7S7+lZI3SAN1Tn0gVseZX4I4Wpq518xjv0KfGpqceoTN4gQM8+kS+n4nJ7i/IhJ1TuCr4tZEZ1rMkb3jn0RImygXROeqcxAkbACBgBI3DuCPDWBWsJ1jad9eQaij9cQ8jKMjY9jaMF2tYbm5XhWEfcVLl1n0V+s5hOOWJc4ZTAWk/l1ynIelLcztbDsiPJba0DR80FgzyUN5H3V/69VAQ+TIqzYe2Q+ggnft/LcbqoJV0zyb+Tj0GiIYUxbnwp13yI/z528Jex+vngnW5ktV5FcnSByrniPka/FfpjMOsYqFOZkf93hYfkj+qdFAjJlGzWABg0qINvS9fPO8orPvRGblvXKY58KM8p5tGmbMq/JenBP1ZD7SnD+8v2dxRjyiO3pM2CT9/gO4RPqB6TxiGZ0pvXrjEUt0ZxyqO4Tr+K8pG3eEN9Ncq3hcwt8kbPgkbbS8GTg08VOEU/yPnbNwKrI6A+xhh0tnTu+p0tcFbMCMxEQH2O9Stu8hMsc7J4NCfRUBopyj92dSZlXbMw4hsrNU/ZNjuNIz1YtKJLjT5Dxb3EOI5t8k8nk+9ci4cTAmwcWWidFaGbnNvZWdXKnjJX0db2SrV9BIM8pz9ubnxSmVn8Yrj8UWEmvIsl6Y8xAsqGWD73QJ3+Ir95+kYZU9zQKWfGXU4EY4T7Qw7/U113xj3F9ekHReS8+/d2Sp/vzdULmehA/bzhYoQm9ZcuPCh6IcdTyUy0Af4FtWO41vWk3ggQX1QmBkPy4vt2fepjHOXNmP5DemTjGN+HnfsHEn29qq+lBw/eaEcZY4yK7+TaE29ZqHgmMRbPojar9Px78IfK87V81oYQ+vXrO1qPO6UNySQj8dK2+njs9asoHzJFk239ni3Mt4XMqI7hvIXRZHtJ5S49+gjre9NMBIT75PpXPBczl0bKMxOqoyST/pewpwyvx48C2saZXHqb2hieixF/ynpcKW/W9thSVqcH79+/b4VK2fw0ks0Ci04WOaFNpPj+FD+LsbzwzQuzp7rHQn8WJZ0wOC36kwSlxxA19GpMo5fuAzCLkUX6NsLO9EdlZJJh8T16dFP3mPRZ6LPhmF1vhyCQXLezK25n1L3q+BraGn0hbzQYz7je22gqbhUSZoyhv8vf5CnNkJLKa3ZfRF4kvXgmcRQP43PeDA6pyrj1wdANx/2FgDD6r65OZiz6S5PbCgl35sxn8jGSmIyAESgQUL9gDuDBxaJ1fCFyUVB6bD7v9RVUnuwxOH2b1xQNi67BJr+y/0xhjOs8AMh7qYaPH8VN7rPEE55Lxcu65nWbwf0ejgdNHWN6yntyHs9ylD6Er/gmy5NlJh3AkNft0YUTtGMnh2HflJT3ZnvKKH7RAkre5Ho8KmsJ39JyRdKLJ9SmxBduz5Q5knfiyw/KuHwi91JpO3vplHe0z4fHh5R/2JZRlImk0A+K2xt37m81GAyOYfl+9iVjki/xjNpcdD9Uj+Qp3lBdJr5J3GvyzmXu+5LBGPi1/Af9e0uvH2UBEs4rO9/I55WCnXwGdD4ezVPKiPGOyQaiwuBHDvI6DVZxYVJaZDE4LloMS04zOcnvyEnxTALozoRJ5V8tqbw8naZOB0+1KL4ZIOQ3hoN0vZMfqf8QbpLldnbl7YyGoHq+6LYm/RkzWDy3i2yFGYvYfGAYacZJyroGSR754Vbra1N6Kc9FfTGSPpUrgiOnytgsDJWfkxp5sp0q1q3fZ7FAm23b7a0DcqTy00ZXHROOpLezMQLHQIDxOxtzjpHfaB6ReWs0sW4sSM8mtbMplizW3LzS3j6sUxic8t6rb0CL7LNq5lIeRHbmCl03fwgovx3PFA6vh8Rbs66IlEdw3JNkgx+HK8jjZJTw4OTw6nvKVLYle/E9XCTz4Hp8L8EGEUvLVZF+sk1JVrg9A0Uk7ySTdkl/bt5sk48N4//k02abta38cJ+P8oqP8oRtGUku4xFr82aMkc817tDDlb0xTPxDNMhXqedkPZJxkjm5xxBfGHeJDeU9VPAijofo6PdCrj+OF2z1wcZwJ6G8MvJYfjlQc8KBa6zH7aRyIAs+nL3a0w/JAmQmsDVOwCGntII3xVAedwo0OivMEwk62bUTOIBHf7JuGrXi6QAZB3gmN8ziZ9PCcezRk3y6v9N9t7PbaWdU+SW3Ndoq30HEUJfHxTz4vta9HEc51yD6H9QM9vfB7X6X9sWK9GEcJXNvnlFcg4v8tfHeDtwTShZOvGZIm2U8HjKCnlC768xaOLNopp2+uc4SulRGYD4CjEVKzXjUWXPOl3i/GUoyD645+3lIh0Vr0LnplY79xRANrcVZg6Mnm9/+KfPQPkvpJ+fSlMeQAeyl8v1NrpxzQ/N4klmzlwyVR7r0KW+s+/HHut5kTzkDv5ryDq7HawTM5V1arsr0kTYVas+UtyJvDGd8Gqz9HJXC6PK74sE+98maPh/iVR61tgzGFk7X5T2NLnesY0bXjOIdG8NI29Ihvko9I/VIvtG6DGGZChLNuy33QIB6h7CnlDg3kYd+hNNBm8rDlBjj1dsBQSxEsRZSoUcj5cdEz+TFaT8MiEz6xM2lL5W+aoKfm9G5p0s4gEe/TpnAaSz42WG4YECYImT15Q2lcTsbQuVK4y68rTEe0vbb9h/sC3NrM49voxPnXMEj6Zb2xWj6KI5jRo+9EwEj5XH0XwhQNywUTcdBgNP6UNXi7D6Jf43A1SPAWLSa0S6hFV1z9sGNzlv9dPm6Or3WDWzcOmuJLEw+a+7/iofylMRYggEsrwvKe1Ph6FzKyZpsTJiSGZ3Hq/GZyvhM72+1p9wMP7Ul9sBDe79jQLy0XEvT98sYbc+ki+bNGzlD63fyKu0oNX2+hhddJ0ntAD0Zkzo2EcVjbxkcDxR/aAxr84zytQnWCUTrcnUsJ9TPbeHQCcYxEQfnt0cpFQXqVGKKzxlzv3z6km6v76nimaiY6KtPfg1pI3ljnWmI/VbiqNcv5do6F079J3tbYOF2tgWq5y3zItua+gML506fSGMJaG9hFHmSqjGPuelyM29pXwylj+Iovr35RXE8IftmMwSuVLBw42EX31L6Wq59+nulxT2HYvFktXnIeA7KWAcjcC4IMAZJF8aiY81rU0UPzVsHhMxJ/5XKz35myHjJOoM/nMCwN0R9g94QTydOsqJzKQY+5ohf5JfftGbe7axxxBNdD83Bp6P/KS9UzsmT6uLZck+5NX576/Ej4b20XEvTd4pZ0Z5JN5l3kgfvO3569Ee6fiafflTT52t4e9mOXjIOsUYcG3OGEh4aw0r+KF+ZZlG4oi63wHJUd/CV437e243y1t54JMGRiSGUsWRxZBF5H8phgOM9fayhe6R4JofcoGlAfJQdXk57IQO/Jd0bmvTa+wcCWJCpsIuhA9isVQbwAJfWcLeW4DE5KpPb2Rg4J4x3W4uBL5yYvJuj1gqP9pvExxM6xjTGQBbskU1L7h+k25SkT87rUD6jY/6S9AmfCI5gx4Ym+qr+HMwPlf+i7wk32pyNdkeoRWFNn72oNcYRYHEWRmCnvnE2Y5B0Ofq8pzwxXHaMYGWz0H3mrSH6hEjd39s/KS68z0oyBudSyfkZJx4MUX8qzFzL6ZDy8yC63Cfx7q2HFDcLX6WrKs++Nn/FJB1YX+RPjlCevT/bSPqzp2SeBB+Ml+jPP3FT/kOnZDbZUybdpcJBGlyXpfJE1kAXt/ebg4vSVLWphN/eujSat/gOGWmwh0C0s514w32+hhfZQcKASDtnjPlKDjvMYD9RPPoeHMPggaJ899yxX8msqsekx97YlOLDuCf+6rxHShUZF0eSDkc/UnQeCO6GWZrYSMbwtH96IMBppBwB50l0u6hN8QySbGibJ0PymdgAFd7OKRfFLSUa6ejEuVT4mumnsFkxLyY1Gvcxye3smGhP5OW21izSJlBqJiMmN/rKczkW0ZyuGSRhysT/WH7zkEE+6XgAERnT8uT+blD4upFL+2J1emERxjEVFSxxB2kh5gdl+6YRMAJGwAhcDQLV81av5FXpNTexB2JDH3lw12aV5krSDj20Cu2zWmH3gdG5VHlx0o79EZtU+NB1dN5Nuo2th6rwUT7QnPLcp+z9Jt1+UjT7yBZzhTHE8RCwMSITFg970A8UvkOMfP7BstmT5jjiR2irPeUc/NCd+oquOy9x71eLS7hNCTvawlh7pvpr8saegaw+kQeEXoOU9Bjr8500NbydhH9dZD2eSVY7xiiM8f6lXGOXgV3h0BgW5ftLhVAoXI9J16m63MtUepNmCPeqvPcE/xVxp2BuQ3/FLgw9DKbPm8pRdgHAYImSDSnMwInBrm80Y8DMT3oaXv3w2ibxWxAVcIzN8Bq6H8RGmDI48/TofT+zdO+V/M5JxT5fugYPGuu5kdvZ8Wqkuq2lNsY3x3C0w+bfxyZUvti2RrlURj5Syj9y8WDhBzn+7Y0n1B1KcZ/LL08GY+Sjzw5N5p308PUjTnw92Rcn9OukFwYhHJEpXsamF/LbBz5Deek+9bAE8yGxjjMCRsAIGIHbRKAzb82AoEzPK8Kjp/MPyMb4xB5p77Si4qL7rEa8+A/OpbrPHMq+jRM3zLfwc+Bib42j+J3iw/M4/ANU4oO8qvIMyCujMm6t0S7dxDjBmjUbUF7rmhNH7X5V15S9MWIo/uC6Q3yn3FP28atdA130elzYj1GLS02bEu/S9ow+Oe+XXEhmu95XmDaX21m/XcKeKbfdvT6fGQq/hrdI1uhG24UwZPfHph8V/73iMw980TEsyofMEEmPqrFB/HPqchDL2rwPFIj+VuJ5gDV+65FYETxG2VLIUco5RENlA9Z8N0A+Rjwmhm9KYYqPnEgpk9SEKUPuODXpDvJKZyqDf16qqRSebr0dEhzBRjw8vcN4AIYtKY7BIQ8WEX2olwhfmwcB5UP95XzKe0070f3SaJHv05kweridZUQqfGG3ajsja8mc7IfiGWprnT8KSHIwFB96reBq2prKy2L6TuXFYNk+rQVT0fdy7dOrJuavfhrpa7kPIX+SlP+SdrG0Ly5KP4EjZWccObTIyfjMwnwhdjnvPV9y9x6o7DE5wggYASNgBBoENGY+OASF7i9Zc/ZFL5q3JCycXnpzgg3dqyiVF6MSa+YodfZZvUSjc2nS8an8vG5no4whiM0sG/hf5UbXI7rXWQ8pTRifno79y0Pl6fM219KFPRD7ojd9Bt1jD0L0V3JvCSwk1mqjuCyQPQe/2jVQaD0uvJasL/sQzClXKWNpemRNtimVudOedU0dh/OGX+4jpclGYgx6tEcc/Qod9khpGCdCfb6Gdy+jbsSQLuzlGLeeydH3Q2NYlK+b/eyryXpEsnQaqstOpuIJ454ShvLuZHJ/0eztBuLRc9b89kgJaWzIpKP2KccNVXLLq/Sc3nki/2kb2Q1kOZysm5oM4GUw4iRFZ1LXNfeQgXV3LC/d3p6U/51yWVOHSWxSqfa+ryBd3uoekxODQ4RmTTySnyf4Th4pX4yzo08LdO8i2pn0pI29TgXMBtKX6N8p9JEuUr5rtjM0n9vWXkkfjFb5iSRH9Inj6Q1tcIguta19QmEGysUpOozXuOZIuXgI0254YlUS8dAYNvd3Z/wqz9njD2nlyBWd+5TjRsf8mvTiDeNYKDK6yMk8kjsbc/SXnLX71E5yO/NV1tW+ETACRsAI1COgMXX2mrOfG+O+HNF5jitZctzieU95sG58LH9UVplxDoufjTL7KNb4e6T46D6rTHtoLmX9hqGhJeXBhpcHsfmVyrzGmZzHU1pkZSxbuUVci8nM8pQyy3BeqzO3j1FTBt1ks8zrs9RR5uceeJyM0EWO/KP4zVkDPZH8XGbyGqSEyyprpNpy9RWqSS/eUB8R32R7lh6N8Ue8qBStE7DtjFlKn9tV2/YRCOnewT5/z3X/W8NbpivDkpHb2KE2wF4+NIZF+UodImHJDdUjssQbqssyX6UZxb0m71LmnLDy6rSVLEPxjNujNpVHiTEfkc7psk8nh/JG/f5q//eZooYs0016KYFRiYaPo0JGSXw0rB/EkAfihldxVA4DFYScKKFXDX9U7mp8UWxShmAwWNkVCoHHUH1ViJjFevbtTKWac6oMMM6+naHkwrZGu8NwVUOX2taaV86FV/9k3VDZeTrOGHfXu8kTXp6k7U3YPT4um/5I/QzIGWBfHLW0L0bT1+C4o/wqGWP/24kSroH5RBa+bQSMgBEwAleEQHTeGityJD37oeeayzi5VhJ7GDZjxLMuaE/oK8xG7WP57Uk7hZs9kPy8fpjcZ5WZKd3oXJruDa41yE8Og13e/yE2Oo9H8MlqVpUnJxrxM0aUeYwyD2tYXhPkVCFx4My6v//qoKIHacu1fg1+c9ZA4NOsNQdLtl1kTbmGtIimj7apaHtGl2jeQ3oTR7/fO7Ck9hbp843MGt4mweGfsfLkVLlPTI5hSoA9Z5JP+rdjXc5kwo/WI2Jq6nIXwLIm74liNLdX72/ZcMdEkq3CpSJY3Ic2pCUP4e9GKgYjE40EsA5aenX/hVzDK/Y1TpWRLZQb4f3VGf5GsRFfMxEXOM0tDRNynsTmypiT7uzbmQr1SvjWnioDi7NvZyipsoX6ofj22pri+gsbDHks8g4ZWS61rd2pbHuTreIY1KE8VhFuxzkuIGHCZI2LPrUkPwi8criJ2OhnaV+Mpq/BkaJmfKcmuzUw3whaizUCRsAIGIEzRCA6b42pHknP2qBcHzSytCbgDxBYU7TGOW6ktQKb3/7m9nPdLtdck/ss5BU0OpcqL9aBOAyJQ3uBx5JTliE6j0fwySrWlien2/NVhrxPZe/YnBLMTLrHWgFCN4jrco3fRFb8bLnWr8FvzhroUtfjUVyibSranmkWobzVzuiv38t9pDDyd/LpR9RTZx+gePYGkT6PjDCvZEaIE6e5L5T86MiYkPt99lse3Rsaw6J8rZxAIFqPiArXZRDLmrynirJJf3tIrioMk8M7+TS8hhSmwfFK3cv7mPtGqPj3co2FM8fL/1ZxNIaWdP11uignKfIprxsW8WI0LDdqNPSDJ/OahLEfjArPA6wfJh6APgVFsAGXQ0aSqN500DXkRPNr+FTPl9DO5pwqo3yX0s7QdXFbU10ymTBedCYkhPfoUtsai+jOGKQyU17GRV7VZ7LYyeeaJ7btuJHimBjhi/azP8QPIW9zkl6L+mI0vQoSwrEocC5/g28R3wYTvmtg3sp04PYQoB3JcdICx0YOl9vf7QHiEhuBK0dA/ftY894QkowtnfFF+jCPsaH/F3AAAB69SURBVFZgLGr2UdlXXLvOUBiK7rPuuf/Ka2wuZS+2N+Ypf14j41XS0qAXmseVJoRvUrC2PLlcY/6nuvGldGBtWhL7S/5gLBsYwIMxn0+8NK8Fyu/US5l4ILx4rU9+cnt7acWF8Ev6zlkDXeR6PIqL6irapkLtmbqvyJv6eEeagvb2AZIX7vM1vEWeH6Zwuycp7lEeDNs8QGgPaylM++/Ye8o0RRi+SF+J8B3SM1qPqBaqS5UxintN3gU0o8FyHB1lqrnx4P37++9pp4qjIu/k2ERi7PpG8Z2Np6759gGV3nldU9eAAoAQDYYGzN9rI68lXefGghyICsbC2fDJ55q/JR78XpDi2Twz6B76IL5Y7kl8DOJMToP8is+WZ4xi5E15AZqJi0H0aKT8prDBOIpFPOPc0U3xIWzEB/YsCvJE1pFTe5HyHX0fu5QnXjA++3aWdZa+6Mr3FgfbT8F3Me0MnVO5CI71w9G2prT0de5/oXDTbxE0RLp/yW2NMaF80EC5GXvafqMwfY4xhPZBGCLM4N8ZO5s7Iz/iZbEMpnwsOpxuRFwoWvks7YvR9JM4ZoWlExjzYIjvSnaenBc8q2Ce5dm/TQTUvuij7TqGayHxQv7Bsf420XKpjcD5IKA+yhwQWnP2tVba6LzF2mVorxNKn/NN4wrzGvMgxLz2RvEYk8iDe0PEKTIMLS3pGt68/h/dZ5Eg8U7NpaxbX8uxX8vEfqld4+RIxYXmcfGF8Uk6hsqDHuLP+mYs0ZN9ZrM5LuTdJb3Bi3G+Ux5dgwuySkIGvKPf6oZZ90m3eE8pObPbl9LOWgOlPFfb+4FHlJR3qF0swQVdlD7UR8QXas9JZlT3vIeHH9rrS6l8oT5fyVtly5BsdM16PlF4z95DASDxsjZC59zv2jGM+5kifOIJ6Sm+UD2St3gn61I84bG2Ju9c9iFfcjCwYd9q15lDfP048R+c31rDXT/hqa6Twq/ldyasrE+6HzbckU5pqDCMDEfZEGdd1/ZTOZoTYQrf9eUrjso+iI146AwYJVfbnKR8Zy2i+mU41nXSebSdoYd4mglaQYwpe3j3dRXPVbSzVHbKstfWVMZmMJXfDETpeid/76lCunfVbU1lZEJ7Jn9wvOq3kbFrpc8LMf7ZrbPIHEtzq/FrYX6r+J263Kq/qnG11Fdp88KYh4vMYcx3e2NPmWYsrHQsqtr+pmvGNsa9oxnPx3QjXvqAE5usvUWf4hbhsDT9Ib2PeW8NjJZisTT9MfG6lryE+cGNzbWU0+VYHwG1Hd6kwEDXPBiUj/ECwwXjLeNqxHh30rW+dK5edyrN6ns/YWUyAkZgAIE0rvBqMQ9o8sOJAc79KPEfnN8e7ic5ecze9+1W0IjBeG/xu4LcY4vAOMamntM5c4kGlBf9c2X002HUmjRs9ROd+PpgOxPOTHLgFDLapbJcSzujOHttLWHCgoGFDa8ZsNChPZVPa3XZ0i20NZ70rGFoy8YH2p3pMAJrYX44F99dDQGNFbwaxNNnxg++A1PdzpWWUxJv5HPCglMRjC88GKiWpXRQ82DiPniWvzyZZkPZoaU4LE3fUeb0F4swWorF0vSnh+9iNbjENefFgn0tiqu/NieG5Len+RXmLab8hxys4dkbTNGp1/pz1kBbrMencPJ9I3CrCOR1KQ+Za+ng/PawVtoR+BmQWIw/XisvyeKVVwwRGci1RB9bDkcu+XZg5+Qg5ZJjInkt14QTX0c/+NL9VV8BllxeJ1hVZkfxbS5G21nCic0hpzGY1ENtR3zX0s5AfKitsXEGN/zsXoGRrjsEZooAt1XbheSdTVuTLoxRlPNNp/AzLiQr9+lFJ/dmZH1RSdbE/KIKfuHKqt4YRzn1jrHsh9riKB0PqzD+dTZciuMaY2A1SVb7iY6UGN2m/mynOp85CaRb/kZwJ/lSHJam7ygz40L5v0g6zEjdTSI5izBaisXS9N3S+KoGAWF/NuuAGr3Ne3IEfpcGrEt56DxEzAGT84nSn2ytr7yr151Ks8l6fAhAxxkBI9AgQJ+D8qGM+6vAr/rrwfntUUDGsVkYVFGaRdnBbw1UKtZ8iFVpIk9TKkUfh12YMKnskeJpGDxNmSImpEEZUwmv8P5gOxOWdDZwwnCXJ3cwi+ALTBffziiEyr7XThT3AfeCdAtt7VnCYo0Td4iiH2eZSbS9HgIZn7Uw74n35ZkiwLiajdulihjNv9bYhFFv7wFCyXgorLSM9byecHLDedKFsgyVZykOS9MfgjFyj00nbhGthNFSLJamX4SBExsBI1CHAHOEHGM8n8n5Sn4+DfOhwoxLfFt3aJ7RrT2i/3Pi99h7yjlroFtYj+9VkCOMwAkRyIa76HgSVvVhmPN4jEMnfXYaTDG0TJ4qG1NTaVkE872Ywae0Y+muJT6Vm/JXW3+vBYNeOQbbmXjCp8p68prLW29ngHBDbY2nt82pzKbyl/8wwOfBfrm065SwNubXidL1lYqTvu8GipXnM+7PIo1X9DnWFjWfRZiVVzDRV9Jp7KTyUhyWpg8WYXO2NTBaisXS9JuD5AyMgBHoIqCxFeMdD+ZxfH+q+QaVfPZH4U22eE+1p6xaA0lP9rze+3Wbga+MwNYIPCcD9b+8Rl0tv0erSVpJkAq5d9IH0anw0VNPg9okGWue4hvM5xwjVfabLPdYXQiPsXZWc6psUPwttzMAuZW2pnKycFvz5BcfTeYfjBedHqIOrpU2wPxaobqactEfAoXZ+xYcaZQW4wonI+irGOjYrLULKYVzXHNqIl3vSh6lORopXzZZnI7YI92bjQPClqRX2oM47im7YYR0WYzREiyWYrkhNBZtBIzAERHQOMJcctS9lfKsWneK/6j6HRF+Z2UEzhkB1kxr7g/bsp6d4a7VzAEjYASMwG0hwJNUiAH/5ybkHyNgBLJRjg3LGO0ZtbRh4RQdRvDmIY18+hUnqpuHM7rGaLfkswhKvh4lfTgN0hoWe9Jn4VDImJVe+hzEsZC/eXBFjGZhURRwafpClINGwAgYASNgBIzANSCQ1imsSTmMsTo9XF2iBRoBI2AEjEA1AhrseU0D4wTfXjEZASMQR4BvFLWkvsT36ji9Wp6sxjCOIQ8DHrToswj3Ilb95XWmsVdkoxl1cIgmKvg66YM4Fsk3Dx4Tow4WM0q2NP2MLJ3ECBgBI2AEjIAROCECeY25yQGMRycsmLM2AkbACBiBLgI/6vLLbpSvjMBNIzD0bbsMSD75lD8ynuO/V6D/aQ1O2EHN6TwZpWZ9FkHpSP9bloPAAPGPuqPfT9I9/jV38BXZQvYcHIrkg98IzPdn45gFlL7KQ1ny4rW81eSj+6VBNd9/q3heax6klTE6BZaD5XKkETACRsAIGAEjcDUI8OmV/2nNMvb2xKKC2nC3CD4nNgJGwAisigD/UvZKA/4LuU2+j7CqthZmBDZGQP2A10fJpTG49bLLce0Cib6TeDGCl5QNSaMGtJJ5LIw+uvd07H5tvORhUOQkYFuGIRnkK8etXOaSLceNyqhNL/7ZOCrtkGFup3hOQvJHY1XfXRL/qhhJ3lGxLCvKYSNgBIyAETACRuBqEWCd039wvFphH60myYKMgBEwAkZgEQLaUP4qx+abkyc23C1C04mvCAH6Qj4xVxYrnxQr+wp9h9NbGNhK4hX0zZ6ClhlVhinXc+mL0b6kT3SBkYt49GYhWINDKSuHa9KfE45bYFSDRcav9JemL2U5bASMgBEwAkbACFwwAlqn8fYEtPSzJ/dSBn5tuBsAxVFGwAgYgRMiwMfgv2WjLtc3PpxQLWdtBE6GAMYr+kWfOPnWN9JxUgyjSkvqRxjBcKudlGuFLwxIN3Tt6ItIxf9JvPzy9dEaHBDTp5r0Z4PjRhjVYNHHkeul6YdkOs4IGAEjYASMgBG4TAR4wPqz1iyb7d0eXiYu1toIGAEjcJ0IaMDnSQ2D/uvrLKFLdeMI5I/259NyLRxq+7wy+l6OP45oKfWJd/J5BaEhhR8rwPcgX97HNMYu4jid1cpOfBhZ+GODRa/J5nyO5FMWXEvSn7EhhIN4bwHHJRgdBcu28hwwAkbACBgBI2AErhIBrbl42Mn6c7PXZAHOJ+5AwWQEjIAROC8EGPg5dfeN3N15qWZtjEA9AmrHGM8gFjfQT4rjtfBf5DevFdDWUxz/ANsnTsv9S/efy+fPKPA/1XVpjMuyOa36te5DH8sd/HOIhutMfqQ3f+zA4g/iW5fg9kZ+/i7cJA7ivWoc18CoQff+BObBNrUClikre0bACBgBI2AEjMCVIpBP241+a3iNcj94//79GnIswwgYASNgBFZEQBtGTh3xqtymT29WVNmijMBJEUgGnWfyz+6V2JMCU5n5VjhK7qw/p6hU3+xGwAgYASNgBIyAETgKAlrb8ND4F7kPFN70sMXDo5TImRgBI2AEjEAtArwC+EqTQD59U5ve/Ebg1hBg8bT3vbhbA2GF8m6FIwvaTRe1K5TdIoyAETACRsAIGAEjEEWANyV4s2Pz9c2jqEbmMwJGwAgYgeMhoAmAj+5/oxyZED47Xs7OyQhcHgLqK4+lNUbuN5en/flovCWOkm2j6vlUtTUxAkbACBgBI2AEFiCgdQ1/nMZ+7ecFYsJJfeIuDJUZjYARMALHRUATAd+14ntVQ/+oeVxlnJsROG8EniX1bBxaVk/GcRl+Tm0EjIARMAJGwAhcOQLam/F2At8i/uJYRbXh7lhIOx8jYASMwAwE0oTAxND+o+YMMU5iBK4dAf7Q4jP1k81fVbhyII3jlVewi2cEjIARMAJGwAjMR0BrzU+UmkMVn86XUp/Sf05Rj5lTGAEjYASOjoAmCT58yr9llv+ieXQ9nKERMAJGwAgYASNgBIyAETACRuDWENA+jE+z/Cb3qcJHfVhsw92ttTaX1wgYASNgBIyAETACRsAIGAEjYASMgBEwAkbgIhC46T+nSBbT16mm+Kg19PLY1tP7bG/z13Vw/vXuOjr/OrKGRsAIGAEjYASMgBEwAkbACBgBI3CdCNy04U5V+i8ZJf6eq1Zh/r3xP3If5zj7myPgOtgc4sUZuI4WQ2gBRsAIGAEjYASMgBEwAkbACBgBI2AE6hG49T+neCVjHf8IkomPDP5NcXxw0HQcBFwHx8F5SS6uoyXoOa0RMAJGwAgYASNgBIyAETACRsAIGIGZCNy64Y7TdvyDmul0CLgOTod9NGfXURQp8xkBI2AEjIARMAJGwAgYASNgBIyAEVgRgZP+OYVOtvGvHJx44/tyv+r6pP+WqPw5cfe5/Jt5VdZ1oBq/Mtq6Tm+xn1xZE3FxjIARMAJGwAgYASNgBIyAETACRuBCEFjtxJ0283uvl2JAkMMot0eKx2DHH0NgrPtO7pnifpJ/ElLe6P+53NOTKLBCpqkMHUmKcx10ELmsi3Or02voJ5fQAmrr/RLKZB2NgBEwAkbACBgBI2AEjIARMAJGoB6Bzp9TaLPIiTPoDzlOnfFR+v8REaDfxMsJunxqjjC0ZwhLfJ/J/8c9S/P7HfFyP8l9UcRvHlR+GBcp+1OF7zbP8EAGyt91cNl1sIvWYZDvbPqV9L2WfrJ2HYXkCT/GxPwv1owCXDPe/cpFj0L1rrTUSR5Hnyn8jmvF53G4Fau4qrGl4M8yfpiQO8inNGEdESB+/iQo0xMFJv/pW2l48PJ3+e2fDWUBigvjHslbPJQn54Pspl8ofq8eE+9k/UT5cpnsGwEjYASMgBEwAkbACBgBI3A7CLSvymrjwL+pfiP/Z4ovnw0JcRjYJo134vmveNlkkQ5+5CDvTn6HFMcG8pXcR+V9hUn7pxwGtL2Np+JXJ+XTbCrlNxuxdL2TP1nmtZVRnq4DgSocqBP8S6yDUB2qbFG+s+hXqU4wCF1DP4livzbftxm/PHbomlPGGMSacbeIn6x3paGfIPOzIh1j69dyjNutIUnhUFmQI17kohf13ciQz/Un8tvPCET4Ek9UxzznwP9v5beTj0HuNznmhNHxQPfA6638vYc+ipvEXTyhvBNf/1+WOakNPl/ofluPCofqJ8on+SYjYASMgBEwAkbACBgBI2AEbhCBh5RZGweMaJx2KzcdGNy4Lk8/6HKU2DR9IPdA7mM5Nn3IGCI2WT/276dr0pT/9DqUfpU45ddsrCSMjR2bUjaJnI7g1Eo1Kf0LObCsppTOdXDBdRCtwyhfakQn71fS92r6SRT7jfgwnvXppSLKU3j5fqTeMdLlk19NOunN+MUYiiGpoWhZMr980mJMbA1/un4s1zecRfhCOko29L3cE+XbGO2IUJgHOL/Ljc5D4sFQOUi6x3gcwT2aN/L4l2WMdZkyTv16jJY9ypfzs28EjIARMAJGwAgYASNgBIzADSHQGO5UXk4pDJ1we6N4jFFs2lYjyeO12M6GsxCO0ex5cb1lkA0dRkL87NiUsfGdQ+A0FyvXweXXQbQOo3xVbXDDfnVN/SSK/dp8nFRrT8ZVVew4M2PXf1Xv/TEHQ1L5bctoWXaShUGKBxh8d7QlxXOCr9U/yicBUR3Ji7z7xkHimZsG5yHpga6M17ghiuIezRtdOvlJh7G8o2WP8g2Vz3FGwAgYASNgBIyAETACRsAIXDkCj1L52Dh0NmopPm+iuN+exkv3Fnna7HByge/f5U3PG8WRB6d7cr4KbkfK74PtpFdLdh1UQ7Z6gqV1EE0f5asu4Bb96sr6SRT7tfl4CPK1sPxFPq9U5nGP01ajp8kmGgAGOk4KZ1l99scpIloW2HmgcndAZhIZ5gvpqPyyrjy46dMfKeKZ/Hy6LfN8pbTtK9w5svAnca/JW7zk35k3FIfRD+rXY6jsShflazLxjxEwAkbACBgBI2AEjIARMAK3hcCjYtNyqORPDt3M9yQLYxwbsA/lMMDxjTtOKHRIcWxe/ye/PXWnMCcqSA8NGu50n01u3sSxWf1dcXvyEXBJpDLkTeshtV0Hh9BZeG9pHUTTR/nK4ijNpv2KvJQHxh1OZtGv6LsYQwb7oe6dhKTPon4STb82H2BJ5s84BTHy/Kkwr7RyGuyXFK9glxR/sN51n/oaIk6h7XSf121rMcM4xtiMjK/kGG/Rs/8nGiE+yZnUUbJ34sNYSHBonGM+gWiXLYmfV2T7xrL2PgHxTOIunuq8cyZKS79pXnVVuPPwS9fRsof4cp72jYARMAJGwAgYASNgBIyAEbgtBB6quHmjxIZ9jCKbP3j4bt2/5diU4v6jMBublnSN0W4nvzXapWtOHeQ4Xs1rSbx/k+O7eJzKQz7fQOK0HhumayDXwelrcWkdRNNH+TIim/WrnIH6E/2IU2D8Kyf9FmNIpw9m3hP7tdj11Y2mX5uv0UPYYqDJxh0wZ2wce/AQqvdGcPGjPD7RZWN4TdHRsmQpeax/RluQY7xlXMZwl0+WwRvly3JbX3L6OuZ7+cR1vs4+/FDOcycZlBGD26RxWTwR3MN5o4hkctIRwyHYUIe/y00S6cRU1s9gmijfYGJHGgEjYASMgBEwAkbACBgBI3BVCDwKliafeBhl10aj/f4RTLrm1AbGOIwAnNggjhMkbFYxug1R3pj92LuJsS+fnMi3vlTgGy4kl3Tfy32u8APioBQPHwaJsTwb3gv4ubg6SPi/TtiyWYVeKv7uPnhxv5N1MFGiaPqWT1ht2a92ko8xhn7T9NGkP0aIx4p7IUcfviRqsZupdDR9NV/CmrYP1oyLjIV8o67zb6SK2ylust7hGyD+MIKxsv2DhwGeflRTFqXJ4y9GqWxgzLyMyd8rvmwPB/nEO9bPx3R8qTz+T+nadqcwhq4spzTSZSNz1m/Ulwza+BTuNXnvJBNjXWN0TfJ5SLVXjwNKjZW9zxrl66fztREwAkbACBgBI2AEjIARMAJXhgCGu3cHyvQk3cuvpx5gHbzFRotNGCfmCHOah1MSzYZnIAXGHTadeaO2U5gNLvGNkS6nUXz7nSH45X5IfA2LrtnwsTGGHt976/xKdt509wU2eOk+pzD6xGtrY69EXWUdCIB/lVgk3DjJVRqJ+jiFrpOsXL9lmlPVQbQOo3xlmfrhxf2qEIjBm35ZEv0NWtRvbqCOwnUpLHho8VR+Hhv4swcMShhoGoOYrttxT3FD1K/3Do/SMy7xwKQcZ8I6FsLIp0/0W8rwTC6fLpviK418jbwRHfM9xvGPdMG4wfiNQfFNcmDV5Kd76EFZJynxTuIuvlDeQxkqbZ6zOJXIP6sP1qPih+pnT2SUby+hI4yAETACRsAIGAEjYASMgBG4SgQescmQo3BDm/QcN7RBawFRek7EPZH/tI3sBh7rHrIaw1z31v2V7mcjTH9Dxom5X3V/cDNUyOKESrtRFP9bXWMsY8O3Kklm3nx35Ka8MFLWnHbZif8q60DgvFLZym9j8XogcZzUoX5mk9KfVR1In1AdRvkARrxb9ivk0+fol/0Trrkvuo6oiC61Y6LwC9V5Sk7bxyjVktJj8MGIzWcAwJzXNXeKm6x3+EpSGoxZjMH9k3phHZU28x4aaxnfGI/J/iAfDCUpzaCOPR5kdvq20oEdhFGSOeSx/INzUsN9/1OD+8G8Ead8MSji9/sGxkzqsK1H+DKJf7Ls8Eb5slz7RsAIGAEjYASMgBEwAkbACFw/Ag9TETF45VM2Zamb00uKaA1i5c0izCmMx8V1DubTT2xysixOUAwRp34wtLV5KYxMHBvZKWLDFOGbknOq+9dYB2zA8+mcU+Fak+/SOoimj/Jt0q8KQDiZRZ/DYFHSV7rASBI1jpRptw5HsRvTI5p+Nb48jg3gvEsYY7DL4yN6R+q9LZ9k8HDiY/ntSTuFMbDlMT1aFmTCOzSWcw/KbSLK1yQK6Njwjfx8ovj88IYyPZc8Hgi0TnHwvEhxjaFP4Wb+kN9v3zvFUY4+7oraozJvbnLykNdiD2HUESLeqfpp+KN8HeG+MAJGwAgYASNgBIyAETACRuDqEXiYSsjrWmwW+8QJuqGNfZ/vO206hl5/xJjGBm+n+2yU9jZQ6R4bG3jbjWeKz/zZJ7olySQNstlENSdB2puXF7i6OlC90C7KusOQh0Gof1rlXGpraR1E00f5NulXBdj0n45hVXWDoQLX6YtFmlMHo9iN6RlNvxpf6gOcZsuGtL5ujF/NOJluTNZ7FpDqC0NW/3VnxtR3iS9aFtg58TykJ3MBZch6Rvl2QR3h41uL/OMueDSUwrTTpny6xoDHt+Q6TvcZZ/K9zEtcCHfJm8xbsiBk5ldjm4j0k+fPjE8TLbn0pan62UX5ygwdNgJGwAgYASNgBIyAETACRuA2EGgMd9o08CHyd/LZ7DWkMJsnXlPlo90NESf3Xo5TByV9q7jOK666/joxlAYANlT916DY2Hwvx3eI8mmOlLTx0K2U0USKl1MVeWPKxu5cjUGNvlM/Ks9V14HKRz3TvsZep56CaPP7S+ugIn2orlXgzfqVdKV/Y6B5koFNcRh5+PD/WfanCowHx6qK9KE6isoTpoxhnBAD95Z0zSuUv8gvx75QvSsN9Ud9UdYmTfYVRx1iZNrJD5Ul8XIKDQNYfj2V9ENzQZQvpCN5i+DNY3oToZ9oe0THDrZJQBT3aN7MYZ2T3cKHcY28W8zJW/Ghskf5kGkyAkbACBgBI2AEjIARMAJG4PYQePD+/fum1No8sPFgs3Yn94fcc7lvFN/ZwOua7zGxsesb4NikNCcd5GMMYAP2D/EhryVds8n5Sq7coO3xtQkUUJq8iSRvCF3b01y6j9GQkxU5f3gaUhz58bHzoROBiWsdL+VV/Y27nLvSX2sd0DaoI07JdNpDLvta/hnUQbQOo3yb9KuEE0YR+gV9BCKMAajT55s7K/5cYR1F6xLj9Wu5cuwrvwHZoiyMJutdPIyH8A0RJ6VbI7nCIR2zIPEz5pIGYjzfmwu4McWn+2EdC3kEc96D+MAAST7jChg0p6/lY1B8o/j2O6MKh3AXX55nDuYtPvIqHyaRP3NM/7RdqOxKF+JTHiYjYASMgBEwAkbACBgBI2AEbhCB1nB3yWVPG5/me2oK35Vl0fXFGO5KvS8tPFYHim8MEPIbQ2+63skvTxitVlzJpb5nG09XU+TMBQknDB7P5LfGnWOp7Do6FtLOxwgYASNgBIyAETACRsAIGAEjYAQuHYGHl16ApH/+vh2vnZ2SMBp2DIenVObIee/VgQw0GO0wEHGKi3+S5eQLpyLLE0e6XJVuuQ5qgOTUUOeEUE3ihbyuo4UAOrkRMAJGwAgYASNgBIyAETACRsAI3AYC13LiDuMQrxvxelhjjEhGI055YaDAYMSrU38ovn2FStemlRAQrkN18KfEP+5nId4H/ThfHw8B4U+dUDe8usyrhSYjYASMgBEwAkbACBgBI2AEjIARMAJG4AwRuArD3RniapWMwNkiIGMdxmw+sP+Bwpx+MxkBI2AEjIARMAJGwAgYASNgBIyAETACZ4jAtbwqe4bQWiUjcLYI/C7NPrPR7mzrx4oZASNgBIyAETACRsAIGAEjYASMgBFoEPh/O++PmfV77uUAAAAASUVORK5CYII=\n",
206
      "text/latex": [
207
       "$\\displaystyle 0.5 c_{l1}^{2} \\left(1 - c_{l1}\\right)^{2} + 0.5 c_{l2}^{2} \\left(1 - c_{l2}\\right)^{2} + 0.3 \\rho \\left(- 0.037 \\rho - \\frac{1.0 \\left(1.7031322 \\rho - 51.0939660000001\\right)}{1.0 \\rho^{2} - 40.0 \\rho + 400.0} + 0.042578305 \\log{\\left(1.0 \\rho \\right)} - 0.0530922164415325\\right) + 0.5 {\\partial c_{l1}}^{2} + 0.5 {\\partial c_{l2}}^{2} + 0.005 {\\partial \\rho}^{2} + 0.00085206629489322$"
208
209
      ],
      "text/plain": [
Markus Holzer's avatar
Markus Holzer committed
210
       "       2          2          2          2         ⎛           1.7031322⋅ρ - 51\n",
211
       "0.5⋅cₗ₁ ⋅(1 - cₗ₁)  + 0.5⋅cₗ₂ ⋅(1 - cₗ₂)  + 0.3⋅ρ⋅⎜-0.037⋅ρ - ────────────────\n",
Markus Holzer's avatar
Markus Holzer committed
212
213
       "                                                  ⎜                   2       \n",
       "                                                  ⎝              1.0⋅ρ  - 40.0\n",
214
       "\n",
Markus Holzer's avatar
Markus Holzer committed
215
216
217
218
       ".0939660000001                                              ⎞              2  \n",
       "────────────── + 0.042578305⋅log(1.0⋅ρ) - 0.0530922164415325⎟ + 0.5⋅D(c_l1)  +\n",
       "                                                            ⎟                 \n",
       "⋅ρ + 400.0                                                  ⎠                 \n",
219
       "\n",
Markus Holzer's avatar
Markus Holzer committed
220
221
222
223
       "            2               2                      \n",
       " 0.5⋅D(c_l2)  + 0.005⋅D(rho)  + 0.00085206629489322\n",
       "                                                   \n",
       "                                                   "
224
225
      ]
     },
Markus Holzer's avatar
Markus Holzer committed
226
     "execution_count": 7,
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "free_energy = free_energy_high_density_ratio(eos, ρ, ρ_g, ρ_l, c_l1, c_l2, λ, κ)\n",
    "free_energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the free energy expressed in the order parameters $\\rho, c_{l1}, c_{l2}$. Next we have to transform it into coordinates $\\rho, \\phi$. "
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
245
   "execution_count": 8,
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
   "metadata": {},
   "outputs": [],
   "source": [
    "transformation_eqs = [ c_l1 - (1 + φ/χ - (ρ - ρ_l)/(ρ_g - ρ_l)) / 2,\n",
    "                       c_l2 - (1 - φ/χ - (ρ - ρ_l)/(ρ_g - ρ_l)) / 2]\n",
    "transform_forward_substitutions = sp.solve(transformation_eqs, [c_l1, c_l2])\n",
    "transform_backward_substitutions = sp.solve(transformation_eqs, [ρ, φ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To do the transformation, we use the substitutions dict.\n",
    "After the substitutions the differentials have to be expanded again."
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
265
   "execution_count": 9,
266
267
268
269
   "metadata": {},
   "outputs": [
    {
     "data": {
270
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAADMAAAAVCAYAAADrVNYBAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADA0lEQVRYCc2W61FTQRSAgUkB0Q6EDjBUoHQQtAJjBzL+Cv8c6ACoQKADtAKFDqQDkQ7i9+3sLpe9ezM3IU48M5tz9rz2vHZvNmez2UYJR0dHQ3gj+dDfSvm69sSyzdm7rFvouzKOrZKB0id43yP/Zylf5z4mYBLH0Ncsi/4Idiat6XQ6Zs1Yw8Trg9Gf9NFbpQ5nXrMumz7LzuyR5h0ZPzym24va6aW1WqVb3DlyGcpksmBB4veC+v9E/dnJ0MW3RGaV1g7PToYMdknov3jxBsuUk+Cd1fcsx2ufPWjDe3MIveh90zYD9nb6gKWfbZY+W88w/BaUnfGpmxsMjk/R+Qj2kBNon0jxMesG2gCWAmz1cQAO/qE966bDmYV82ZTlZGIQ7xDqoAroXCrwsIitYhgxeFZP2oAWBuzHGI2T7+jA79wQnueUcAVDmXYBBmyspAGIrUp1/uFPkOv0FSuBI3aYNmC7Wju4odJJniNp+lIxddmJeQKc6yfkBcxz8Gfwly2ZEiyr4Oi0DOELHnSBfN4Yepd6zbcOE+DTAnjuReJFnArT9VqahAkb99WWRhD+1wmjw9YK1UCjPL/oe/ivpBj3Ht45pkm3gr3wxlAWykfGDrQKBM/x8q/XG+gwTSEZGAlaX9UkACszgQQG3hxJi3CGY2c5AHRTP7FrWF9ORgZs7bLLRGvgv5UH9HIBymRqRon3ASJ1T94ejkLFwD4M9+AshzaRP+DcTY1KiHp2Pb9MkadPX7WuEdPVvT8JNv2jlgBDHwJflOp/LfhWyoAdr31W6tZp7VB4aQxfQ+cKYpcBvuNi4J6ZXibpqk/4AbBrxTpIwj4YBwZvtaykLT6bZ4d8h5UC7FINRUHPLvu9WhoWGbPmIXaoeV+aspJ2HKtdiYrl3Svte+/LZFpf1Q5P+b50yAObJLw3+qxClNvlH1WF+Ux9PylSmUzrq9rhrzPAQn9CwPNGZxT1+3Y5qMcitP6tPHkA1ETRSvmB9HX5yj4/texXCjGoEbhXMjE2HyBjbMX2FyjLoHFGrT5SAAAAAElFTkSuQmCC\n",
271
      "text/latex": [
272
       "$\\displaystyle \\left\\{\\phi, \\rho\\right\\}$"
273
274
275
276
277
      ],
      "text/plain": [
       "{φ, ρ}"
      ]
     },
Markus Holzer's avatar
Markus Holzer committed
278
     "execution_count": 9,
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "free_energy_transformed = free_energy.subs(transform_forward_substitutions)\n",
    "free_energy_transformed = expand_diff_full(free_energy_transformed, functions=(ρ, φ))\n",
    "free_energy_transformed.atoms(sp.Symbol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the free energy depends only on ρ and φ. This transformed form is later used to derive expressions for the chemical potential, pressure tensor and force computations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Data setup"
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
305
   "execution_count": 10,
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
   "metadata": {},
   "outputs": [],
   "source": [
    "dh = create_data_handling(domain_size, periodicity=True, default_target=target)\n",
    "\n",
    "# Fields for order parameters\n",
    "ρ_field = dh.add_array(\"rho\")\n",
    "φ_field = dh.add_array(\"phi\")\n",
    "c_field = dh.add_array(\"c\", values_per_cell=2)\n",
    "\n",
    "# Chemical potential, pressure tensor, forces and velocities\n",
    "μ_phi_field = dh.add_array(\"mu_phi\", latex_name=r\"\\mu_{\\phi}\")\n",
    "pbs_field = dh.add_array(\"pbs\")\n",
    "pressure_tensor_field = dh.add_array(\"p\", len(symmetric_tensor_linearization(dh.dim)))\n",
    "force_field = dh.add_array(\"force\", values_per_cell=dh.dim, latex_name=\"F\")\n",
    "vel_field = dh.add_array(\"velocity\", values_per_cell=dh.dim)\n",
    "\n",
    "# PDF fields for lattice Boltzmann schemes\n",
    "pdf_src_rho = dh.add_array(\"pdf_src_rho\", values_per_cell=len(stencil))\n",
    "pdf_dst_rho = dh.add_array_like(\"pdf_dst_rho\", \"pdf_src_rho\")\n",
    "\n",
    "pdf_src_phi = dh.add_array(\"pdf_src_phi\", values_per_cell=len(stencil))\n",
    "pdf_dst_phi = dh.add_array_like(\"pdf_dst_phi\", \"pdf_src_phi\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3a: Compute kernels and time loop\n",
    "\n",
    "We define one function that takes an expression with derivative objects in it, substitutes the spatial derivatives with finite differences using the strategy defined in the `fd_discretization` function and compiles a kernel from it."
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
342
   "execution_count": 11,
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_kernel(assignments):\n",
    "    # assignments may be using the symbols ρ and φ\n",
    "    # these is substituted with the access to the corresponding fields here\n",
    "    field_substitutions = {\n",
    "        ρ: ρ_field.center,\n",
    "        φ: φ_field.center\n",
    "    }\n",
    "    \n",
    "    processed_assignments = []\n",
    "    for a in assignments:\n",
    "        new_rhs = a.rhs.subs(field_substitutions)\n",
    "        \n",
    "        # ∂∂f representing the laplacian of f is replaced by the explicit carteisan form\n",
    "        # ∂_0 ∂_0 f + ∂_1 ∂_1 f     (example for 2D)\n",
    "        # otherwise the discretization would not do the correct thing\n",
    "        new_rhs = replace_generic_laplacian(new_rhs)\n",
    "        \n",
    "        # Next the \"∂\" objects are replaced using finite differences\n",
    "        new_rhs = discretize_spatial(new_rhs, dx=1, stencil=fd_discretization)\n",
    "        processed_assignments.append(Assignment(a.lhs, new_rhs))\n",
    "        \n",
    "    return create_kernel(processed_assignments, target=target, cpu_openmp=threads).compile()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Chemical Potential\n",
    "\n",
    "In the next cell the kernel to compute the chemical potential is created. First an analytic expression for μ is obtained using the free energy, which is then passed to the discretization function above to create a kernel from it. We only have to store the chemical potential of the φ coordinate explicitly, which enters the Cahn-Hilliard lattice Boltzmann for φ."
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
381
   "execution_count": 12,
382
383
384
385
   "metadata": {},
   "outputs": [
    {
     "data": {
386
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArQAAAAyCAYAAACtSEg2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2d7ZUcNRaG23MmADNEAGQAJoI1GSzrCGAzgON//PNZMgBHYCADdiMwdgaQgc1kMPs+al2NqlpVdVXd0z3tuTpHLZXqfunV1y2VuvvRzc3NJkIgcEoEfvjhh0+l/6nie0XyXyq+UPlbpRECgUDggSGgsf9YVX6eq82cQPhG5dfbbHwGAoFAIDBE4HJ4GVeBwEkQ+I+0/qXF6me0K/1Oyf8UP+I6QiAQCDw4BP6jeeDfVmvlf1L+jeJnVhZpIBAIBAI1Ahf1ReQDgRMh8EJ6X1W6P1ae3doIgUAg8DAR+FZOLG9tLPDQ+6nKPreCSAOBQCAQqBGIHdoajcifBAEtUuOjBf+UISxgEQKBQOBhIsDu7B8Ps+pR60AgEFiDwKM4Q7sGtuC5CwTk2HLU4JniK+V/vAsdIXOLgPC1c8u/KH+Qc4mS87ni4OFE15yFvFL6111iL/kHr89d2huy+xBQ+/KA+0+lceRAQDzk/q66n2ye6eu1QX1sBMKhPTbiWZ8Gpe1AvlMRkzRnxtyLvpf/0HRjuCSfV4D/VlrOu41p6mvR8SoxnZWtyy2vezhAvyr+rvxJnVrpP/s2Uh2a7aNydsHBeSpciyadYVaKs2jtS/twTX/9b82s6791zX1zaskTvtC96212++mVCbVoF9tBNO76eGV20oHJ9/AoPFHkyMz3ssuwoJy6uOgy7WK9Mx042xeoKOL6V+kat4+LDgEWJMNssCIeNgd1sht1KprZcV7TLuUliz5MX93pR0u8h75f4XHSeVt2dPV3cKhsN1h22rKima2f6FxzQtbL+WcLV8o0v9wnmV461zyTbVwckxgm2u6xkfma8yv36iD5s3S6b+NsEnfR9GDunmews9JvZu/0DbtRp+I72Div5a7Nh0O7Frk9+NQJ3oidb/H/hhilDCbKvlJ+0akVjYv/0HTYOg7S8afK3ir9enyvdS06HCGbZFokG91n8FNHFrDFxbMpZM9C6XVhPKXGy39ourE9kt9sH5UzgYJzq79xdhFn7DdF+iZtZg7tRnlbTL9WPvVh0VCOritFeJDLPfr5tdISdN0j09UOkumqD0aI1ivTS8fi8ZPkfmWVzPbwxoExnRxLpS46ZIjWpTvToru0Ty7DAWRRqtvHRZf5sRUZ9AOzn2t2xxZ3SUWzOM7RsxQkJ2EmOvraoB8t8R76vvS726Sl28vvoRNNT393taVHL/USnWv8Zjowo9+lzQmlzDl84Ze5Pc09XjrxpCD6xXlGNO6xhlDRu8fG1ortZ7Zlcf2bo9O9xX4lGhfmuS7uukuuq2/Uda7z4j/IOK9l7pO/2Ic5ePsRUAf4VlyPldYLDRM11/UTalO4l//QdC1jpIMFuze8qxkkAyz+VmSis2ALFztdq4LkPVUE6+6Q+c6+jVSP2fbRfZyttLtuqcBiofyvrq1/giFP4TixFmznr94V5B4T+0eKjxQ/U8QZsrY0XlKXTPF2jRXRL9bHK9NLlysFZmOHkoc26o4TaMFF16M707IgjsM3Kijt46WrhJhDbG3NrceKrQegiq1kB+O8lHZkZDOLLX2IduWNAV8Ko+zoQXq7+uLYQC+/lw75ol3s79mOxbbs0SuZrvErupeKHDcqb9qUZ4PiD8V6rfPSiS0FzzzjGmtIy3VfHENb1bef4pudX41yji7r9qw1XsxR6667aBf7htVjIt17nE/IbRYLr9l1/aLJFYUbOppi7WQdChV2Mlu7jq9VTmOxaMwFL/+h6QY2ZWxYsImuIB52/gZ1Vxn87xXrhRI6yn9RXBsei5G4Jnixm5Lt5T80XbHH0T70t1YY7MaKgPaiLYgp5DazyzWpV6YXH2zw1scr00uHbvrrn8Jl3N9wBlmszAnz0vXoZre07AwrPxW8dBvZy8MLc9/gaJDKcaAWdYlmZ5xPGTVVLhlghtPDzhm7wtjDQwJzxSlCT5u07PPye+lc/V24edvSq5e6eccvuut53XCBv17rvHTG70m9Yw1Z7rFhinN/HMyLdq9OHXRe3L2Yo95V946+UVep5MW/9zgvwvwZ5tjxPFu4L0ouMmME2HFpDcYxXe81naA1KZsu7s8FL/+h6cY2PVOHHix4Y4LGNQtTveNjJAzq57r3HVF5Fs1TnpfzYmf2j1Mv/6Hpajtm20c42w5s4VEZT/YvSoEyKmO3ll3X0m7KswAR6l2WbYnjs0OmF5+NZLrqI/O8Mr101Bhs+B1lFrhWsAnYS9ejG8eGL0tx5tz0YANtWbePlw5e5j52RKfqA81cmBrnczzje+yYgQOpxW/3sGksv/e6p01asr38Ljrh4O3v3rZ06aVi0r04J4jG+mJrrbNdvSdeuhagC2XesYaYnrFhamfnVyNSukTnwt2DeaXTW3dv36hED7KHGOcDgfteXNYCBBpPxW+UDn7QXteAjvHl1UHN96Hlc33nFqhVVZZcG+Rz/FdTN738h6Yb2yP5OJ31YjkmKdeiZWflmSKTGDs8StIT8ffKpwVTKU+fxJMH2XL2baQ6uNvHABcPY58xvnS+mbkgvdIS7c4DjcrSq1nR8FvCyOQM7Wzb6v6OTJXt2w479fHK9NKpbimIngeyVqDvb3Q/1d9DJ5queouec844NzxkcHSH9mPHCQe3OD1eOvERnigy/9VjF5k7XzSDmDCibY7zLeXtp3hod7BjHqC9mBPSg73SwRqkeycLsqWrTcaGevm9dGP5XIt3p79nusW23Edv1t0avzwQcbu1njE3EDhCgnNMfpYOgjqIZ3ae0X3XmESmaF1jyPSL3jW/LtHp/up+Jd4dzCv7vHVf7Bsm01LprecE1ziHN9vbHOsm+xDpxUgIHvv7URmXLGBMaB98EPB0MurLGbRDBxu0TOBTYa6Te/kPTVdsFT5MnExWtqNc7o0zosHp5YwmCxUPQ/bLBeDLgxOy7lvwYjdlt5f/0HTJnoypq31GFaBNiM0guTi7TOTMEThofzQI6bv8DNiPijhWRNqZyXcnqHxOphefHbm5oFUfr0wv3ZTuDXXTTfr30gPCmK5bt3SxUNjDBfUG752HCC+deG0OYgctjV2ltDsOLY7zIKise5yLBzv5oleaH5RHBjux9zF0t8moEl5+L91IfLps9XdueNpylV613dz4RTcPVK05/nNuKphtXrot15bPPc8YE/Yq3xyTuucdQ/Avzq+S56Hrxt2BuVV3kE7U3fC/s3GOEdJ9tLF+Maj1diJku7oEGUOl6Qi/l8KJjGh5ajrbIPupJ9/A/Eb56xNVxJ5e16r38q+lYwGyxXPSRtH8yk2lLISkLLKpbymPM0yejn6OwYvdVN28/GvoXO1TG6b2YPJ9qnQw9kc0fBEDR5WJ/5UijurAudE1T+xl3Chv7dzczdf9RZm1DY18Ex/JXaxPQ5YVNWXazSpdoqP/s/Oz9FbLS1epTrvf5Vo6aAdwZ9OBNqT+nOkdt88inXiY7wk4K+Nxzpn2lxXNRvnucS4e7OCYRJobUKbAA9JjlTUffhLF/f5Y6g9L1nv5d+iEWbO/q7yrLRcMbOldGr9pU6huU+VZY22OsE0RL10yUTK65pmqXpNjTTIXx0aW451fvXSVec3sAHfZuYR5U4gKB3WXnK6+IfrucY4h4jvqWL8c1Z7O9mJUZhPM5GJX0Z/tLq6Ax7liB4rBxk6EklWBHQ2eOFuhtfttdFc5Y+eLrLxOvfyHpks2qF48sDSdk9rITEe/+aQqZxKqd6vA2fpWRdaXlUzsaclJeOp+vWiacCYFHLNW8GLX4qXMy39ouo3q5GqfhuFgZItL4/awSHpw1GyccL6W/FRALs4yrxcndejeQKZ4vPi09E7VxyvTS9fSvVFd6JO8sp/qY4lvgq5Lt2TQ5pw3t37OOGMRYQHC+eSVLjtKLrpk2Paj1VbsoCLniSJyya8Z53yrvZ4LdJmccNLHfOwTZBcy2JjokcVu8dsJvV1t0pDh5ffSjVVM9XejW2pLHiamgmdd2gi7wfjVNX2OyBrAF03xLXDOXudIH012eelEPxeQNTnPSMfkmNQ919jIdN71b5FO9q5t74SD7NnBvAXQXN1Fv9Q39hnnmLNqrGebu9f1SwNAAox57LjyBR3vedI5Z8xU3ctU9ccR5fftmBz+ofzcIr2qDshUhLc10VpZq4MlfV7+Q9OhXDLZBXisdNK+ZOT2g8WKV0JzGDLBeWRVYnez0mEL+eCmypkwcaKWdsjGfGfZRqpnT/sM6qyLsriMb0gu7bRR+nZ0j0WQOYPIxMobnCulXyhtBevfyFqUKRpkIqfwVUKtbKr/NOsjea629dJV9pSseFkcwYF5czJM0a3QzYN4/eC4kQywY3OB+Sy1j1IXXaV/buzS1wjd41zysYf2Y7e3DpQTxv1sW9rxSR1EPtUPOyRtSStMrN/VMqxsqi9uvPxeulp5zi/1d/CYCvU5VqtLTWtlpX6yc3H8SkDa1KFOyg/maJXRFwm1TC+de57Zqkhz19KYXBwbspl+ubj+ic49D4ONImYaxuQtWFnCSHRuzE0AqfiadVe56Qb3qbB6nCNQOlaPdfEO+owZqPLZdf3SCJVOOa4YZZ2TBt0BQGVUHEX89uTSTgx07JCVzqzrgwXJ5ZvyXU6MKRffz4os1DzdH2xCNPk5tVeCo+JyKH78QDGm8/Ifmo42/lL4pFcPlVEMNNqcch58WOSgLefhVMbgZHFNIV/Tr2Z3sDL5KRIvdlO2efkPSfdExnjbp9id24L2mnIkUjuKbmknFv3vi+DbzBVZ8dfyvTK9+BRtjvp4ZXrpat02B5Z+LXvAlvqX+U75JTqXbslhXE3NyYxF5m2caxcdduYwpd/uW13WjHOwYf4fryPPVI7NJtt03Zd0CpPUv2Uk9+eCl99Ll3Tltp0bv1PyzFbDe4quVT/v+DUd45Q1I705GN8YXbfoeuaZjfCZHWsZv8UxJLvAeHF+FQ0O9yKd9LJOEry4d2O+VPcZ3ckwfVjfWDPOkXH0sV47tDgY9aKzESBUhEgjEZ4rWkNwnw5H2StFXifwpMPZqI8VC53KUlAZTwuDSUtl9rT2TveafFvu2885Ht3jnN/Ov37ccs/nxMtkC/9d/QMGjp/VuTYGB7o10dc05L38B6UTHgw84iConL8hZHIqi7iu6UcsohboWzXvS13z8JAelIzoHqVe7KZM9vIfkq6nfWq7WSAILWeUchyP1uJjfNautOfOmBfvuO17ZHrxQaYFs2uqPl6ZXrqkV3VnLmQhG2PAglrOojrpXLoli10W4tQmAmMwtZ2HLlVk+8HrUmwYB+Yo9FmbrxnnO/1B8sCOeFebCON6rLl2tcmMYC+/l85ULfV3b1v26HXNCWpX+j5z/Sf0PwxWSp+kD5S29tLBr+CeZyR3cUyKxjuGcO6s3ydD+BB/a/3z0iHCi7sLcwQSZNdi3UXm7RtrxjlmHH2sX6BVlaeT2aRCkQWbnHFCcWzrXTboOaP1tSKOyTul0LE7+lQp90vQNTo431UaW3kcXPjqb8Gb81x464yThzOwLaexFjWZFy8LEIPx4CHLfq+0yFcebP6lmA7Ho5QyxRvF9GRGGUHX2ObhPyhdUt7+wHZiHahH/cqAhT497SllAGN/fb/mPXm+A+NzaKNW+9QYW9td14VVnjlgMCaFD30XPr74YHw8BDJBlqBrzqQT6ocdrl0yxe/qwwiswmx9vDK9dOgVLXMj/Zr+kHCwVGUFow66nnqDLfOd1VuXySbmVn5VxHZZvHQb8TCf4wiXOVR55A/mKF13jfMsA6yuFFPIZWAHTiyc9zLINlebUB/FO5+3K5Cs3a+rspKVLa629NYvC3aNX9HS1uMHy1Zbe+lQ75pnVB/XmESggntsbMkHn+BvbTC4Mbpo0nXg7sV846276Fx9Q/XoGufUW7Kp79HH+qObmxuU40mzaDFomWxeK1LObieL0seKOJ7lVb7y0NO5AGWjtLzqVx4nDCe1dl7T08CojKcbzquWiUz5G5VxdMEmYl3eBpW7eESHDavPwoqfyZwFodTh1or9cpJJYyP/WpGd6S8VX6i84KDrja55gGBhGTh/uvbyH5QOmyzIBtqTDkvfIdAPXqs89RGlPNBgN3X4SpG6YQ99ZlBPlR08SAdOFztXpc/2KBGfF7t72Uayf7Z9DAvR0YaMlW+UT2PZ7lmqctqYSd8CPLzBGIyNLIuJl8A88l6ReeCagjp0yHS1g8nONizVxyVTsrx09AEwaQXeuqTdKKUuOoR4dWdaxtpzRfC2gJM7bh8XnQkQP3MUGBBoz9YchUzXOJc8xiQODed77YGe/FHmBOnZK3jbRHRHmROojHQtjt9M52lLV3/P8rxzAnoJ1o92+iU3VQ8XXaalzrPzjOS5x1qW2Ts2vPPrIp1sdeEuOi/mvXX39A33OM943slYFwaz67o5tIDOb5GV1wAYNRdEi2OZXiUoT2U5q4XzlQaY0sEPY+uan5Epv4KQ6QB+cDZP5chtLq49PKKlkXCwmou07kU4AgK5zZ4q5WHpaEH6Zjv+0QwJRYHAA0DAM85F073OPADoooqBwNkg4BnnVOauxrrkzq7rFxlJPP/B03wun0ugh4+Aw4Izy5MGk9Y/KLSgcpzca7vOqT2xjYrTLsPVuHAFD7vM7AxGOC0CPOz09q1DWEx/G/e5Q8gNGYFAILCLgGecr1lndjVFSSAQCJwKAc84x7a7Guuz6/pldkJxOAfn5JbQEh9nZ3ntyOvyz5UaC+VjRwL59eswo22lOLOPWzdmylo82IDeCKdFgPOzR98ll85TONGnRTq0BwKnQ2B2nGs8MqczH6fjbKczMzQHAoHAHgjMjnPk3uVYX1rXL6X/Sa7cHzl1JxKezrEoLednJ5iZzHAw6zC+tnvQNs/PqryHx+tAm95I7wYBzghHCAQCgQ8bgaVxbutMPGh+2P0gavdhI7A0zqn9ycb6hZTjyPLrA1POIgZOBvGxBf12kmB7AwcVR7UE8VGGztYualNeJw9ypxzjYkdk7hYBtdmqL2XdrVUhPRAIBA6JgGOc77XOHNLWkBUIBALrEHCMcwSfbKxz5ACncp+nZr5MNvuFH93nm74tx/WFdHPWIjmwoiHPv9wkR1QpTvBzpfaNRl2mv+ad5IEgB/TxpbMIgUAgEAgEAidEQHP4vuvMCa0P1YFAIOBF4JRj/cJr5BSdjJ91Zis++y3bUiRedu/4MwXO4vLzYPw0EL95ZgGn9Fvd45ttKTh4jJQvhP1mF5EGAoFAIBAIBAKBQCAQCHyYCKSf7TpG1eSIsqvK0YZ6t9WlGodW0e2cipadXX7vLn7lwIVwEAUCgUAgEAgEAoFAIHC+COy9Q+utupxLjjXwLyrsurqD6HGEe49E8Bu0gz8jcCsMwkAgEAgEAoFAIBAIBAKBs0LgaA4tqMg5xcnkLw7ZQfUG/u2J81euIFr+8pHd2fhCmAuxIAoEAoFAIBAIBAKBQOC8ETiqQwtUcjS7jhyI3ntG11riF/H07ugab6SBQCAQCAQCgUAgEAgEAmeGwNHO0J4ZLmFuIPBgEdADIW9QOOrD8SD+AbD5M3o1QGt4av7IBwKBQCAQCAQC+yBw9B3afYxt8Woh5XdwB4HFVZHFOEIgEAh0IKBxgyP7XBEnlrcj/Czfr0onwxqeSWFxIxAIBAKBQCAQWIHAYIdWCxNfpiLwbxCfKfJzWu6zqB5+0eBo2rED/lHiPdcq39kFcsr7W/zsKBm/nc/9QvyDs7c9uiVvU+nnkvBKZaaH+9TFvnyGXq7BrHnkQeU/6b4F/q73G5UNbLSbKre2sKKBbgpFg06cDwtcc364qd+IdJ+HAM4ym+12a6My07vYB0S7WB/R9GKEbThQO+1XjMwZyXbT1rzi+1ax9yhLLeIgedngxrql0MMvmp7xlvqTeGx8JrW65if1+MtDflZvEFTWzTMQoAvJuBftMbYrrgOBQCAQCATOB4Hi0GpReSOzXyhNP4+llIWKMn5qa9Gp9fCLhsX1J6Xl57SUZ1FnwURPccSUd9kjOv48AecQe7ET+6nHtdISdN2jG1qcKhztZJNSrj9XiqO/UYo+nNfiFCrP7+VC97Xy5WfGMi31oe7pn7OU4oz9TxHHreCr/KJu8aQgWuQV/RTqGv04v0V/Iq4+dA/M+LOLgYOiay/m1jdm6yN5Lowy3UvZxMMNDzlg85HKB22osk0PLfStIBm028Bpa9HdZZn0u7CessHDLxp3n0eP6BmLfKnyE+UL9srTjjw40lfLA52uN7ru5oGvDshQPGl71PZEPhAIBAKBQOD8ELjAZC0mLGK8pi9OkPIsaFzXu3C63A0d/Cx+YweMhQxdOGIpdMiDHscM5+eR4meKOKFlMd5KTJ8u3ZnenMLiYKucRb04nsqDGTtL5U8fdG309a6pijc4a1eiLX8DqzyOwR+KY3w9ujfiRz9O0TjwxxRj/YVGfDw87IQsz9sHvPVxYSTd14o8BNA3Xu0YVxX00FZs4+y7ccExr1WHezXeqrrzoMOXKgfjJ19TxnGEcVjDM5Zx0vYYGxPXgUAgEAgEAueHQHJoZTY7dYOdl1yV10qfakHDmZsLXn4WxD8b8nAE63OvXnlzNo3vuXTLNhxUdggHr6RVzg5y2VnWffBikS+Lv+6XvMrrgMzaGbZ7yCj4duiGn53i2h7KZoPkU6+BzRVDD+au+kh2D0aVKXeXFQb0g1Zfvzulu5J7sN7l9o9XV583BcLmZ8XBA6fdU/pe8cvqOmXX8NQy7kl71CZFPhAIBAKBQOAMEbjMNrPwDRy4XG5OGPfL7m2+VydefhxXXttPOX6Ps1CvvNqGpbxXNws6O4ZTNiY9uo+8j2qlKrPd2rLrqjKrEw7BONjO1BPdQJ5LdxbCw8Z3kv+7UnY3zV52oov+TGvJM9Gxg91yWlyYi9ddH9G6MDLjjpTS/8pO+ZF0jtW4sB4zVddefm+fL6KFDbvHXyhaf3qtMsb+p4o2Hyh7G9bw3HKn+eDU7VGZE9lAIBAIBAKBc0TgUouROShz9l9N3ezhFy07U63AzuFG9zk+0G2PeNIrXIn4WJGFlzO0g104XS/qFh8B5/Iv0WPTM0WcTnZDZ79sJXqcjHSsQfnycKA8zrFupXO+pHXAXgI2E9y6JfM3onhwov9WnqMb2Pl7Llf2NqiMowZNR1f33JiLtqc+twYoJ94mRgOiO7iQ3rot2WlHC1hNHU/h/p0E6XZj3TKgh1+03j6fVImehyP6fnngUZ43CIwvwo5Du5Ln3rTHtlrxGQgEAoFAIHDuCFyoAle5EtczlZlbhPfi14LI4oZDZ18K6ZWHbZz7+1ERGcQ3yuM8zYaGbuitrvxcEQ4PclngcWhxHgdBZez44SxCgxPNudhxsB2ucTl1J5hOS726cVjMecaZps4DR17XG9kHvjiiOw4J9xV6MffWJwmXXg9GifbQH9KNE88vOqS2VB6Hnx1B8KKfgM0xQy/WY9v24ld96XP1eEvyVY4zu1FanNl8bW8OuHzDh4WVPPetPaw6kQYCgUAgEAicMQIXTtttJ9FJvkM2x8+XoNht7HntWOSJjx23a9OoPE4bi3BzN9LocjrQLV5zKHHAzFE0ll+UeVnRpHJds6uM04tz+UoRJ2ns+PJFrY3Ki5OtPI6F2c2u2Brd6EEGu43UGUeFM8pj/Th04/qItCsUzMW1WJ9asnR7MKpZDpKXXtp3ozQ5aUrBH5wos36CY3vfQo31Gtvm+Ad9HuHCgh1YsLGHSorrYH2TMZDCSp5zbQ+rdqSBQCAQCAQC9xSBS9nVOttp5tpukJ31tPI6Xc2vRRGnE2eufjW6Wl5lFM7KU8n9VLG5Kzmh20S0eNidYuF/opicIiO2VDJxzHEw2c0tPztFmeInKufniXBkcThe54jzWeur87qVwo5uycEWfkbJdtRw7JGF04DjzT88oRe6Jee+C/Ms11sfqb8N4m1idEtxmFyuN04adloAo9ppo62gmQyix5nj59XMqZukrW5wpvltdV1nu7CuGXN+Nb9sao03xIIJfWXKZh6UrN2gJ3TxSLY5zXu1x1Z1fAYCgUAgEAgEAkMELrXQsJBR2lqwrazlZCVJa/nzAneldPBN/R55ouU1KTK+SMbsfpj9gzuiZ3Gd0309YBhesLhvJAPHlHTsBHDkACeJyKv5FESHTHM+rcx2B3HqrR0WdW8lplfmtXOALTgd7NbyU0o49Nj2WOlk+yFL9013Cy8rG8iAR6yT9clyuzCC54ABh4ujKHN4Yt+gXmP9mX+qf43JF6/XYF0LXcsvvqk+T/smh7XWY3nxmcNfHopU1s0jeQdpD7Mr0kAgEAgEAoFAoEbgMl/Y6+r6Hvmr6v74Xn3dxa8FkZ1EfjO27Mwqb44iDoZXHrulrR2rZLdkjp3NjcrW6rb6mgOUzhNKXtmJNYKOFIcq7aRmnql6m8ikWzpxKHBUd5w1leEc40iDAZh+qev0qld5C+hl95py6HE2pnQnLPN9JbNhXJ9DYDSrcOYmdU/6oVEdwQxHP4V8jbNW+mC+dYxkX6y7+FXXyT5fVZY3Bq1A3+DICDotWJ/o4bnP7WH1ijQQCAQCgUDgTBEwhxbHxnYL66qwM8VituM41UTKu/klC6cHJ4uFsg4sunbO0yuP380cy0Emjkq9AFO2cepmJwr94wAW7GSaXDCpnVGjx8kmGB16qdtLxfIPTCrDwcLOevfPpVu82EGcOlKBbGyzhwNd3gaV869P3K+dOS/mG/F56+PG6Na6g+XS7nQlbdwnaA/6T9lFr2jvOuvGesIQN7/qNzvedD/1pZae3M7gxq5/CSpPbxRKQZWZ4hHJfW6PqgaRDQQCgUAgEDhHBC4wWosQjuT7vBileiiPU/QvxfQFIAopU7xRLDtfudzLzy4NizFy+NvUElXGF5dwgDZKXfJEmvjhsSBe+yes2llDplc3Dg7OXnHwld/BQjQ40umb4ZVuHM65TNIAAAHkSURBVD1oS13yPXSPd5LBAbqyi6y8VzdiqR9nddFXgq55tcw3+W0nudyrMvCM+byYI8ZVH9H1YIRcwsfbpLwdyJfNZI6WflsfieAhyna4wZ7+Xt9vKriLQul1YS26o4w31ZF2GmAh3TjCOP2Dv2bWtYVennvbHlahSAOBQCAQCATOF4FHNzc3yXoWT2Vw4q4V3yl+qfhC5cXh0vVG17y2xeEbL4CL/JkXZ6gV2Akuu5XKL8pDiOiQx+JK4FUojuPO74uKDrtdukWHXLDABgJyW1iwe1U7zsjni19ld1bXKWR55E3m5O/aenQjSHQ4Hc8Va2d5Ti47wNiI3QQc6NeSk35hQqkLcxhFaw7/bH1E58JIdDiZBOiRSb/DAcU5xwEswUsrOvChn9L2nNVGJrJ5ECJ/siD9LqxFd6zxxsPYM8W6L+2Moxow2dbFI/p72x51vSIfCAQCgUAgcH4IFIf2/EwPiwMBHwJypJITr3TgGPu4g+rQCER7HBrRkBcIBAKBQCBwERAEAg8AAXYGd3bNH0C972sVoz3ua8uEXYFAIBAInCkC4dCeacOF2V0IlPOzXVxBfFcIRHvcFbIhNxAIBAKBB4pAOLQPtOEfWLU5Ex7h/iAQ7XF/2iIsCQQCgUDgg0Dg/5cbX+BvwiqcAAAAAElFTkSuQmCC\n",
387
      "text/latex": [
388
       "$\\displaystyle {{\\mu_{\\phi}}_{(0,0)}} \\leftarrow 0.0004 \\phi^{3} + 0.000473530700220886 \\phi \\rho^{2} - 0.00760399528440326 \\phi \\rho + 0.0205263968409311 \\phi - 0.02 {\\partial {\\partial \\phi}}$"
389
390
391
392
393
394
395
396
397
      ],
      "text/plain": [
       "                 3                           2                                \n",
       "μ_φ_C := 0.0004⋅φ  + 0.000473530700220886⋅φ⋅ρ  - 0.00760399528440326⋅φ⋅ρ + 0.0\n",
       "\n",
       "                                  \n",
       "205263968409311⋅φ - 0.02⋅D(D(phi))"
      ]
     },
Markus Holzer's avatar
Markus Holzer committed
398
     "execution_count": 12,
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "μ_ρ, μ_φ = chemical_potentials_from_free_energy(free_energy_transformed, \n",
    "                                                order_parameters=(ρ, φ))\n",
    "μ_phi_assignment = Assignment(μ_phi_field.center, μ_φ)\n",
    "μ_kernel = make_kernel([μ_phi_assignment])\n",
    "μ_phi_assignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pressure tensor and force computation\n",
    "\n",
    "For the pressure tensor a trick for enhancing numerical stability is used: the bulk component is not stored directly in the pressure tensor field, but the related quantity called `pbs` is stored in a separate field.\n",
    "\n",
    "$ pbs = \\sqrt{|ρ  c_s^2 - p_{bulk} |} $\n",
    "\n",
    "The force is then calculated as $ \\nabla \\cdot  P_{if} + 2  (\\nabla pbs) pbs$\n",
    "\n",
    "In the following kernel the pressure tensor field is filled with $P_{if}$ and the pbs field with above expression."
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
428
   "execution_count": 13,
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bulk part\n",
    "pressure_assignments = [\n",
    "    Assignment(pbs_field.center, \n",
    "               pressure_tensor_bulk_sqrt_term(free_energy_transformed, (ρ, φ), ρ)),\n",
    "]\n",
    "\n",
    "# Interface part\n",
    "P_if = pressure_tensor_from_free_energy(free_energy_transformed, (ρ, φ), \n",
    "                                        dim=dh.dim, include_bulk=False)\n",
    "index_map = symmetric_tensor_linearization(dh.dim)\n",
    "\n",
    "pressure_assignments += [\n",
    "    Assignment(pressure_tensor_field(index_1d), P_if[index_2d])\n",
    "    for index_2d, index_1d in index_map.items()\n",
    "]\n",
    "pressure_kernel = make_kernel(pressure_assignments)\n",
    "\n",
    "\n",
    "# Force kernel\n",
    "pressure_tensor_sym = sp.Matrix(dh.dim, dh.dim, \n",
    "                                lambda i, j: pressure_tensor_field(index_map[i, j] \n",
    "                                                                   if i < j else index_map[j, i]))\n",
    "force_term = force_from_pressure_tensor(pressure_tensor_sym, \n",
    "                                        functions=[ρ, φ], \n",
    "                                        pbs=pbs_field.center)\n",
    "force_assignments = [\n",
    "    Assignment(force_field(i), \n",
    "               force_term[i] + external_force[i] *  ρ_field.center / ρ_l)\n",
    "    for i in range(dh.dim)\n",
    "]\n",
    "force_kernel = make_kernel(force_assignments)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Lattice Boltzmann schemes for time evolution of ρ and φ\n",
    "\n",
    "- ρ is handled by a normal LB method (compressible, entropic equilibrium)\n",
    "- stream and collide are splitted into separate kernels\n",
    "- macroscopic values are computed after the stream, but inside the stream kernel\n",
    "- velocity field stores the velocity which was not corrected for the forces yet\n",
    "- the φ collision kernel corrects the velocity itself, because then the updated forces are used for the correction. When u is computed, the updated forces are not computed yet\n",
    "- when ρ and φ are updated, they are clipped to a valid region, this clipping should be only necessary during equilibration of the system\n",
    "- exact difference method is used to couple the force into the ρ-LBM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell handles the clipping of the order parameters:"
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
489
   "execution_count": 14,
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
   "metadata": {},
   "outputs": [],
   "source": [
    "if clipping:\n",
    "    def clip(ac, symbol, min_value, max_value):\n",
    "        \"\"\"Function to clip the value of a symbol which is on one of lhs of the assignments \n",
    "        in an assignment collection\"\"\"\n",
    "        assert symbol in ac.bound_symbols\n",
    "        for i in range(len(ac.subexpressions)):\n",
    "            a = ac.subexpressions[i]\n",
    "            if a.lhs == symbol:\n",
    "                new_assignment = Assignment(symbol, sp.Piecewise((max_value, a.rhs > max_value), \n",
    "                                                                 (min_value, a.rhs < min_value), \n",
    "                                                                 (a.rhs, True)))\n",
    "                ac.subexpressions[i] = new_assignment\n",
    "                break\n",
    "\n",
    "    # TODO: how can this 'densgin' be derived automatically?            \n",
    "    tred = reduced_temperature\n",
    "    densgin = -67.098 \\\n",
    "              + 549.69 * tred \\\n",
    "              - 1850.6 * tred * tred \\\n",
    "              + 3281 * tred * tred * tred \\\n",
    "              - 3237.3 * tred * tred * tred * tred \\\n",
    "              + 1687.6 * tred * tred * tred * tred * tred \\\n",
    "              - 361.51 * tred * tred * tred * tred * tred * tred\n",
    "    ρ_clip_min, ρ_clip_max = densgin * 0.5, 1.2 * ρ_l \n",
    "    φ_clip_min, φ_clip_max = -χ * 1.5, χ * 1.5          "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the collide and stream kernels for the ρ lattice Boltzmann are created"
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
529
   "execution_count": 15,
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
   "metadata": {},
   "outputs": [],
   "source": [
    "force_model = EDM(force_field.center_vector)\n",
    "\n",
    "ρ_lbm_params = {\n",
    "    'stencil': stencil,\n",
    "    'method' : 'trt-kbc-n2',\n",
    "    'compressible': True,\n",
    "    'relaxation_rate': ρ_relaxation_rate,\n",
    "    'optimization': {\n",
    "        'symbolic_field': pdf_src_rho,\n",
    "        'symbolic_temporary_field': pdf_dst_rho,\n",
    "        'openmp': threads, \n",
    "        'target': target\n",
    "    }\n",
    "}\n",
    "\n",
    "# Standard collision step, that does not compute ρ and u from pdfs, but reads\n",
    "# them from fields - this is necessary because ρ may have been clipped before\n",
    "# the velocity field is not force corrected, which is the correct for the EDM model\n",
    "# but might be wrong for other force models\n",
    "ρ_collide = create_lb_function(kernel_type='collide_only',\n",
    "                               density_input=ρ_field,\n",
    "                               force_model=force_model,\n",
    "                               velocity_input=vel_field,\n",
    "                               **ρ_lbm_params)\n",
    "\n",
    "\n",
    "# First the assignments are created, then the density is clipped\n",
    "# then a kernel is created from the clipped assignments\n",
    "ρ_stream_ur = create_lb_update_rule(kernel_type='stream_pull_only',\n",
    "                                    force_model=None, #save uncorrected velocity\n",
    "                                    output={'density': ρ_field, \n",
    "                                            'velocity': vel_field},\n",
    "                                    **ρ_lbm_params)\n",
    "\n",
    "if clipping:\n",
    "    clip(ρ_stream_ur, sp.Symbol(\"rho\"), ρ_clip_min, ρ_clip_max)\n",
    "ρ_stream = create_lb_function(update_rule=ρ_stream_ur, **ρ_lbm_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The φ lattice Boltzmann solve the Cahn-Hilliard equation. "
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
581
   "execution_count": 16,
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
   "metadata": {},
   "outputs": [],
   "source": [
    "φ_lb_method = cahn_hilliard_lb_method(stencil=stencil, \n",
    "                                      mu=μ_phi_field.center, \n",
    "                                      relaxation_rate=φ_relaxation_rate, \n",
    "                                      gamma=1)\n",
    "φ_lbm_params = {\n",
    "    'lb_method' : φ_lb_method,\n",
    "    'compressible': True,\n",
    "    'optimization': {\n",
    "        'symbolic_field': pdf_src_phi,\n",
    "        'symbolic_temporary_field': pdf_dst_phi,\n",
    "        'openmp': threads, \n",
    "        'target': target\n",
    "    }\n",
    "}\n",
    "\n",
    "ch_method = cahn_hilliard_lb_method(stencil=stencil, \n",
    "                                    mu=μ_phi_field.center, \n",
    "                                    relaxation_rate=φ_relaxation_rate, \n",
    "                                    gamma=1)\n",
    "\n",
    "corrected_vel = vel_field.center_vector + sp.Matrix(force_model.macroscopic_velocity_shift(ρ_field.center))\n",
    "φ_collide = create_lb_function(kernel_type='collide_only',\n",
    "                               density_input=φ_field,\n",
    "                               velocity_input=corrected_vel,\n",
    "                               **φ_lbm_params)\n",
    "\n",
    "φ_stream_ur = create_lb_update_rule(kernel_type='stream_pull_only',\n",
    "                                    output={'density': φ_field},\n",
    "                                    **φ_lbm_params)\n",
    "if clipping:\n",
    "    clip(φ_stream_ur, sp.Symbol(\"rho\"), φ_clip_min, φ_clip_max)\n",
    "φ_stream = create_lb_function(update_rule=φ_stream_ur, **φ_lbm_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Time loop\n",
    "\n",
    "Now we can put all kernels together into a time loop function"
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
630
   "execution_count": 17,
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
   "metadata": {},
   "outputs": [],
   "source": [
    "op_sync = dh.synchronization_function([ρ_field.name, φ_field.name])\n",
    "p_sync = dh.synchronization_function([pbs_field.name, pressure_tensor_field.name])\n",
    "pdf_sync = dh.synchronization_function([pdf_src_phi.name, pdf_src_rho.name])\n",
    "\n",
    "def time_loop(steps):\n",
    "    for t in range(steps):\n",
    "        op_sync()\n",
    "        dh.run_kernel(μ_kernel)\n",
    "        dh.run_kernel(pressure_kernel)\n",
    "        \n",
    "        p_sync()\n",
    "        dh.run_kernel(force_kernel)\n",
    "        \n",
    "        dh.run_kernel(ρ_collide)\n",
    "        dh.run_kernel(φ_collide)\n",
    "        \n",
    "        pdf_sync()\n",
    "        dh.run_kernel(ρ_stream)\n",
    "        dh.run_kernel(φ_stream)\n",
    "        dh.swap(pdf_dst_phi.name, pdf_src_phi.name)\n",
    "        dh.swap(pdf_dst_rho.name, pdf_src_rho.name)\n",
    "    return dh.cpu_arrays[φ_field.name][1:-1, 1:-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3b: Compiling getter & setter kernels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The setter kernel computes ρ, φ from C and sets the pdfs to equilibrium using the values in the order parameter and velocity fields."
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
674
   "execution_count": 18,
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
   "metadata": {},
   "outputs": [],
   "source": [
    "init_assignments = [ \n",
    "    Assignment( φ_field.center, transform_backward_substitutions[φ].subs({ c_l1: c_field(0), c_l2: c_field(1)} )),\n",
    "    Assignment( ρ_field.center, transform_backward_substitutions[ρ].subs({ c_l1: c_field(0), c_l2: c_field(1)} )),\n",
    "]\n",
    "\n",
    "init_rho = pdf_initialization_assignments(ρ_collide.method, \n",
    "                                          density=ρ_field.center, \n",
    "                                          velocity=vel_field.center_vector, \n",
    "                                          pdfs=pdf_src_rho.center_vector)\n",
    "init_rho = init_rho.new_without_subexpressions()\n",
    "init_assignments += init_rho.all_assignments\n",
    "\n",
    "init_phi = pdf_initialization_assignments(φ_collide.method,\n",
    "                                          density=φ_field.center,\n",
    "                                          velocity=(0,0,0), \n",
    "                                          pdfs=pdf_src_phi.center_vector)\n",
    "init_phi = init_phi.new_without_subexpressions().new_with_substitutions({μ_phi_field.center: 0})\n",
    "init_assignments += init_phi.all_assignments\n",
    "\n",
    "init_pdfs_assignments = init_rho.all_assignments + init_phi.all_assignments\n",
    "\n",
    "init_pdfs_kernel = create_kernel(init_pdfs_assignments).compile()\n",
    "init_kernel = create_kernel(init_assignments).compile()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4: Geometry initialization & plotting"
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
712
   "execution_count": 19,
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_status():\n",
    "    plt.figure(figsize=(25, 6))\n",
    "    plt.subplot(1, 4, 1)\n",
    "    ρ_arr = dh.cpu_arrays[ρ_field.name][1:-1, 1:-1]\n",
    "    plt.scalar_field(ρ_arr)\n",
    "    plt.title(\"ρ ({:.2f}, {:.2f})\".format(np.min(ρ_arr), np.max(ρ_arr)))\n",
    "    plt.colorbar()\n",
    "    \n",
    "    plt.subplot(1, 4, 2)\n",
    "    φ_arr = dh.cpu_arrays[φ_field.name][1:-1, 1:-1]\n",
    "    plt.scalar_field(φ_arr)\n",
    "    plt.title(\"φ ({:.2f}, {:.2f})\".format(np.min(φ_arr), np.max(φ_arr)))\n",
    "    plt.colorbar()\n",
    "\n",
    "    plt.subplot(1, 4, 3)\n",
    "    f_arr = dh.cpu_arrays[force_field.name][1:-1, 1:-1]\n",
    "    plt.vector_field_magnitude(f_arr)\n",
    "    plt.title(\"F ({:.2f}, {:.2f})\".format(np.min(f_arr), np.max(f_arr)))\n",
    "    plt.colorbar()\n",
    "    \n",
    "    plt.subplot(1, 4, 4)\n",
    "    μ_arr = dh.cpu_arrays[μ_phi_field.name][1:-1, 1:-1]\n",
    "    plt.scalar_field(μ_arr)\n",
    "    plt.title(\"μ_φ ({:.2f}, {:.2f})\".format(np.min(μ_arr), np.max(μ_arr)))\n",
    "    plt.colorbar()\n",
    "\n",
    "def init_drop():\n",
    "    radius = dh.shape[0] // 5\n",
    "    mid1 = [dh.shape[0] // 2, dh.shape[1] // 2]\n",
    "\n",
    "    for block in dh.iterate(ghost_layers=True):\n",
    "        x, y = block.midpoint_arrays\n",
    "        mask1 = (x - mid1[0]) ** 2 + (y - mid1[1])**2 < radius ** 2\n",
    "\n",
    "        block[force_field.name].fill(0)\n",
    "        block[vel_field.name].fill(0)\n",
    "        block[μ_phi_field.name].fill(0)\n",
    "\n",
    "        c_arr = block[c_field.name]\n",
    "        c_arr[:, :].fill(0.0)\n",
    "        c_arr[mask1, 0] = 1.0\n",
    "\n",
    "        gaussian_filter(c_arr[..., 0], sigma=3, output=c_arr[..., 0])\n",
    "        gaussian_filter(c_arr[..., 1], sigma=3, output=c_arr[..., 1])\n",
    "\n",
    "    dh.run_kernel(init_kernel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 5: Putting it all together"
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
773
   "execution_count": 20,
774
775
776
777
   "metadata": {},
   "outputs": [
    {
     "data": {
Markus Holzer's avatar
Markus Holzer committed
778
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABakAAAF1CAYAAADiCWaNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9WklEQVR4nO3dfbisd13f+8+XRBAUSSAxpkk4iRLb4hNqirYerYcgRKWEc4oYSjF6haL1sbXWYj0FRT0HWi3oKbWNEAmKPBisxJqahgCt9VxGglAUkLKJYhKDhCQgRQGT/T1/zB3OsFk7e+01M2vu39qv13XNtWfuue+Z3+yV671XfnM/VHcHAAAAAAC24X7bHgAAAAAAACcuk9QAAAAAAGyNSWoAAAAAALbGJDUAAAAAAFtjkhoAAAAAgK0xSQ0AAAAAwNaYpCZJUlWnV9UfVNUDtz2WTamqn6qqf7jtcQCsqqpeUVVP2vY4NqWqvriq/t9tjwNgFVX1gKp6R1Wdue2xbEpVfU9VPX/b4wDYFHMlsH+qu7c9Bmagqn4qye3d/bxtj2VTpv9B+J0kn9fdH9/2eAD2oqq+OMkrk3xBH+B/xKvqmiQ/292/tu2xANyrqv4oyRlJ7lla/Pnd/Sc7rPs9WbT6O/ZpePuuqj49yaEkX9bd79/2eADWzVwJ7B97UpOqekCSS5P84rbHskndfVuSP0jyxG2PBWAF357k5Qd5gnry8iw+K8Dc/J3u/syl26dMUE++I8kv7OfA9lt3fzTJf0ryLdseC8C6mSuB/WWS+gCrqn9SVbdX1V1VdcW0p8NOviLJB7v7lqVt31hVP1ZVv1VVH66q/1xVpy09/8tV9b6q+lBV/deq+oKl515aVS+qql+ftr2hqj5vl2P+1qr6b1X1k9O4/7Cqvn7p+W+rqndOr3tTVX370nNfW1W3TJ/7/VV1W1V92xFv8cYk37ibsQDsh6r6W1X136vqa6Zmv7mqHnkfm3x9kv+ytP2mu3lfY++q+o6qendVfXBqf03PfV5Vvb6q7qiqD1TVy6vqlKVt/6iqfqCq3jb9W/KqI/6demOSC6f/OQAYSlU9PMnnJrlhadl9/o5cVT9dVTdX1Z9N/xZ89dJzP1JVr66ql03bvr2qLtjlWO6z9VX1jVX1lul9b66qH1l67typ9ZdW1R9PPf/hI97ijfH7NbDPpjY9YunxS6vqx3ex3XOq6l9V1Qun318vr6pPO8rq5kpgH5mkPqCq6v9I8g+SfFmSz0vyV5P8n0dZ/YuSvGuH5X8vybcl+ewk90/yA0vP/ack50/P/W4We7wtuyTJjyY5NYtDAH/iOIb/FdN4TkvyL5O85N5JjyTvT/KEJJ81je0FVfVlS9t+TpKHJDkryWVJXlRVpy49/84kX3IcYwHYmKp6UJJfTvKiJE/LoqWvTPLype4tr/8ZSc7LpzZ7k908lick+RtJvjjJU5I8/t7hJvm/k/yVJH89yTlJfuSIbZ+S5KLpM31xkm+994nuvjXJX2bx7xfAaL4oyU3dffcRy+/rd+Q3JXlUkocm+aUkv3zEl3dPzOLfiFOSXJ3k3xzHeO6r9R/JYk/oU7KYoPiH9anXPfhfs+jxhUmeXVV/fek5v18DQ6iqi5N8Z5L/J4ujXb4hizmT7z7KJuZKYB+ZpD64npHkhd19c3ffmUX4nnqUdU9J8uEdlv98d/+P7v6LJK/O4pfmJEl3X9HdH+7uj2Ux6fAlVfWQpW3/Q3f/zvSL+cuXt92F93b3z3X3PUmuTHJmFuf+S3f/ene/pxf+S5L/nOSrl7b9yyTP7e6/7O5rkvzPfPIEx4enzwswB4/OYjL355J8RpIPZvFL8xclefgO658y/XlkszfZzWN5Xnd/sLv/OMkbMvW+uw9193Xd/bHuvj3Jv07yt4/Y9me6+0+mf6d+LZ/6b4VmA3P0q9Pedx+sql89yjqnZOffr4/6O3J3/2J339Hdd3f3TyV5QD65x/+tu6+ZWv8LOb7JhKO2vrvf2N2/192Hu/ttSV6RT+31j3b3X3T3f0/y34947w9nMfEBMHf/e5KfT/InWTT2T5L8uxz9NBenxFwJ7BuT1AfXeUluXnr83iz2ZtvJXUkevMPy9y3d//Mkn5kkVXVSVT2vqt5TVX+W5I+mdU471ra79Iltu/vPp7v3vvfXV9VvV9WdVfXBLL75XH7fO47YY+XI935wFpNAAHNwRpJbl88vPZ3f864s9nY40genP49s9ia7eSxH+7fijKp6ZVXdOv1b8YtHvO9Rt12i2cAcPam7T5luTzrKOsf1+3WS1OIUSO+cDhH/YBYTv/f1+/WnV9XJuxzzUVtfVV9RVW+oxSmnPpTF3oXH0+sHJ/nQLscBsE1n5JPnSZLktuz8e3dirgT2lUnqg+uBWRxafa+HZ/Et4U7eluTzj+O1/16Si5M8Notfns+dln/KoenrVIvzkr4myU8mOaO7T0lyzXG+71/PYu8PgDm4NcnZy6f2qKoHZnH43y1HrtzdH0nynuyy2Wvq5l79X0k6yRd192cl+fvH875VdVYWh0/udIglwNy9Lcl5u51ErsX5p38wi9MgnTr1+kPZn17/UhanDzmnux+SxV6Ffr8GRvDApfu7OaLj1nzq0YpnZYffuyfmSmAfmaQ+2L6rqs6uqocm+eEkrzrKer+T5JRpQmA3HpzkY0nuSPKgLCYidm260MCPHM82k/tncUjO7Ununi4S8LjjfI2/ncU5ogDm4IYszgX6XVn8EnlSFueo+63pnMw7uSafehj20azUzekCLX+02/WP8OAsDiP80PTvyz89zu3/dpLXT4dKAgxlusjWoSxO67QbD05ydxa9Prmqnp3FeUV3ZboY10uPd5xL731nd3+0qh6dxSTL8fD7NbAt3zrtvfyoLCaGH3wfF0FMFqfX+NYk915w8cwk35fFl3U7MVcC+8gk9cH2S1mch+imLPa82/FKt9398SQvzWIvt914WRanD7k1yTuS/PZxjuucJL91nNukuz+c5HuzOOfTXVn8An31brevqjOTPDLJrx7vewNsQnf/ZRZ7Wzw9i73n/mkWF295+n1sdnmSp+10YcUdXn+lbmaPvZ78aBaf5UNJfj3Jrxzn9k/LYm8+gFH9+9x3z5ddm+Q3kvyPLH7P/mg+9ZD0+7JKr78zyXOr6sNJnp3Fvxm7Ml3Y8RuyODcqwH57UBan6/i5LPr1LUkec7SVu/sNSX4syXXTotdOtx0bZq4E9lctnQaTA2Ta8+0Z3f26Xa5/epLfTPKl08n/NzWus5O8urv/1qbe4z7e+6eSvKe7/+1+vzfAsVTVLyY51N0/sot1fymLlv7qhsf0n5N8X3e/c5Pvs8P7fnGSf9/df3M/3xdgnabDr9+S5MLuvm2D73P/LA7R/uLpy899U1Xfk8VpQn5wP98XoKo6yfndfWgP256cxYUEz+vuPzrGuuZKYJ+YpD6gjneSGoDtOp5JagAAOJHt1yQ1sH+c7gM4YVTVRVX1rqo6VFXP2uH5r6mq362qu6vqyUc8d2lVvXu6Xbq0/Mur6vem1/yZ3ZyCAYD7ptcAY9Br5qaq/nlV/c8dbs63zAlthF7bkxo4IVTVSVmc5/Hrsrh685uSPLW737G0zrlZXKToB5Jc3d1XTcsfmuTGJBck6SRvTvLl3X1XVf1OFuf/uiGLC9r9THf7BQhgj/QaYAx6DTCGUXptT2rgRPHoLE6lcNN0AYxXZnHBuk/o7j/q7rclOXzEto9Pcl1339ndd2VxoY2LpgtMfFZ3/3YvvvF7WZInbfqDABxweg0wBr0GGMMQvTZJDZwozkpy89LjW6Zlq2x71nR/L68JwM70GmAMeg0whiF6ffIqG6/Laaed1ueee+62hwGs6M1vfvMHuvv0vWz7+P/tM/qOO+/Z+3u/7WNvT/LRpUWXd/fle35BdqTXcDDo9cGn13Aw6PXBp9cwvlVanej1vWYxSX3uuefmxhtv3PYwgBVV1Xv3uu0H7rwnN1x79p7f+9POfM9Hu/uC+1jl1iTnLD0+e1q2G7cm+dojtn3jtPzsI5bv9jWHpNdwMOi1XgNj0Gu9BuZvlVYnen0vp/sAZqJzTx/e820X3pTk/Ko6r6run+SSJFfvcnDXJnlcVZ1aVacmeVySa7v7tiR/VlVfOV3F9luSvPb4PzvASPQaYAx6DTAGvU5MUgMz0UkOp/d8O+brd9+d5LuzCOw7k7y6u99eVc+tqicmSVX9jaq6Jck3Jfn3VfX2ads7k/xYFmF/U5LnTsuS5DuTvDjJoSTvSeLK48CBptcAY9BrgDHo9cIsTvcBsB+6+5ok1xyx7NlL99+UTz5cZXm9K5JcscPyG5N84XpHCnBi02uAMeg1wBhG6LVJamA2DmdXh6kAsGV6DTAGvQYYg16bpAZmotO5p499mAoA26XXAGPQa4Ax6PWCSWpgNnZzLiUAtk+vAcag1wBj0GsXTgQAAAAAYIvsSQ3MQie5xzeHALOn1wBj0GuAMej1gklqYDYc3gIwBr0GGINeA4xBr01SAzPRiQsFAAxArwHGoNcAY9DrBZPUwGwc3vYAANgVvQYYg14DjEGvXTgRAAAAAIAtsic1MAuddqEAgAHoNcAY9BpgDHq9YJIamIdO7tFkgPnTa4Ax6DXAGPQ6iUlqYCY6zsEEMAK9BhiDXgOMQa8XTFIDM1G5J7XtQQBwTHoNMAa9BhiDXicunAgAAAAAwBbZkxqYhU5y2DmYAGZPrwHGoNcAY9DrBZPUwGw4vAVgDHoNMAa9BhiDXpukBmaiI8oAI9BrgDHoNcAY9HrBOakBAAAAANgae1IDs3G4fXMIMAK9BhiDXgOMQa9NUgMz4fAWgDHoNcAY9BpgDHq9YJIamIVO5R5nIAKYPb0GGINeA4xBrxdMUgOz4fAWgDHoNcAY9BpgDHrtwokAAAAAAGyRPamBWXAOJoAx6DXAGPQaYAx6vWCSGpiJyj3t4A6A+dNrgDHoNcAY9DoxSQ3MRCc57AxEALOn1wBj0GuAMej1gklqYDYc3gIwBr0GGINeA4xBr104EQAAAACALbInNTAL3c7BBDACvQYYg14DjEGvF0xSA7Nx2OEtAEPQa4Ax6DXAGPTaJDUwE53kHmcgApg9vQYYg14DjEGvF/wNAAAAAACwNcecpK6qK6rq/VX1+0vLHlpV11XVu6c/T52WV1X9TFUdqqq3VdWXbXLwwEGyOAfTXm+7eoeqi6rqXVOjnrXD8w+oqldNz99QVedOy59WVW9duh2uqkdNz71xes17n/vsNf6lHBe9BvaHXq9Kr4H9oder0mtgf+h1srs9qV+a5KIjlj0ryfXdfX6S66fHSfL1Sc6fbs9M8rOrDA44cXSSw7nfnm/HUlUnJXlRFp16ZJKnVtUjj1jtsiR3dfcjkrwgyfOTpLtf3t2P6u5HJXl6kj/s7rcubfe0e5/v7vev+FexipdGr4EN0+u1eGn0GtgwvV6Ll0avgQ3T64VjfpLu/q9J7jxi8cVJrpzuX5nkSUvLX9YLv53klKo6c5UBAieOe7r2fNuFRyc51N03dffHk7wyi2YtW27bVUkurKojX/yp07azo9fAftHr1eg1sF/0ejV6DewXvd77OanP6O7bpvvvS3LGdP+sJDcvrXfLtOxTVNUzq+rGqrrx9ttv3+MwgIOiU7kn99vzLclp9zZluj3ziLfYTZ8+sU53353kQ0kedsQ635zkFUcs+/np0JZ/sUPEt02vgbXS643Ra2Ct9Hpj9BpYK71eOHmVjZOku7uqeg/bXZ7k8iS54IILjnt7gCN8oLsv2OQbVNVXJPnz7v79pcVP6+5bq+rBSV6TxeEvL9vkOPZKr4GZ0Otj0GtgJvT6GPQamIkD0eu97kn9p/cetjL9ee85R25Ncs7SemdPywCO6XDfb8+3XdhNnz6xTlWdnOQhSe5Yev6SHPGtYXffOv354SS/lMVhNHOi18Da6fVG6DWwdnq9EXoNrJ1e732S+uokl073L03y2qXl3zJd1fYrk3xo6TAYgKPqZNXDW47lTUnOr6rzqur+WQT26iPWWW7bk5O8vrs7SarqfkmekqXzL1XVyVV12nT/05I8IcnvZ170Glgrvd4YvQbWSq83Rq+BtdLrhWOe7qOqXpHka7M4v8ktSZ6T5HlJXl1VlyV57zTQJLkmyTckOZTkz5N82yqDA04cnV2f8H9vr999d1V9d5Jrk5yU5IrufntVPTfJjd19dZKXJPmFqjqUxQVSLll6ia9JcnN337S07AFJrp2CfFKS1yX5uY19iGPQa2A/6PXq9BrYD3q9Or0G9oNeLxxzkrq7n3qUpy7cYd1O8l2rDAhgU7r7mix+eVxe9uyl+x9N8k1H2faNSb7yiGUfSfLlax/oHuk1cFDo9Setq9fAbOn1J62r18BsjdDrlS+cCLAuh/d8BiIA9pNeA4xBrwHGoNcmqYGZ6E7u2d0J/wHYIr0GGINeA4xBrxdMUgMzUTmczZ2DCYB10WuAMeg1wBj0OjFJDcxExzeHACPQa4Ax6DXAGPR6wd8AAAAAAABbY09qYDbu8b0ZwBD0GmAMeg0wBr02SQ3MRKdyuJ2DCWDu9BpgDHoNMAa9XjBJDcyGbw4BxqDXAGPQa4Ax6LVJamAmOslhFwoAmD29BhiDXgOMQa8X/A0AAAAAALA19qQGZqJyT5yDCWD+9BpgDHoNMAa9TkxSAzPh8BaAMeg1wBj0GmAMer1gkhqYDd8cAoxBrwHGoNcAY9Br56QGAAAAAGCL7EkNzEJ3ObwFYAB6DTAGvQYYg14vmKQGZuMeUQYYgl4DjEGvAcag1yapgZnoJIedgwlg9vQaYAx6DTAGvV4wSQ3MRPnmEGAIeg0wBr0GGINeJy6cCAAAAADAFtmTGpiFTnK4Hd4CMHd6DTAGvQYYg14vmKQGZuMeB3cADEGvAcag1wBj0GuT1MBMdMo3hwAD0GuAMeg1wBj0esEkNTAbh31zCDAEvQYYg14DjEGvXTgRAAAAAIAtsic1MAvdyT0ObwGYPb0GGINeA4xBrxdMUgOz4RxMAGPQa4Ax6DXAGPTaJDUwE4sLBTgDEcDc6TXAGPQaYAx6veBvAAAAAACArbEnNTAb98ThLQAj0GuAMeg1wBj02iQ1MBMd52ACGIFeA4xBrwHGoNcLJqmBmXAOJoAx6DXAGPQaYAx6nZikBmbksMNbAIag1wBj0GuAMei1CycCJ5Cquqiq3lVVh6rqWTs8/4CqetX0/A1Vde60/Nyq+ouqeut0+3dL23x5Vf3etM3PVJV/WQBWpNcAY9BrgDGM0Gt7UgOz0J3cs8FzMFXVSUlelOTrktyS5E1VdXV3v2NptcuS3NXdj6iqS5I8P8k3T8+9p7sftcNL/2ySf5DkhiTXJLkoyX/azKcA2D69BhiDXgOMQa8X7EkNzMbhvt+eb7vw6CSHuvum7v54klcmufiIdS5OcuV0/6okF97XN4FVdWaSz+ru3+7uTvKyJE86zo8NMBy9BhiDXgOMQa9NUgMz0akc7r3fkpxWVTcu3Z55xFucleTmpce3TMt2XKe7707yoSQPm547r6reUlX/paq+emn9W47xmgAHil4DjEGvAcag1wtO9wHMxooXCvhAd1+wrrEc4bYkD+/uO6rqy5P8alV9wYbeC2D29BpgDHoNMAa9tic1cOK4Nck5S4/PnpbtuE5VnZzkIUnu6O6PdfcdSdLdb07yniSfP61/9jFeE4Djo9cAY9BrgDEM0WuT1MAsdLLq4S3H8qYk51fVeVV1/ySXJLn6iHWuTnLpdP/JSV7f3V1Vp08XGkhVfW6S85Pc1N23JfmzqvrK6VxN35LktSv/ZQDMmF4DjEGvAcag1wsrne6jqv5xkmdk8ff5e0m+LcmZWZyA+2FJ3pzk6dNJuQHu0y5P+L8n3X13VX13kmuTnJTkiu5+e1U9N8mN3X11kpck+YWqOpTkzizCnSRfk+S5VfWXSQ4n+Y7uvnN67juTvDTJA7O4iu0srzyu18A66fXm6DWwTnq9OXoNrJNerzBJXVVnJfneJI/s7r+oqldn8QG+IckLuvuVVfXvklyW5GdXGSRwAtj9N4B7f4vua5Jcc8SyZy/d/2iSb9phu9ckec1RXvPGJF+43pGul14Da6XXG6PXwFrp9cboNbBWep1k9dN9nJzkgdO5Sh6Uxcm0H5Pkqun5K5M8acX3AGB1eg0wBr0GGINeA6zRniepu/vWJD+Z5I+ziPGHsjic5YPdffe02i1Jzlp1kMDB11lczXavN45Or4F10uvN0WtgnfR6c/QaWCe9XtjzJHVVnZrk4iTnJfkrST4jyUXHsf0zq+rGqrrx9ttv3+swgANkwxcKOGHpNbBuer0Zeg2sm15vhl4D66bXq53u47FJ/rC7b+/uv0zyK0m+Kskp0+EuSXJ2klt32ri7L+/uC7r7gtNPP32FYQAHwT5czfZEptfA2uj1Ruk1sDZ6vVF6DayNXi+sMkn9x0m+sqoeVFWV5MIk70jyhiRPnta5NMlrVxsicKIQ5Y3Ra2Ct9Hpj9BpYK73eGL0G1kqvVzsn9Q1ZXBDgd5P83vRalyf5Z0m+v6oOJXlYkpesYZwA7JFeA4xBrwHGoNcA63fysVc5uu5+TpLnHLH4piSPXuV1gRNP52B9Azg3eg2si15vll4D66LXm6XXwLro9cJKk9QA63SQrkoLcJDpNcAY9BpgDHptkhqYi45vDgFGoNcAY9BrgDHodRKT1MBM3Hs1WwDmTa8BxqDXAGPQ64U9XzgRAAAAAABWZU9qYDZ8cwgwBr0GGINeA4xBr01SAzPharYAY9BrgDHoNcAY9HrBJDUwGy3KAEPQa4Ax6DXAGPTaOakBAAAAANgie1IDs3E4vjkEGIFeA4xBrwHGoNcmqYGZ6HahAIAR6DXAGPQaYAx6vWCSGpgN52ACGINeA4xBrwHGoNcmqYHZcDVbgDHoNcAY9BpgDHqduHAiAAAAAABbZE9qYDYc3gIwBr0GGINeA4xBr01SAzPRcaEAgBHoNcAY9BpgDHq9YJIamIdeXNEWgJnTa4Ax6DXAGPQ6iUlqYEYOxzeHACPQa4Ax6DXAGPTahRMBAAAAANgie1IDs9BxoQCAEeg1wBj0GmAMer1gkhqYiXKhAIAh6DXAGPQaYAx6nZikBmbEhQIAxqDXAGPQa4Ax6LVzUgMAAAAAsEX2pAZmwzmYAMag1wBj0GuAMei1SWpgJrpFGWAEeg0wBr0GGINeL5ikBmbDhQIAxqDXAGPQa4Ax6LVzUgMzsvj2cG+33aiqi6rqXVV1qKqetcPzD6iqV03P31BV507Lv66q3lxVvzf9+Zilbd44veZbp9tnr+mvA2C29BpgDHoNMAa9tic1cIKoqpOSvCjJ1yW5Jcmbqurq7n7H0mqXJbmrux9RVZckeX6Sb07ygSR/p7v/pKq+MMm1Sc5a2u5p3X3jvnwQgANOrwHGoNcAYxil1/akBmaju/Z824VHJznU3Td198eTvDLJxUesc3GSK6f7VyW5sKqqu9/S3X8yLX97kgdW1QPW8JEBhqTXAGPQa4Ax6LVJamAmOnsP8hTl06rqxqXbM494i7OS3Lz0+JZ88rd/n7ROd9+d5ENJHnbEOn83ye9298eWlv38dGjLv6gqJ5ICDjS9BhiDXgOMQa8XnO4DmI1dnkrpaD7Q3ResZyQ7q6ovyOKQl8ctLX5ad99aVQ9O8pokT0/ysk2OA2Db9BpgDHoNMAa9tic1MBe98cNbbk1yztLjs6dlO65TVScneUiSO6bHZyf5D0m+pbvf84lhd986/fnhJL+UxWE0AAeXXgOMQa8BxqDXSUxSAyeONyU5v6rOq6r7J7kkydVHrHN1kkun+09O8vru7qo6JcmvJ3lWd//WvStX1clVddp0/9OSPCHJ72/2YwAceHoNMAa9BhjDEL12ug9gPlY8vuU+X7r77qr67iyuRHtSkiu6++1V9dwkN3b31UlekuQXqupQkjuzCHeSfHeSRyR5dlU9e1r2uCQfSXLtFOSTkrwuyc9t7lMAzIReA4xBrwHGoNcmqYH52OVhKiu8fl+T5Jojlj176f5Hk3zTDtv9eJIfP8rLfvk6xwgwAr0GGINeA4xBr01SAzPSG/zmEID10WuAMeg1wBj02jmpAQAAAADYIntSA7PQ2fzhLQCsTq8BxqDXAGPQ6wWT1MA8dBJRBpg/vQYYg14DjEGvk5ikBmbEOZgAxqDXAGPQa4Ax6LVJamBORBlgDHoNMAa9BhiDXq924cSqOqWqrqqqP6iqd1bV36yqh1bVdVX17unPU9c1WAD2Rq8BxqDXAGPQa4D1WmmSOslPJ/mN7v5rSb4kyTuTPCvJ9d19fpLrp8cAx1Dp3vuNY9JrYE30esP0GlgTvd4wvQbWRK+TFSapq+ohSb4myUuSpLs/3t0fTHJxkiun1a5M8qTVhgicMHqFG0el18Da6fVG6DWwdnq9EXoNrJ1er7Qn9XlJbk/y81X1lqp6cVV9RpIzuvu2aZ33JTljp42r6plVdWNV3Xj77bevMAzgQOj45nBz9BpYH73eJL0G1kevN0mvgfXR6ySrTVKfnOTLkvxsd39pko/kiENZuvuoc/rdfXl3X9DdF5x++ukrDAOAY9BrgDHoNcAY9BpgzVaZpL4lyS3dfcP0+KosIv2nVXVmkkx/vn+1IQInDIe3bIpeA+ul15ui18B66fWm6DWwXnq990nq7n5fkpur6q9Oiy5M8o4kVye5dFp2aZLXrjRC4ARSK9w4Gr0G1k+vN0GvgfXT603Qa2D99PrkFbf/niQvr6r7J7kpybdlMfH96qq6LMl7kzxlxfcAThQH6BvAGdJrYH30epP0Glgfvd4kvQbWR69Xm6Tu7rcmuWCHpy5c5XWBE5Qob4xeA2ul1xuj18Ba6fXG6DWwVnq90jmpAQAAAABgJaue7gNgPTpJH5xzKQEcWHoNMAa9BhiDXicxSQ3MSDu8BWAIeg0wBr0GGINem6QG5kSUAcag1wBj0GuAMei1SWpgRhzeAjAGvQYYg14DjEGvXTgRAAAAAIDtsSc1MBvl8BaAIeg1wBj0GmAMem2SGpiLjnMwAYxArwHGoNcAY9DrJCapgdko52ACGIJeA4xBrwHGoNeJc1IDAAAAALBF9qQG5sPhLQBj0GuAMeg1wBj02iQ1MCOiDDAGvQYYg14DjEGvTVIDMyLKAGPQa4Ax6DXAGPTaJDUwEx0XCgAYgV4DjEGvAcag10lcOBEAAAAAgC2yJzUwG+XwFoAh6DXAGPQaYAx6bU9qYE56hdsuVNVFVfWuqjpUVc/a4fkHVNWrpudvqKpzl577oWn5u6rq8bt9TYADSa8BxqDXAGPQa5PUwImhqk5K8qIkX5/kkUmeWlWPPGK1y5Lc1d2PSPKCJM+ftn1kkkuSfEGSi5L826o6aZevCcBx0GuAMeg1wBhG6bVJamA2qvd+24VHJznU3Td198eTvDLJxUesc3GSK6f7VyW5sKpqWv7K7v5Yd/9hkkPT6+3mNQEOHL0GGINeA4xBr01SAyeOs5LcvPT4lmnZjut0991JPpTkYfex7W5eE4Djo9cAY9BrgDEM0WsXTgTmo2uVrU+rqhuXHl/e3ZevOCIAdqLXAGPQa4Ax6LVJamAmjuOE/0fxge6+4D6evzXJOUuPz56W7bTOLVV1cpKHJLnjGNse6zUBDha9BhiDXgOMQa+TON0HMCebvZrtm5KcX1XnVdX9szjx/9VHrHN1kkun+09O8vru7mn5JdPVbs9Lcn6S39nlawIcPHoNMAa9BhiDXtuTGjgxdPfdVfXdSa5NclKSK7r77VX13CQ3dvfVSV6S5Beq6lCSO7OIbKb1Xp3kHUnuTvJd3X1Pkuz0mvv92QAOEr0GGINeA4xhlF6bpAZmY5dXpd2z7r4myTVHLHv20v2PJvmmo2z7E0l+YjevCXDQ6TXAGPQaYAx6bZIamJMNRxmANdFrgDHoNcAY9NokNTAjogwwBr0GGINeA4xBr01SA/NQvfnDWwBYnV4DjEGvAcag1wv32/YAAAAAAAA4cdmTGpiPrm2PAIDd0GuAMeg1wBj02iQ1MCMObwEYg14DjEGvAcag1yapgflwDiaAMeg1wBj0GmAMem2SGpgTUQYYg14DjEGvAcag1y6cCAAAAADA9tiTGpiHdngLwBD0GmAMeg0wBr1OYpIamBNRBhiDXgOMQa8BxqDXJqmBGRFlgDHoNcAY9BpgDHrtnNQAAAAAAGyPPamB2XAOJoAx6DXAGPQaYAx6vYY9qavqpKp6S1X9x+nxeVV1Q1UdqqpXVdX9Vx8mAKvSa4Ax6DXAGPQaYH3WcbqP70vyzqXHz0/ygu5+RJK7kly2hvcATgS9wo3d0GtgPfR60/QaWA+93jS9BtZDr1ebpK6qs5N8Y5IXT48ryWOSXDWtcmWSJ63yHsAJoheHt+z1xn3Ta2Bt9Hqj9BpYG73eKL0G1kavk6y+J/ULk/xgksPT44cl+WB33z09viXJWSu+BwCre2H0GmAEL4xeA4zghdFrgLXZ8yR1VT0hyfu7+8173P6ZVXVjVd14++2373UYwEHi8JaN0Gtg7fR6I/QaWDu93gi9BtZOr3PyCtt+VZInVtU3JPn0JJ+V5KeTnFJVJ0/fHp6d5NadNu7uy5NcniQXXHDBAforBfZMCTZFr4H1UoJN0WtgvZRgU/QaWC8l2Pue1N39Q919dnefm+SSJK/v7qcleUOSJ0+rXZrktSuPEjjwKs7BtCl6DayTXm+OXgPrpNebo9fAOun1wqrnpN7JP0vy/VV1KItzMr1kA+8BHEQOb9lveg3sjV7vN70G9kav95teA3uj1yud7uMTuvuNSd443b8pyaPX8boArJdeA4xBrwHGoNcA67GWSWqAlR2ww1QADiy9BhiDXgOMQa+TmKQG5kSUAcag1wBj0GuAMei1SWpgRkQZYAx6DTAGvQYYg15v5MKJAAAAAACwK/akBmbDOZgAxqDXAGPQa4Ax6LVJamBORBlgDHoNMAa9BhiDXpukBmaiI8oAI9BrgDHoNcAY9DqJSWpgRhzeAjAGvQYYg14DjEGvXTgRAAAAAIAtMkkNzEevcFtBVT20qq6rqndPf556lPUundZ5d1VdOi17UFX9elX9QVW9vaqet7T+t1bV7VX11un2jNVGCjATeg0wBr0GGINem6QG5qN677cVPSvJ9d19fpLrp8efPLaqhyZ5TpKvSPLoJM9ZivdPdvdfS/KlSb6qqr5+adNXdfejptuLVx4pwAzoNcAY9BpgDHptkhqYky19c5jk4iRXTvevTPKkHdZ5fJLruvvO7r4ryXVJLuruP+/uNyRJd388ye8mOXvlEQHMmV4DjEGvAcag1yapgZlYJciLKJ9WVTcu3Z55HO9+RnffNt1/X5IzdljnrCQ3Lz2+ZVr2CVV1SpK/k8W3j/f6u1X1tqq6qqrOOY4xAcyTXgOMQa8BxqDXSZKTj2PQAHP2ge6+4GhPVtXrknzODk/98PKD7u6q4z9gpqpOTvKKJD/T3TdNi38tySu6+2NV9e1ZfCv5mON9bYADRq8BxqDXAGM4EL02SQ3MQk23Tenuxx71vav+tKrO7O7bqurMJO/fYbVbk3zt0uOzk7xx6fHlSd7d3S9ces87lp5/cZJ/efwjB5gXvQYYg14DjEGvF5zuA5iP7Z2D6eokl073L03y2h3WuTbJ46rq1OkCAY+blqWqfjzJQ5L8o+UNpsDf64lJ3rnySAHmQK8BxqDXAGPQa3tSA/OxhqvS7tXzkry6qi5L8t4kT0mSqrogyXd09zO6+86q+rEkb5q2ee607OwsDpH5gyS/W1VJ8m+mK9d+b1U9McndSe5M8q37+aEANkWvAcag1wBj0GuT1AD3HoZy4Q7Lb0zyjKXHVyS54oh1bslRjszp7h9K8kNrHSzACUyvAcag1wBjmFOvTVID87G9bw4BOB56DTAGvQYYg16bpAZmRJQBxqDXAGPQa4Ax6LVJamAmeqvnYAJgt/QaYAx6DTAGvU5ikhqYE1EGGINeA4xBrwHGoNe537YHAAAAAADAicue1MBsOLwFYAx6DTAGvQYYg16bpAbmRJQBxqDXAGPQa4Ax6LVJamA+fHMIMAa9BhiDXgOMQa+dkxoAAAAAgC2yJzUwDx2HtwCMQK8BxqDXAGPQ6yQmqYE5EWWAMeg1wBj0GmAMem2SGpiHinMwAYxArwHGoNcAY9DrBZPUwHyIMsAY9BpgDHoNMAa9duFEAAAAAAC2x57UwGxU++oQYAR6DTAGvQYYg16bpAbmwtVsAcag1wBj0GuAMeh1EpPUwIy4UADAGPQaYAx6DTAGvTZJDcyJKAOMQa8BxqDXAGPQaxdOBAAAAABge+xJDcyGw1sAxqDXAGPQa4Ax6LVJamBORBlgDHoNMAa9BhiDXpukBmaifXMIMAS9BhiDXgOMQa+TrHBO6qo6p6reUFXvqKq3V9X3TcsfWlXXVdW7pz9PXd9wATheeg0wBr0GGINeA6zfKhdOvDvJP+nuRyb5yiTfVVWPTPKsJNd39/lJrp8eAxxbr3Djvug1sF56vSl6DayXXm+KXgPrpdd7n6Tu7tu6+3en+x9O8s4kZyW5OMmV02pXJnnSimMETgCVxeEte71xdHoNrJNeb45eA+uk15uj18A66fXCWs5JXVXnJvnSJDckOaO7b5ueel+SM46yzTOTPDNJHv7wh69jGMDo+gDVdab0GlgLvd44vQbWQq83Tq+BtdDrlU73kSSpqs9M8pok/6i7/2z5ue4+6o7n3X15d1/Q3Recfvrpqw4DOAB8c7hZeg2si15vll4D66LXm6XXwLro9YqT1FX1aVkE+eXd/SvT4j+tqjOn589M8v7VhgjAqvQaYAx6DTAGvQZYrz1PUldVJXlJknd2979eeurqJJdO9y9N8tq9Dw84YfSKN45Kr4G10uuN0WtgrfR6Y/QaWCu9TrLaOam/KsnTk/xeVb11WvbPkzwvyaur6rIk703ylJVGCJww6vC2R3Bg6TWwVnq9MXoNrJVeb4xeA2ul1ytMUnf3f8viApQ7uXCvrwucwA7QN4BzotfA2un1Rug1sHZ6vRF6DaydXq+0JzXAWh2kE/4DHGR6DTAGvQYYg16veOFEgIOgqh5aVddV1bunP089ynqXTuu8u6ouXVr+xqp6V1W9dbp99rT8AVX1qqo6VFU3VNW5+/SRAA4kvQYYg14DjGFOvTZJDcxDJ+ne+201z0pyfXefn+T66fEnqaqHJnlOkq9I8ugkzzki3k/r7kdNt3uv4n1Zkru6+xFJXpDk+asOFGDr9BpgDHoNMAa9TmKSGpiR6r3fVnRxkiun+1cmedIO6zw+yXXdfWd335XkuiQXHcfrXpXkwulK4ABD02uAMeg1wBj02iQ1MCe9wi05rapuXLo98zje+Yzuvm26/74kZ+ywzllJbl56fMu07F4/Px3a8i+WwvuJbbr77iQfSvKw4xgXwDzpNcAY9BpgDHrtwonAgfGB7r7gaE9W1euSfM4OT/3w8oPu7qrj/i7yad19a1U9OMlrkjw9ycuO8zUAThR6DTAGvQYYw4HotUlqYBYqm72abXc/9qjvXfWnVXVmd99WVWcmef8Oq92a5GuXHp+d5I3Ta986/fnhqvqlLM7R9LJpm3OS3FJVJyd5SJI7Vv80ANuj1wBj0GuAMej1gtN9APOwykUCVr9QwNVJLp3uX5rktTusc22Sx1XVqdMFAh6X5NqqOrmqTkuSqvq0JE9I8vs7vO6Tk7y+e/XBAmyVXgOMQa8BxqDXSexJDczIJr85PIbnJXl1VV2W5L1JnpIkVXVBku/o7md0951V9WNJ3jRt89xp2WdkEedPS3JSktcl+blpnZck+YWqOpTkziSX7N9HAtgcvQYYg14DjEGvTVIDc7KlKHf3HUku3GH5jUmesfT4iiRXHLHOR5J8+VFe96NJvmmtgwWYA70GGINeA4xBr53uAwAAAACA7bEnNTAbWzy8BYDjoNcAY9BrgDHotUlqYC46yWFVBpg9vQYYg14DjEGvk5ikBuZEkwHGoNcAY9BrgDHotUlqYD4c3gIwBr0GGINeA4xBr104EQAAAACALbInNTAf7atDgCHoNcAY9BpgDHptkhqYD4e3AIxBrwHGoNcAY9Brk9TAXHRcKABgBHoNMAa9BhiDXidxTmoAAAAAALbIntTALFSScg4mgNnTa4Ax6DXAGPR6wSQ1MB+Htz0AAHZFrwHGoNcAY9Brk9TAfPjmEGAMeg0wBr0GGINem6QG5sKFAgDGoNcAY9BrgDHodRIXTgQAAAAAYIvsSQ3MRCcObwEYgF4DjEGvAcag14lJamBGSpMBhqDXAGPQa4Ax6LVJamBOfHMIMAa9BhiDXgOMQa9NUgMz0Ukd3vYgADgmvQYYg14DjEGvk7hwIgAAAAAAW2RPamA+HN4CMAa9BhiDXgOMQa9NUgMzoskAY9BrgDHoNcAY9NokNTAf5ZtDgCHoNcAY9BpgDHrtnNQAAAAAAGyRPamB+fDNIcAY9BpgDHoNMAa9NkkNzEQnObztQQBwTHoNMAa9BhiDXicxSQ3MRKWdgwlgAHoNMAa9BhiDXi+YpAbmQ5QBxqDXAGPQa4Ax6LULJwIAAAAAsD0bmaSuqouq6l1VdaiqnrWJ9wAOoO6939gTvQb2RK/3nV4De6LX+06vgT3R6/Wf7qOqTkryoiRfl+SWJG+qqqu7+x3rfi/gAHGhgH2n18Ce6PW+02tgT/R63+k1sCd6nWQz56R+dJJD3X1TklTVK5NcnESUgfvkQgH7Tq+BPdHrfafXwJ7o9b7Ta2BP9Hozp/s4K8nNS49vmZYB3LctHd5SVQ+tquuq6t3Tn6ceZb1Lp3XeXVWXTsseXFVvXbp9oKpeOD33rVV1+9Jzz1hpoOun18De6PV+02tgb/R6v+k1sDd6vb0LJ1bVM6vqxqq68fbbb9/WMACS5FlJru/u85NcPz3+JFX10CTPSfIVWewh8ZyqOrW7P9zdj7r3luS9SX5ladNXLT3/4o1/kg3Qa2BG9Po+6DUwI3p9H/QamJHZ9HoTk9S3Jjln6fHZ07JP0t2Xd/cF3X3B6aefvoFhAGNZ4VvD1Q+LuTjJldP9K5M8aYd1Hp/kuu6+s7vvSnJdkouWV6iqz0/y2Ul+c9UB7RO9BvZAr7dAr4E90Ost0GtgD/Q62cwk9ZuSnF9V51XV/ZNckuTqDbwPcJB0Vo3yaffujTDdnnkc735Gd9823X9fkjN2WGc3h+5dksU3hcv/SvzdqnpbVV1VVedkXvQaOH56vQ16DRw/vd4GvQaOn14n2cCFE7v77qr67iTXJjkpyRXd/fZ1vw9wAK12NdsPdPcFR3uyql6X5HN2eOqHlx90d1fVXr+KvCTJ05ce/1qSV3T3x6rq27P4VvIxe3zttdNrYM/0el/pNbBner2v9BrYM71e/yR1knT3NUmu2cRrA+xFdz/2aM9V1Z9W1ZndfVtVnZnk/TusdmuSr116fHaSNy69xpckObm737z0nncsrf/iJP9yb6PfHL0G5kavd6bXwNzo9c70GpibUXq9tQsnAhypuvd8W9HVSS6d7l+a5LU7rHNtksdV1anT1W4fNy2711OTvOKTPs8i8Pd6YpJ3rjpQgDnQa4Ax6DXAGPR6Q3tSA+zJ6nHdq+cleXVVXZbF1WifkiRVdUGS7+juZ3T3nVX1Y1mcZy5Jntvddy69xlOSfMMRr/u9VfXEJHcnuTPJt27wMwDsH70GGINeA4xBr01SAzPRSQ5vJ8rTYSgX7rD8xiTPWHp8RZIrjvIan7vDsh9K8kPrGynADOg1wBj0GmAMep3EJDUwG73Nbw4B2DW9BhiDXgOMQa8T56QGAAAAAGCL7EkNzIdvDgHGoNcAY9BrgDHotUlqYEZEGWAMeg0wBr0GGINem6QGZmKLFwoA4DjoNcAY9BpgDHqdxCQ1MBud9OFtDwKAY9JrgDHoNcAY9Dpx4UQAAAAAALbIntTAfDgHE8AY9BpgDHoNMAa9NkkNzIRzMAGMQa8BxqDXAGPQ6yQmqYE58c0hwBj0GmAMeg0wBr12TmoAAAAAALbHntTAfPjmEGAMeg0wBr0GGINem6QG5qJFGWAIeg0wBr0GGINeJyapgbnoJIcPb3sUAByLXgOMQa8BxqDXSUxSA3Pim0OAMeg1wBj0GmAMeu3CiQAAAAAAbI89qYH58M0hwBj0GmAMeg0wBr02SQ3MRSeHRRlg/vQaYAx6DTAGvU5MUgNz0Um3CwUAzJ5eA4xBrwHGoNdJnJMaAAAAAIAtsic1MB8ObwEYg14DjEGvAcag1yapgRlxoQCAMeg1wBj0GmAMem2SGpiJ7uSwczABzJ5eA4xBrwHGoNdJTFIDc+KbQ4Ax6DXAGPQaYAx67cKJAAAAAABsjz2pgdloh7cADEGvAcag1wBj0GuT1MBstMNbAIag1wBj0GuAMeh1YpIamItOcliUAWZPrwHGoNcAY9DrJCapgTlph7cADEGvAcag1wBj0GsXTgQAAAAAYHvsSQ3MQidph7cAzJ5eA4xBrwHGoNcL9qQG5qF7cXjLXm8rqKqHVtV1VfXu6c9Tj7Leb1TVB6vqPx6x/LyquqGqDlXVq6rq/tPyB0yPD03Pn7vSQAHmQK8BxqDXAGPQ6yQmqYEZ6cO959uKnpXk+u4+P8n10+Od/KskT99h+fOTvKC7H5HkriSXTcsvS3LXtPwF03oAw9NrgDHoNcAY9NokNUCSXJzkyun+lUmetNNK3X19kg8vL6uqSvKYJFftsP3y616V5MJpfQD2Rq8BxqDXAGOYTa+dkxqYj+1dzfaM7r5tuv++JGccx7YPS/LB7r57enxLkrOm+2cluTlJuvvuqvrQtP4HVh8ywBbpNcAY9BpgDHo9j0nqN7/5zR+oqo9k/H9YTsvYn2H08Sfjf4bRx/+/7HXDD+eua1/XV522wnt/elXduPT48u6+/N4HVfW6JJ+zw3Y/vPygu7uqXLHgKPR6NkYffzL+Zxh9/Hp9wOn1bIw+/mT8zzD6+PX6gNPr2Rh9/Mn4n2Hk8e+51Yle32sWk9TdfXpV3djdF2x7LKsY/TOMPv5k/M8w+vhX0d0Xbfj1H3u056rqT6vqzO6+rarOTPL+43jpO5KcUlUnT98enp3k1um5W5Ock+SWqjo5yUOm9Yel1/Mw+viT8T/D6ONfhV6PQa/nYfTxJ+N/htHHvwq9HoNez8Po40/G/wyjj38Ver3gnNQAydVJLp3uX5rktbvdsLs7yRuSPHmH7Zdf98lJXj+tD8De6DXAGPQaYAyz6bVJaoDkeUm+rqreneSx0+NU1QVV9eJ7V6qq30zyy1mc8P+Wqnr89NQ/S/L9VXUoi3MsvWRa/pIkD5uWf3+OfpVcAHZHrwHGoNcAY5hNr2dxuo/J5cdeZfZG/wyjjz8Z/zOMPv4hdfcdSS7cYfmNSZ6x9Pirj7L9TUkevcPyjyb5pvWNdDYOwn+no3+G0cefjP8ZRh//kPT6uB2E/05H/wyjjz8Z/zOMPv4h6fVxOwj/nY7+GUYffzL+Zxh9/EOaU6/LkTEAAAAAAGyL030AAAAAALA1s5ikrqqLqupdVXWoqmZ/TqmqOqeq3lBV76iqt1fV903LH1pV11XVu6c/T932WI+lqk6qqrdU1X+cHp9XVTdMP4tXVdX9tz3Go6mqU6rqqqr6g6p6Z1X9zdF+BlX1j6f/hn6/ql5RVZ8+0s+AE49eb8fIrU70GrZBr7dDr7dPrxmNXm+HXm+XVrOTrU9SV9VJSV6U5OuTPDLJU6vqkdsd1THdneSfdPcjk3xlku+axvysJNd39/lJrs8YF3H4viTvXHr8/CQv6O5HJLkryWVbGdXu/HSS3+juv5bkS7L4HMP8DKrqrCTfm+SC7v7CJCcluSRj/Qw4gej1Vo3c6kSvYV/p9Vbp9RbpNaPR663S6y3Rao5m65PUWZxc+1B339TdH0/yyiQXb3lM96m7b+vu353ufziLGJyVxbivnFa7MsmTtjLAXaqqs5N8Y5IXT48ryWOSXDWtMtvPUFUPSfI1ma4a2t0f7+4PZrCfQRYXL31gVZ2c5EFJbssgPwNOSHq9BSO3OtFr2BK93gK9ng29ZiR6vQV6PQtazaeYwyT1WUluXnp8y7RsCFV1bpIvTXJDkjO6+7bpqfclOWNb49qlFyb5wSSHp8cPS/LB7r57ejznn8V5SW5P8vPTITovrqrPyEA/g+6+NclPJvnjLIL8oSRvzjg/A048er0dL8y4rU70GrZBr7fjhdHrrdJrBqTX2/HC6PXWaDVHM4dJ6mFV1WcmeU2Sf9Tdf7b8XHd3kt7KwHahqp6Q5P3d/eZtj2WPTk7yZUl+tru/NMlHcsShLAP8DE7N4pvO85L8lSSfkeSirQ4KDqhRe30AWp3oNXAc9Hqr9BrYNb3eqqF7rdUczRwmqW9Ncs7S47OnZbNWVZ+WRZBf3t2/Mi3+06o6c3r+zCTv39b4duGrkjyxqv4oi0OKHpPFOY1OmQ63SOb9s7glyS3dfcP0+KosIj3Sz+CxSf6wu2/v7r9M8itZ/FxG+Rlw4tHr/Td6qxO9hm3Q6/2n1/Og14xGr/efXm+fVrOjOUxSvynJ+dNVPO+fxcnSr97ymO7TdL6ilyR5Z3f/66Wnrk5y6XT/0iSv3e+x7VZ3/1B3n93d52bxd/767n5akjckefK02mw/Q3e/L8nNVfVXp0UXJnlHBvoZZHFoy1dW1YOm/6bu/QxD/Aw4Ien1Phu91Ylew5bo9T7T69nQa0aj1/tMr2dBq9lRLY4A2PIgqr4hi3MCnZTkiu7+ie2O6L5V1f+a5DeT/F7+/3MY/fMszsP06iQPT/LeJE/p7ju3MsjjUFVfm+QHuvsJVfW5WXyb+NAkb0ny97v7Y1sc3lFV1aOyuNDB/ZPclOTbsvjiZZifQVX9aJJvzuIKyW9J8owszrs0xM+AE49eb8+orU70GrZBr7dHr7dLrxmNXm+PXm+PVrOTWUxSAwAAAABwYprD6T4AAAAAADhBmaQGAAAAAGBrTFIDAAAAALA1JqkBAAAAANgak9QAAAAAAGyNSWoAAAAAALbGJDUAAAAAAFtjkhoAAAAAgK35/wABuAE/8aaljwAAAABJRU5ErkJggg==\n",
779
780
781
782
      "text/plain": [
       "<Figure size 1800x432 with 8 Axes>"
      ]
     },
783
784
785
     "metadata": {
      "needs_background": "light"
     },
786
787
788
789
790
791
792
793
794
795
796
     "output_type": "display_data"
    }
   ],
   "source": [
    "init_drop()\n",
    "time_loop(1000)\n",
    "plot_status()"
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Markus Holzer committed
797
   "execution_count": 21,
798
   "metadata": {},
Markus Holzer's avatar
Markus Holzer committed
799
800
801
802
803
804
805
806
807
808
809
810
811
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-21-3f4503d1cf76>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mρ_arr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdh\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgather_array\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mρ_field\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misnan\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mρ_arr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0;32massert\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misnan\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mρ_arr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAssertionError\u001b[0m: "
     ]
    }
   ],
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
   "source": [
    "ρ_arr = dh.gather_array(ρ_field.name)\n",
    "assert not np.isnan(np.min(ρ_arr))\n",
    "assert not np.isnan(np.max(ρ_arr))"
   ]
  }
 ],
 "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",
Markus Holzer's avatar
Markus Holzer committed
835
   "version": "3.8.2"
836
837
838
839
840
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}