Commit 458f6d6b authored by Michael Kuron's avatar Michael Kuron
Browse files

correct importorskip in notebook tests

parent 3ac2bf38
%% Cell type:code id: tags:
``` python
import pytest
pytest.importorskip('scipy.ndimage')
```
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.analytical import *
from lbmpy.phasefield.eos import *
from lbmpy.phasefield.high_density_ratio_model import *
from pystencils.fd.derivative import replace_generic_laplacian
from pystencils.fd.spatial import discretize_spatial, fd_stencils_standard, fd_stencils_isotropic
from lbmpy.phasefield.fd_stencils import fd_stencils_isotropic_high_density_code
from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
from lbmpy.forcemodels import EDM
from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
from scipy.ndimage.filters import gaussian_filter
```
%% Cell type:markdown id: tags:
# Implementation of high density difference model
According to *"Ternary free-energy entropic lattice Boltzmann model with high density ratio" by Wöhrwag, Semprebon, Mazloomi, Karlin and Kusumaatmaja*
Up front we define all necessary parameters in one place:
%% Cell type:code id: tags:
``` python
a = 0.037
b = 0.2
reduced_temperature = 0.61
gas_constant = 1
κ = (0.01, 1, 1)
λ = (0.6, 1, 1)
χ = 5
φ_relaxation_rate = 1
ρ_relaxation_rate = 1
external_force = (0, 0)
clipping = False
domain_size = (100, 100)
stencil = get_stencil('D2Q9', ordering='uk')
fd_discretization = fd_stencils_isotropic_high_density_code
target = 'cpu'
threads = 4
```
%% Cell type:markdown id: tags:
## Part 1: Free Energy
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.
Here we use the Carnahan-Starling EOS, with the parametrization from the paper
%% Cell type:code id: tags:
``` python
ρ, φ, c_l1, c_l2 = sp.symbols("rho, phi, c_l1, c_l2")
critical_temperature = carnahan_starling_critical_temperature(a, b, gas_constant)
temperature = reduced_temperature * critical_temperature
eos = carnahan_starling_eos(ρ, gas_constant, temperature, a, b)
eos
```
%%%% Output: execute_result
![](data:image/png;base64,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)
$\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}}$
⎛ 3 2 ⎞
2 0.042578305⋅ρ⋅⎝- 0.000125⋅ρ + 0.0025⋅ρ + 0.05⋅ρ + 1⎠
- 0.037⋅ρ + ──────────────────────────────────────────────────────
3
(1 - 0.05⋅ρ)
%% Cell type:markdown id: tags:
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 id: tags:
``` python
ρ_g, ρ_l = maxwell_construction(eos, tolerance=1e-3)
(ρ_g, ρ_l, ρ_l / ρ_g)
```
%%%% Output: execute_result
![](data:image/png;base64,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==)
$\displaystyle \left( 0.0695273860315274, \ 8.02904149705209, \ 115.480272671423\right)$
(0.0695273860315274, 8.02904149705209, 115.480272671423)
%% Cell type:markdown id: tags:
With this information we call a function that assembles the full free energy density:
%% Cell type:code id: tags:
``` python
free_energy = free_energy_high_density_ratio(eos, ρ, ρ_g, ρ_l, c_l1, c_l2, λ, κ)
free_energy
```
%%%% Output: execute_result
![](data:image/png;base64,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=)
$\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$
2 2 2 2 ⎛ 1.0⋅(1.7031322⋅ρ
0.5⋅cₗ₁ ⋅(1 - cₗ₁) + 0.5⋅cₗ₂ ⋅(1 - cₗ₂) + 0.3⋅ρ⋅⎜-0.037⋅ρ - ────────────────
⎜ 2
⎝ 1.0⋅ρ - 4
- 51.0939660000001) ⎞
──────────────────── + 0.042578305⋅log(1.0⋅ρ) - 0.0530922164415325⎟ + 0.5⋅D(c_
0.0⋅ρ + 400.0 ⎠
2 2 2
l1) + 0.5⋅D(c_l2) + 0.005⋅D(rho) + 0.00085206629489322
%% Cell type:markdown id: tags:
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 id: tags:
``` python
transformation_eqs = [ c_l1 - (1 + φ/χ - (ρ - ρ_l)/(ρ_g - ρ_l)) / 2,
c_l2 - (1 - φ/χ - (ρ - ρ_l)/(ρ_g - ρ_l)) / 2]
transform_forward_substitutions = sp.solve(transformation_eqs, [c_l1, c_l2])
transform_backward_substitutions = sp.solve(transformation_eqs, [ρ, φ])
```
%% Cell type:markdown id: tags:
To do the transformation, we use the substitutions dict.
After the substitutions the differentials have to be expanded again.
%% Cell type:code id: tags:
``` python
free_energy_transformed = free_energy.subs(transform_forward_substitutions)
free_energy_transformed = expand_diff_full(free_energy_transformed, functions=(ρ, φ))
free_energy_transformed.atoms(sp.Symbol)
```
%%%% Output: execute_result
![](data:image/png;base64,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)
$\displaystyle \left\{\phi, \rho\right\}$
{φ, ρ}
%% Cell type:markdown id: tags:
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 id: tags:
## Part 2: Data setup
%% Cell type:code id: tags:
``` python
dh = create_data_handling(domain_size, periodicity=True, default_target=target)
# Fields for order parameters
ρ_field = dh.add_array("rho")
φ_field = dh.add_array("phi")
c_field = dh.add_array("c", values_per_cell=2)
# Chemical potential, pressure tensor, forces and velocities
μ_phi_field = dh.add_array("mu_phi", latex_name=r"\mu_{\phi}")
pbs_field = dh.add_array("pbs")
pressure_tensor_field = dh.add_array("p", len(symmetric_tensor_linearization(dh.dim)))
force_field = dh.add_array("force", values_per_cell=dh.dim, latex_name="F")
vel_field = dh.add_array("velocity", values_per_cell=dh.dim)
# PDF fields for lattice Boltzmann schemes
pdf_src_rho = dh.add_array("pdf_src_rho", values_per_cell=len(stencil))
pdf_dst_rho = dh.add_array_like("pdf_dst_rho", "pdf_src_rho")
pdf_src_phi = dh.add_array("pdf_src_phi", values_per_cell=len(stencil))
pdf_dst_phi = dh.add_array_like("pdf_dst_phi", "pdf_src_phi")
```
%% Cell type:markdown id: tags:
## Part 3a: Compute kernels and time loop
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 id: tags:
``` python
def make_kernel(assignments):
# assignments may be using the symbols ρ and φ
# these is substituted with the access to the corresponding fields here
field_substitutions = {
ρ: ρ_field.center,
φ: φ_field.center
}
processed_assignments = []
for a in assignments:
new_rhs = a.rhs.subs(field_substitutions)
# ∂∂f representing the laplacian of f is replaced by the explicit carteisan form
# ∂_0 ∂_0 f + ∂_1 ∂_1 f (example for 2D)
# otherwise the discretization would not do the correct thing
new_rhs = replace_generic_laplacian(new_rhs)
# Next the "∂" objects are replaced using finite differences
new_rhs = discretize_spatial(new_rhs, dx=1, stencil=fd_discretization)
processed_assignments.append(Assignment(a.lhs, new_rhs))
return create_kernel(processed_assignments, target=target, cpu_openmp=threads).compile()
```
%% Cell type:markdown id: tags:
#### Chemical Potential
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 id: tags:
``` python
μ_ρ, μ_φ = chemical_potentials_from_free_energy(free_energy_transformed,
order_parameters=(ρ, φ))
μ_phi_assignment = Assignment(μ_phi_field.center, μ_φ)
μ_kernel = make_kernel([μ_phi_assignment])
μ_phi_assignment
```
%%%% Output: execute_result
![](data:image/png;base64,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)
$\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}}$
3 2
μ_φ_C := 0.0004⋅φ + 0.000473530700220886⋅φ⋅ρ - 0.00760399528440326⋅φ⋅ρ + 0.0
205263968409311⋅φ - 0.02⋅D(D(phi))
%% Cell type:markdown id: tags:
#### Pressure tensor and force computation
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.
$ pbs = \sqrt{|ρ c_s^2 - p_{bulk} |} $
The force is then calculated as $ \nabla \cdot P_{if} + 2 (\nabla pbs) pbs$
In the following kernel the pressure tensor field is filled with $P_{if}$ and the pbs field with above expression.
%% Cell type:code id: tags:
``` python
# Bulk part
pressure_assignments = [
Assignment(pbs_field.center,
pressure_tensor_bulk_sqrt_term(free_energy_transformed, (ρ, φ), ρ)),
]
# Interface part
P_if = pressure_tensor_from_free_energy(free_energy_transformed, (ρ, φ),
dim=dh.dim, include_bulk=False)
index_map = symmetric_tensor_linearization(dh.dim)
pressure_assignments += [
Assignment(pressure_tensor_field(index_1d), P_if[index_2d])
for index_2d, index_1d in index_map.items()
]
pressure_kernel = make_kernel(pressure_assignments)
# Force kernel
pressure_tensor_sym = sp.Matrix(dh.dim, dh.dim,
lambda i, j: pressure_tensor_field(index_map[i, j]
if i < j else index_map[j, i]))
force_term = force_from_pressure_tensor(pressure_tensor_sym,
functions=[ρ, φ],
pbs=pbs_field.center)
force_assignments = [
Assignment(force_field(i),
force_term[i] + external_force[i] * ρ_field.center / ρ_l)
for i in range(dh.dim)
]
force_kernel = make_kernel(force_assignments)
```
%% Cell type:markdown id: tags:
#### Lattice Boltzmann schemes for time evolution of ρ and φ
- ρ is handled by a normal LB method (compressible, entropic equilibrium)
- stream and collide are splitted into separate kernels
- macroscopic values are computed after the stream, but inside the stream kernel
- velocity field stores the velocity which was not corrected for the forces yet
- 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
- when ρ and φ are updated, they are clipped to a valid region, this clipping should be only necessary during equilibration of the system
- exact difference method is used to couple the force into the ρ-LBM
%% Cell type:markdown id: tags:
The following cell handles the clipping of the order parameters:
%% Cell type:code id: tags:
``` python
if clipping:
def clip(ac, symbol, min_value, max_value):
"""Function to clip the value of a symbol which is on one of lhs of the assignments
in an assignment collection"""
assert symbol in ac.bound_symbols
for i in range(len(ac.subexpressions)):
a = ac.subexpressions[i]
if a.lhs == symbol:
new_assignment = Assignment(symbol, sp.Piecewise((max_value, a.rhs > max_value),
(min_value, a.rhs < min_value),
(a.rhs, True)))
ac.subexpressions[i] = new_assignment
break
# TODO: how can this 'densgin' be derived automatically?
tred = reduced_temperature
densgin = -67.098 \
+ 549.69 * tred \
- 1850.6 * tred * tred \
+ 3281 * tred * tred * tred \
- 3237.3 * tred * tred * tred * tred \
+ 1687.6 * tred * tred * tred * tred * tred \
- 361.51 * tred * tred * tred * tred * tred * tred
ρ_clip_min, ρ_clip_max = densgin * 0.5, 1.2 * ρ_l
φ_clip_min, φ_clip_max = -χ * 1.5, χ * 1.5
```
%% Cell type:markdown id: tags:
Next, the collide and stream kernels for the ρ lattice Boltzmann are created
%% Cell type:code id: tags:
``` python
force_model = EDM(force_field.center_vector)
ρ_lbm_params = {
'stencil': stencil,
'method' : 'trt-kbc-n2',
'compressible': True,
'relaxation_rate': ρ_relaxation_rate,
'optimization': {
'symbolic_field': pdf_src_rho,
'symbolic_temporary_field': pdf_dst_rho,
'openmp': threads,
'target': target
}
}
# Standard collision step, that does not compute ρ and u from pdfs, but reads
# them from fields - this is necessary because ρ may have been clipped before
# the velocity field is not force corrected, which is the correct for the EDM model
# but might be wrong for other force models
ρ_collide = create_lb_function(kernel_type='collide_only',
density_input=ρ_field,
force_model=force_model,
velocity_input=vel_field,
**ρ_lbm_params)
# First the assignments are created, then the density is clipped
# then a kernel is created from the clipped assignments
ρ_stream_ur = create_lb_update_rule(kernel_type='stream_pull_only',
force_model=None, #save uncorrected velocity
output={'density': ρ_field,
'velocity': vel_field},
**ρ_lbm_params)
if clipping:
clip(ρ_stream_ur, sp.Symbol("rho"), ρ_clip_min, ρ_clip_max)
ρ_stream = create_lb_function(update_rule=ρ_stream_ur, **ρ_lbm_params)
```
%% Cell type:markdown id: tags:
The φ lattice Boltzmann solve the Cahn-Hilliard equation.
%% Cell type:code id: tags:
``` python
φ_lb_method = cahn_hilliard_lb_method(stencil=stencil,
mu=μ_phi_field.center,
relaxation_rate=φ_relaxation_rate,
gamma=1)
φ_lbm_params = {
'lb_method' : φ_lb_method,
'compressible': True,
'optimization': {
'symbolic_field': pdf_src_phi,
'symbolic_temporary_field': pdf_dst_phi,
'openmp': threads,
'target': target
}
}
ch_method = cahn_hilliard_lb_method(stencil=stencil,
mu=μ_phi_field.center,
relaxation_rate=φ_relaxation_rate,
gamma=1)
corrected_vel = vel_field.center_vector + sp.Matrix(force_model.macroscopic_velocity_shift(ρ_field.center))
φ_collide = create_lb_function(kernel_type='collide_only',
density_input=φ_field,
velocity_input=corrected_vel,
**φ_lbm_params)
φ_stream_ur = create_lb_update_rule(kernel_type='stream_pull_only',
output={'density': φ_field},
**φ_lbm_params)
if clipping:
clip(φ_stream_ur, sp.Symbol("rho"), φ_clip_min, φ_clip_max)
φ_stream = create_lb_function(update_rule=φ_stream_ur, **φ_lbm_params)
```
%% Cell type:markdown id: tags:
#### Time loop
Now we can put all kernels together into a time loop function
%% Cell type:code id: tags:
``` python
op_sync = dh.synchronization_function([ρ_field.name, φ_field.name])
p_sync = dh.synchronization_function([pbs_field.name, pressure_tensor_field.name])
pdf_sync = dh.synchronization_function([pdf_src_phi.name, pdf_src_rho.name])
def time_loop(steps):
for t in range(steps):
op_sync()
dh.run_kernel(μ_kernel)
dh.run_kernel(pressure_kernel)
p_sync()
dh.run_kernel(force_kernel)
dh.run_kernel(ρ_collide)
dh.run_kernel(φ_collide)
pdf_sync()
dh.run_kernel(ρ_stream)
dh.run_kernel(φ_stream)
dh.swap(pdf_dst_phi.name, pdf_src_phi.name)
dh.swap(pdf_dst_rho.name, pdf_src_rho.name)
return dh.cpu_arrays[φ_field.name][1:-1, 1:-1]
```
%% Cell type:markdown id: tags:
## Part 3b: Compiling getter & setter kernels
%% Cell type:markdown id: tags:
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 id: tags:
``` python
init_assignments = [
Assignment( φ_field.center, transform_backward_substitutions[φ].subs({ c_l1: c_field(0), c_l2: c_field(1)} )),
Assignment( ρ_field.center, transform_backward_substitutions[ρ].subs({ c_l1: c_field(0), c_l2: c_field(1)} )),
]
init_rho = pdf_initialization_assignments(ρ_collide.method,
density=ρ_field.center,
velocity=vel_field.center_vector,
pdfs=pdf_src_rho.center_vector)
init_rho = init_rho.new_without_subexpressions()
init_assignments += init_rho.all_assignments
init_phi = pdf_initialization_assignments(φ_collide.method,
density=φ_field.center,
velocity=(0,0,0),
pdfs=pdf_src_phi.center_vector)
init_phi = init_phi.new_without_subexpressions().new_with_substitutions({μ_phi_field.center: 0})
init_assignments += init_phi.all_assignments
init_pdfs_assignments = init_rho.all_assignments + init_phi.all_assignments
init_pdfs_kernel = create_kernel(init_pdfs_assignments).compile()
init_kernel = create_kernel(init_assignments).compile()
```
%% Cell type:markdown id: tags:
## Part 4: Geometry initialization & plotting
%% Cell type:code id: tags:
``` python
def plot_status():
plt.figure(figsize=(25, 6))
plt.subplot(1, 4, 1)
ρ_arr = dh.cpu_arrays[ρ_field.name][1:-1, 1:-1]
plt.scalar_field(ρ_arr)
plt.title("ρ ({:.2f}, {:.2f})".format(np.min(ρ_arr), np.max(ρ_arr)))
plt.colorbar()
plt.subplot(1, 4, 2)
φ_arr = dh.cpu_arrays[φ_field.name][1:-1, 1:-1]
plt.scalar_field(φ_arr)
plt.title("φ ({:.2f}, {:.2f})".format(np.min(φ_arr), np.max(φ_arr)))
plt.colorbar()
plt.subplot(1, 4, 3)
f_arr = dh.cpu_arrays[force_field.name][1:-1, 1:-1]
plt.vector_field_magnitude(f_arr)
plt.title("F ({:.2f}, {:.2f})".format(np.min(f_arr), np.max(f_arr)))
plt.colorbar()
plt.subplot(1, 4, 4)
μ_arr = dh.cpu_arrays[μ_phi_field.name][1:-1, 1:-1]
plt.scalar_field(μ_arr)
plt.title("μ_φ ({:.2f}, {:.2f})".format(np.min(μ_arr), np.max(μ_arr)))
plt.colorbar()
def init_drop():
radius = dh.shape[0] // 5
mid1 = [dh.shape[0] // 2, dh.shape[1] // 2]
for block in dh.iterate(ghost_layers=True):
x, y = block.midpoint_arrays
mask1 = (x - mid1[0]) ** 2 + (y - mid1[1])**2 < radius ** 2
block[force_field.name].fill(0)
block[vel_field.name].fill(0)
block[μ_phi_field.name].fill(0)
c_arr = block[c_field.name]
c_arr[:, :].fill(0.0)
c_arr[mask1, 0] = 1.0
gaussian_filter(c_arr[..., 0], sigma=3, output=c_arr[..., 0])
gaussian_filter(c_arr[..., 1], sigma=3, output=c_arr[..., 1])
dh.run_kernel(init_kernel)
```
%% Cell type:markdown id: tags:
# Part 5: Putting it all together
%% Cell type:code id: tags:
``` python
init_drop()
time_loop(1000)
plot_status()
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABbAAAAF1CAYAAAA5ssyvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzde5w0aVnf/891V3XPzHPYEwsLLODu6uIBFFQEfxoUAyoSdfGAP0TlEA2SgPhS81I0RtBIQtT4QhMUNwYBQRBBDiYoKpFI/AmygKILooAcdllglz09p5nuqrp+f9x3Vdf09Glmume6Z77v16t3prurq6r7efbq+7nquq/b3B0RERERERERERERkWUTDvsERERERERERERERERGUQJbRERERERERERERJaSEtgiIiIiIiIiIiIispSUwBYRERERERERERGRpaQEtoiIiIiIiIiIiIgsJSWwRURERERERERERGQpKYF9jJnZPc3sg2a2ftjncljM7K/M7EGHfR4iIqOY2avM7PGHfR7LwsyuMLMPmNnaYZ+LiEjNzNbM7P1mdu/DPpdlYWbfamavPuzzEBGZlfIjO5nZs83sBYd9HiKgBPZx9xzgt9x9E5rB90vM7G4z+5SZ/eikF5vZk8zsY2Z2zszeYGaXtZ77BTP7RNrXx8zs3816Uhb9vJndbGZ3mdnbJiWZzez69EVTmdlTh557YnruLjP7jJm9zMwuam3yS8DPzXpuIiIHxcy+BHgI8MYJ21xlZn9mZufN7O/N7DETtjUz+89m9tl0+wUzs/TcvVKy/JMpXv6FmT1iF+f6UjPrmdnZ1i2bsP3Y74/0/GPM7D3p+U+Y2XcBuPungT8Dnj7ruYmI7IWZfdTMLgzFtfuO2fzpwJ+7+6cm7O/RKU6fT3H7cyZs+zYz22wd94Ot575n6JzOm5mb2ZfP+L52873xPDPrDx3vmvTcA4YeP5vO48cA3P1NwIPTd5mIyCrYlh8ZZdoYdmjbifF2N/sase+Hmtm7077fbWYPnbDtxDzPlH1dD3yvmd1r1nMTWRQlsI8pi9VrTwFe0Xr4ecC1wOcAXwf8uJk9dszrHwT8BvB9wBXAeeDXWpv8D+AL3P0i4KuAJ5nZt894ek8A/iXwSOAy4C+B356w/d8A/wZ4z4jn/gL4ane/GLgGyIGfbz3/JuDrzOw+M56biMhB+UHgle7uE7Z5FfBe4B7AvwNea2b3HLPt04HHE5PiXwJ8czoGwCngXcCXE+Puy4D/ZWandnG+v+Dup1q3ctRG074/zOyLgN9J7+di4KHAu1u7eGXrvEVEFulbhuLaJ8ds94NMGKua2eXA7wP/nhhjbwB+d8qxn9U67ufXD7r7K9vnRBwDf4TR4+BRdvO9AfC7Q5/BR9J5fHzoPL4YqIDXDR1LFxxFZOmNyY8MbzMtBzJsbLzdw77a59ElFri8AriUOG5/Y3p8lOcxJs8zbV8pmf+HwJNnOTeRRVIC+4gxsx8zs1vN7I50lW3c9JdHAHe6+02tx54M/Ad3v8PdPwD8d+CpY17/PcAfuPufu/tZ4oD8283sNIC7f9Ddz7W2r4DPm/FtXA38X3f/SEqAvAL4onEbu/uL3P2twI4rpe7+CXe/rfVQ2T6PFJDfDXzDjOcmIrInZvZVZvY3ZvY1KU6/OyVrx/km4P9M2N8DgS8DnuvuF9z9dcDfAt8x5iVPAf6Lu9/k7jcD/4UU41O8/WV3v8XdS3e/HugCnz9mX/sx8fsD+GngN9z9D929cPfPuvuHW69/J3DNpOpFEZGDYmYPAD6XGJvG+XbgRnf/vTT2fB7wEDP7gjmcwlOAl0+52Fmf626/N3bjycQq9I+2Hnsb8C/msG8ROcbS7I7Pa91/qZn9/KTXpO2ea2a/aGYvNLM7Lc7c7ozZfFR+ZNi0MWz72NPi7cz7GuFRxMK8F7r7lrv/KmDAPx+z/aQ8zyz7ehuK5bIElMA+QlKF878iBsrPJSYefnrM5l8MtKcjXgrcl1jNXPsbYFzrjge1t03JhR7wwNY+n2NmZ4GbgJPEirpZvBr4PDN7YPqCeQrwRzO+dgcz+2dmdhdwhviF8cKhTT5ArEgUEVkIMzsB/B7wIuKA9ZXEWPdKs9jGY2j7k8SLeR8cfq7lQcBH3P1M67GZ4/akbdPUwS7woQnHH/ZvzOz2lJiflAyZ9v3xlekc/tbMbjGzV7SnVLp7kc5LcVtElsEXE2NxMWGb4bh3Dvgw4+M1wH8ys9sstnR61KgN0oW8rwFePuO57vZ7A+BbUmy/0cz+9YTtnkys3Gv7AHCVbW/fJyKycGZ2HXGGyn8FngE8jpgnedaYl2zLj4wxNQcytO2keLubfY3a9/uGLly+jxGxfIY8zyz7Ur5EloIS2EfLDxCvnH3C3W8Hng9895htLyEmdGv1NPG7Wo/dBYy7AnhqaNsd27v7C9L9LyNOqxzefpxbgLcTv0AuEFuK/MiMr93B3f9vaiFyP+AXgY8ObXKG+HmIiCzKw4nVDP+deEHvTuKA+ouBB4zYvo5JZ0Y8V5sah6dsfxdwajiBnhINvw38rLvPGrd/lTg18V7ECpKXmtlX7/G870ecTvkdaZ8bxM+qTXFbRA7CG1LV3p1m9oYx2wyPqUfZbbz+CWLruyuJ/Uf/wMw+d8R2Twbe7u7/NOX4ez2P1wBfCNyTWCTzM2a2498WZvZI4hT41w49VX8uitcictC+Dfgt4JPAWvr5YuBbx2w/71g+bdvdxuP9nEf9/F7P4wyxrZ/IoVIC+2i5GvhE6/7HiFfbRrmD7UHpbPrZrpC4iPFB/OzQtiO39+i9xET0z4498+2eC3wFcH9gPb3uf6cKxj1LU+b/iFj12HaamEwSEVmUK4Cb29UNaRr5HcC9R2xfx6QmTqfqt3qhrEcyYxxuGd7+IuBs+5zMbAP4A+Ad7v6fZnpn8b28J7X6KNz9zcQK83HrHkw77wvEBXT+IU2p/I/Eqpk2xW0ROQiPd/dL0u3xY7bZNqa2ocUN08O7itfu/k53P5Omc7+MuKbLcByE0VXPk+z2PN7v7p9MraX+P+BXgO8cselTgNelmN1Wfy6K1yJy0K5ge24EYqHcqHE37Izlj2zF8hvTw7uJodO23e04fjf7Ht62fn6v53Ga2YsRRRZGCeyjZYOY9K09gHilcZT30Zqe4u53EAN6e2rIQ4AbGe3G9rYWVyRfA/5hzPY5sa3JLB5CXDDmppQMeSlxQYFJvWJnNeo8vpDtU2pERObtZuB+7WrnlCy+lNhmaZvW9PJ2nH5Qa8GstxPj8DVDvfJmjtvD26bFa96QznW/iyQ6seJ86nmM+P54X3r9SGaWE9cyUNwWkWXwPmIszmHk4oawM+6dJI5Hx8XrYTtiaprlcl92Vj1PstvvjVnOY4M4W3JUIv0LgY+6+927OEcRkVHaa3vNMqvjZnbOcrySEePuZDg/8vZWLK/baewmBzIt3u42nzK87y8ZmkX5JYyI5TPkeWbZl/IlshSUwD56nmlm90v9Qn+K8Suc/xVwiZld2Xrs5cBPm9mlaVGZfwW8dMzrX0nsiffINAj/OeD33f2MmQUz+8G0HzOzhwPPBN5av9jMPmpmTx2z73cBTzCzK9K+vg/oMKYXq5l102KVBnTMbN3MQnrue1IljKU+gc8fOo814MuBPxlzLiIi8/BO4BwxFhqQEWeX/EWaHTLKm4GvHbdDd/8H4K+B56a4923EAefrxrzk5cCPmtmVZnZf4MdIMT6tN/BaYvXzk929ar/QzK5KC+hcNWrHZvadZnYqxexvAL4XeNOY8xj7/ZGe/y3gaWZ2TZp58xPA/2y9/uHEhMjHxuxfROTApAW//pEYm8Z5PfBgM/uONGb9GWLP0b8f3tDMLjGzb0xxPTez7yH2uX7L0KZ11fOZodc/1cw+OuZcd/W9YWbXDY3nnw28cWizbyNWWP/ZiF18LfCHo/YtIrJLTzOzzOI6LY8GTtv4BRkhjjefSix6ALgP8MOMX5drVH5k1D4njWEbM8Tbifsys+eZ2dvGnMfbgBJ4tpmtmVnd1/t/j9l+Up5nln0plstycHfdjsiN2Nv5J4H3EweSLwNOTNj+F4GfaN1fA14C3A18GvjRoe3PAo9s3X8S8HFiUuaNwGXp8UBs1XF7es0/EJPplp7vEqekfMGY81onLnR2SzqX9wCPbT3/h8BPte6/jVgR0r49Kj33fOJV1nPp5/XAPVqvfQLxi+LQ//x00023o30jLobyTqAPbAJ/CjxgwvYPJlY/2IRtrkox8AJx3YDHtJ57JLFFSH3fgF9Isfn29Hsdl782xc7zKW7Xt0e29vVRoDPmPN5OnFp4N7FC44lDz8/0/dF6/meBW9Ptt4FLW8+9CHj2Yf956qabbkf7lmLeY2bc9pnAr0/Z5jHA36d4/TbgqtZzPwX8Yfr9nsRijjPE8fw7gK8f2td6eu7RI47z74FXTjiP3XxvvAr4bIrhfz8q9hIT6/9hzLH+FnjIYf9Z6qabbqt9S2PUXwc+k+Ljj6QY+I1TXvcsYhsRT6/9j0CYsP22/MiYbcaOYYk9tl/cuj823s6wr/8BPH/CeXwp8O607/cAX9p67nuAG1v3p+V5Ju1rnZhHueKw/x7oplv9D1c5AlK1xQ+4+5/OuP09iUmHL3X3C4s8t6Hj/jPgme4+boHJA2Nm7wS+393/7rDPRUSOBzN7BfAhd3/eDNv+DvAadx+3eNiBMLOfBm5199845PO4F/B/iN9bm4d5LiIitTSj773EhPIth3wufwz8sLt/4JDP41uA73P37zrM8xCR1WdmDlzr7iNnZE95bU4sHrna3T86ZdtDyY+MOZe/Jn6nfPaQz+OHgPu7+48f5nmIAEpgHyW7TWCLiMjB200CW0RERETkODuoBLaILDf1wBYRERERERERkZVhZj9lZmdH3NSvWeQIUgW2iIiIiIiIiIiIiCwlVWCLiIiIiIiIiIiIyFJSAltEREREREREREREllJ+2CcAcPnll/tVV1112KchIvv07ne/+zZ3v+deXvuNX3fSP3t7ufdjv2/rLe7+2D3vQGaieC1yNCx7vE4LU58BSqBw94ft+YDHlOK1yNGw7PFa9k/xWuRoULxerKVIYF911VXccMMNh30aIrJPZvaxvb72tttL3vmW++352J37fPjyaduY2Y8APwA48LfA09x9c88HPYYUr0WOhmWP18nXufttez7QMad4LXI0rEi8ln1QvBY5GhSvF0stRETkWDCzK4FnAw9z9wcDGfDEwz0rERERERERERGZZCkqsEVEwCm9WvRBcmDDzPrACeCTiz6giMjRcyDx2oE/NjMHfsPdr1/0AUVEjp4DidciIrJvitfTKIEtIkvBgQpf3P7dbzazXwI+DlwA/tjd/3hhBxQROaLmEK8vN7P2XOnrRySov9rdP2lm9wL+xMz+3t3/fD8HFRE5bhY9vhYRkflQvJ5OCWwRWRoV+7riODEhYmaXAtcBVwN3Ar9nZt/r7q/Yz0FFRI6jfcbr26Ytyujun0w/P2NmrwceDiiBLSKyS/uM1yIickAUrydTAltEjoppCZHHAP/k7rcCmNnvA18FKIEtIrJEzOwkENz9TPr9G4CfO+TTEhERERGRQ6IEtogsBccpfaFTZj4OfKWZnSC2EHk0oOW+RUR26QDi9RXA680M4lj1d9z9jxZ5QBGRo+gA4rWIiMyB4vV0SmCLyNJYcA/sd5rZa4H3AAXwXkCLgomI7MGC4/VHgIcs7AAiIseIeqqKiKwGxevJlMAWkaXgQLnggO3uzwWeu9CDiIgccQcRr0VEZP8Ur0VEVoPi9XRKYIvI0tAVRxGR1aB4LSKyGhSvRURWg+L1ZOGwT0BEREREREREREREZBRVYIvIUnDQogUiIitA8VpEZDUoXouIrAbF6+mUwBaRpVEd9gmIiMhMFK9FRFaD4rWIyGpQvJ5MCWwRWQqOa9ECEZEVoHgtIrIaFK9FRFaD4vV0SmCLyHJwKBWvRUSWn+K1iMhqULwWEVkNitdTaRFHEREREREREREREVlKqsAWkaXgqOeTiMgqULwWEVkNitciIqtB8Xo6JbBFZEkYJXbYJyEiIlMpXouIrAbFaxGR1aB4PY0S2CKyFByo1PNJRGTpKV6LiKwGxWsRkdWgeD2demCLiIiIiIiIiIiIyFJSBbaILA1NmRERWQ2K1yIiq0HxWkRkNSheT6YEtogsBUcBW0RkFShei4isBsVrEZHVoHg9nRLYIrI0KlfAFhFZBYrXIiKrQfFaRGQ1KF5PpgS2iCwFXXEUEVkNitciIqtB8VpEZDUoXk+nRRxFREREREREREREZCmpAltEloJjlLqmJiKy9BSvRURWg+K1iMhqULyeTglsEVka6vkkIrIaFK9FRFaD4rWIyGpQvJ5MCWwRWQrq+SQishoUr0VEVoPitYjIalC8nk4JbBFZEkbpmjIjIrL8FK9FRFaD4rWIyGpQvJ5Gn46IiIiIiIiIiIiILCVVYIvIUnCg0jU1EZGlp3gtIrIaFK9FRFaD4vV0SmCLyNJQzycRkdWgeC0ishoUr0VEVoPi9WRKYIvIUnBXzycRkVWgeC0ishoUr0VEVoPi9XT6dERERERERERERERkKakCW0SWRqUpMyIiK0HxWkRkNShei4isBsXryZTAFpGl4ECpSSEiIktP8VpEZDUoXouIrAbF6+mmfjpm9hIz+4yZ/V3rscvM7E/M7B/Tz0tbz/2kmX3IzD5oZt+4qBMXkaMm9nza600Ur0XkoChe75fitYgcDMXr/VK8FpGDoXg9zSzv8qXAY4ceew7wVne/F