{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Use examples for numerical calculations\n", "This jupyter notebook can be found at:\n", "[https://github.com/environmentalscience/essm/blob/master/docs/examples/examples_numerics.ipynb](https://github.com/environmentalscience/essm/blob/master/docs/examples/examples_numerics.ipynb)\n", "\n", "Below, we will import variable and equation defintions that were previously exported from api_features.ipynb\n", "by running the file `test_equation_definitions.py`:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from IPython.display import display\n", "from sympy import init_printing, latex\n", "init_printing() \n", "from sympy.printing import StrPrinter\n", "StrPrinter._print_Quantity = lambda self, expr: str(expr.abbrev) # displays short units (m instead of meter)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/stan/Programs/essm/essm/variables/_core.py:89: UserWarning: \"test_variable_definitions:p_Dva1\" will be overridden by \"__main__:\"\n", " instance[expr] = instance\n", "/home/stan/Programs/essm/essm/variables/_core.py:89: UserWarning: \"test_variable_definitions:p_Dva2\" will be overridden by \"__main__:\"\n", " instance[expr] = instance\n", "/home/stan/Programs/essm/essm/variables/_core.py:89: UserWarning: \"test_variable_definitions:p_alpha2\" will be overridden by \"__main__:\"\n", " instance[expr] = instance\n", "/home/stan/Programs/essm/essm/variables/_core.py:89: UserWarning: \"test_variable_definitions:p_alpha1\" will be overridden by \"__main__:\"\n", " instance[expr] = instance\n", "/home/stan/Programs/essm/essm/variables/_core.py:89: UserWarning: \"test_variable_definitions:p_ka2\" will be overridden by \"__main__:\"\n", " instance[expr] = instance\n", "/home/stan/Programs/essm/essm/variables/_core.py:89: UserWarning: \"test_variable_definitions:p_ka1\" will be overridden by \"__main__:\"\n", " instance[expr] = instance\n", "/home/stan/Programs/essm/essm/variables/_core.py:89: UserWarning: \"test_variable_definitions:p_nua1\" will be overridden by \"__main__:\"\n", " instance[expr] = instance\n", "/home/stan/Programs/essm/essm/variables/_core.py:89: UserWarning: \"test_variable_definitions:p_nua2\" will be overridden by \"__main__:\"\n", " instance[expr] = instance\n", "/home/stan/Programs/essm/essm/variables/_core.py:89: UserWarning: \"test_variable_definitions:p_CC2\" will be overridden by \"__main__:\"\n", " instance[expr] = instance\n", "/home/stan/Programs/essm/essm/variables/_core.py:89: UserWarning: \"test_variable_definitions:p_CC1\" will be overridden by \"__main__:\"\n", " instance[expr] = instance\n" ] } ], "source": [ "%run -i 'test_equation_definitions.py'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Numerical evaluations\n", "See here for detailed instructions on how to turn sympy expressions into code: https://docs.sympy.org/latest/modules/codegen.html\n", "\n", "We will first list all equations defined in this worksheet:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "eq_Le: Eq(Le, alpha_a/D_va)\n", "eq_Cwa: Eq(C_wa, P_wa/(R_mol*T_a))\n", "eq_Nu_forced_all: Eq(Nu, -Pr**(1/3)*(-37*Re**(4/5) + 37*(Re + Re_c - Abs(Re - Re_c)/2)**(4/5) - 664*sqrt(Re + Re_c - Abs(Re - Re_c)/2))/1000)\n", "eq_Dva: Eq(D_va, T_a*p_Dva1 - p_Dva2)\n", "eq_alphaa: Eq(alpha_a, T_a*p_alpha1 - p_alpha2)\n", "eq_ka: Eq(k_a, T_a*p_ka1 + p_ka2)\n", "eq_nua: Eq(nu_a, T_a*p_nua1 - p_nua2)\n", "eq_rhoa_Pwa_Ta: Eq(rho_a, (M_N2*P_N2 + M_O2*P_O2 + M_w*P_wa)/(R_mol*T_a))\n", "eq_Pa: Eq(P_a, P_N2 + P_O2 + P_wa)\n", "eq_PN2_PO2: Eq(P_N2, P_O2*x_N2/x_O2)\n", "eq_ideal_gas_law: Eq(P_g*V_g, R_mol*T_g*n_g)\n", "eq_Pwa_CC: Eq(P_wa, p_CC1*exp(-M_w*lambda_E*(-1/p_CC2 + 1/T_g)/R_mol))\n", "eq1: Eq(P_wa, Piecewise((0, T_a < 0), (p_CC1*exp(-M_w*lambda_E*(-1/p_CC2 + 1/T_g)/R_mol), True)))\n", "eq_Pwa_Delta: Eq(P_wa, P_wa1 + Integral(Delta_Pwa, (T_g, T_a1, T_a2)))\n", "eq_PO2: Eq(P_O2, (P_a*x_O2 - P_wa*x_O2)/(x_N2 + x_O2))\n", "eq_PN2: Eq(P_N2, (P_a*x_N2 - P_wa*x_N2)/(x_N2 + x_O2))\n", "eq_rhoa: Eq(rho_a, (x_N2*(M_N2*P_a - P_wa*(M_N2 - M_w)) + x_O2*(M_O2*P_a - P_wa*(M_O2 - M_w)))/(R_mol*T_a*x_N2 + R_mol*T_a*x_O2))\n", "eq_Pg: Eq(P_g, R_mol*T_g*n_g/V_g)\n", "eq_Pwa_nw: Eq(P_wa, R_mol*T_g*n_w/V_g)\n" ] } ], "source": [ "for eq in Equation.__registry__.keys():\n", " print(eq.definition.name + ': ' + str(eq))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Substitution of equations and values into equations\n", "The easiest way is to define a dictionary with all variables we want to substitute as keys. We start with the default variables and then add more. First, however, we will define a function to display the contents of a dictionary:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def print_dict(vdict, list_vars=None):\n", " \"\"\"Print values and units of variables in vdict.\"\"\"\n", " if not list_vars:\n", " list_vars = vdict.keys()\n", " for var1 in list_vars:\n", " unit1 = var1.definition.unit\n", " if unit1 == 1:\n", " unit1 = ''\n", " if vdict[var1] is not None:\n", " print('{0}: {1} {2}'.format(var1.name, str(vdict[var1]), str(unit1)))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "c_pa: 1010.0 J/(K*kg)\n", "c_pamol: 29.19 J/(K*mol)\n", "c_pv: 1864 J/(K*kg)\n", "g: 9.81 m/s**2\n", "lambda_E: 2450000.0 J/kg\n", "M_air: 0.02897 kg/mol\n", "M_N2: 0.028 kg/mol\n", "M_O2: 0.032 kg/mol\n", "M_w: 0.018 kg/mol\n", "R_mol: 8.314472 J/(K*mol)\n", "sigm: 5.67e-08 J/(K**4*m**2*s)\n", "T0: 273.15 K\n", "x_N2: 0.79 \n", "x_O2: 0.21 \n", "p_Dva1: 1.49e-07 m**2/(K*s)\n", "p_Dva2: 1.96e-05 m**2/s\n", "p_alpha1: 1.32e-07 m**2/(K*s)\n", "p_alpha2: 1.73e-05 m**2/s\n", "p_ka1: 6.84e-05 J/(K**2*m*s)\n", "p_ka2: 0.00563 J/(K*m*s)\n", "p_nua1: 9e-08 m**2/(K*s)\n", "p_nua2: 1.13e-05 m**2/s\n", "p_CC1: 611.0 Pa\n", "p_CC2: 273.0 K\n" ] } ], "source": [ "vdict = Variable.__defaults__.copy()\n", "print_dict(vdict)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can substitute a range of equations into each other by using the custom function `subs_eq`:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJwAAAAuBAMAAAAvlG0AAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdrtEVN3vqxDNIomZZjLe39VDAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACmElEQVRIDe1WPYgTQRT+ks0lm2x2k9hcY5FKwSYRrhPJemhhdRGtBM01opUXC9Er5AYtFYJVFA5NYSleOAjc2bgKovgDV4mFkoggWMgdyaHxCMS3kz0vm012Jyt2ebDzZr73vS8zLzPwAE+L/Ty41Lz2sOpJFCMoVYTziG+IsT1ZEWDFgJr1JIoRfgB1HQEmxvZkNYAlBsmTJ064YFGl1Sep6860I4e+TDMnPAqJdayIEjmDVWjpAWL1qq6klRevB+BRS3XbijyKFlHX3uXtxJhehrr1FIt2eOQq0rZCLKHjGZO4nHyfWwOQsQ2lfRablZEKtgBdO8tyDAfQk9uFyEttBOdPY0XwLtG1A+5MPwBqkDtOOSWPElHmGNEEjK4d5PMaHXkZ6iWnXNTABwacEpAyKRnixg2NjtzEt4ZTLpHWLgOKYXI9TVnslpMI64EipFbqDZXqbymt3Nzh/RXgU2+1r5tFqHCRWbERrlQJfmRUJDKHXI2jValhemTSQJSZMxcLVzaPgopEFhrc3bKJfj8+y0wvze7Qwc2Zm0lXjr1HSSdK7OWvpI0p81u+0O1yVLlFb2jGRhixCBXuDouc+N0vH0cBWBvG84XFUd+A4St1WNIMEnlUgZNNGixb5++Rl3gXEvVrUFpSFlDnRTNceQZiTfOPjaSH07pjGCnQETOPyYeTw+XGQyUdyL2lnBIdGLdTXLRXu/J4SpytMGCqRdM5+uR7cpGjvocpylS3aFinL3Ljs28hnnhz4Rx5c0t0Lb6Gs/+m1pdNL6+YqID1QXzqr0fRdhAoBnVn7fz1KLXOq0JWfp4a3BwmPYqtJJMexVYO22LSo/SV43/2KH0/4zIV71FcRPZCwj3KXorbTLhHcROxx7x7lD8fURZXwH0iVgAAAABJRU5ErkJggg==\n", "text/latex": [ "$$N_{Le} = \\frac{T_a p_1 - p_2}{T_a p_1 - p_2}$$" ], "text/plain": [ " Tₐ⋅pₐₗₚₕₐ₁ - pₐₗₚₕₐ₂\n", "Le = ────────────────────\n", " Tₐ⋅p_Dva1 - p_Dva2 " ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from essm.variables.utils import subs_eq\n", "subs_eq(eq_Le, [eq_alphaa, eq_Dva])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also use subs_eq to substitute equations into each other and a dictionary with values. We will first add an entry for T_a into the dictionary and then substitute:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOwAAAARBAMAAAAs8V3gAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdrtEVN3vqxDNIomZZjLe39VDAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADtUlEQVQ4Eb1VXWhURxT+9vfevfuXlKJghW6tpT4Ucx/6UhBy6YM/2Hb3rbRSdlFsUrG6D22hTTRTxYooZN9WQdt9aR9aC5uCqKzFtQoRspIVSmuVkAUfFRJsajapevuduZrcmO2jHphv58x8Z745Z+bOAnjBtRHO9yn8j43casmM2X3B9gGQqsvwXbaRl1owG2ta4gOxG6krXg8/1dYCx2pfiavpAsd+W6+nezLkKt3tADHH/FCGE8BnPgAiRRl+HzhiJ0o4jPjfQLgfiLpulRMbFaw6sk3zU/yg6Gs6wSziok0/9PY8EGens/0BnJeZd7m+D4AVIhvoA95AMIOrwMc4fLbNDM6uYm0GhhQCM0jnjAKMOpmaTkgoRAv04wcfAG+y09leB7KKU78Ah3wA6y+R/bmMEHcN7ASGKkhS1hAXmFQI7UC6lC4hPg2PLjHRHAyJSCAPnNPUTjAL9Dqc6F2HX32AeExkq2UEC/zFdwpDradkOZxtDmcQoIymC8TaMB5xJoHJJursdDSTpejt4pTh3nJ8gJsim6yUkf525LIOzStPdmJ1iz6zZa23IcpseeqaroHHzJqwvukcqsCWe4THdv6E2El61r/AlMji4pzygdklsl+jjOFTiDicD5EpRQ4qizJa1jpzVSocewhN92KYhgid425CNkMKdJaZxWynhBa+nN/hgzhEtiqy00jwAuurKrK019h0tlhh4xVsneGVJd2L0bcfrK91Ty6ykSEss4Ui/4nAfbUI78k64RZlWaukpNdgeyw7XnkiG+1H+PqdR9B0DTzVOplS357T/I12ibfMeKWyDkeZ0JSzAJurInsUlI2UkHzIT6lOksh+yaOwtaxVQXCGo0bbFLoGuqvZEHJY7DF2hsnF0W4t7p1tmQNY423cnGb4jQX4fbAxvrtrU6Nx/xI/TMn2HVBFZD8BxpUOSrc92VgxLPQXdQzCdfzIm01KRIqUZTOP8xFZak+ei128LJVFYHU0sw+pNhIFpKoIeLJ14AMuwbONZnhpk/2YqtDXdIE7wCpKcijJLPRbZOybYG+JRRxzD4ZLeEvhug+4TkZ424EBfOPg+9qFz1lpZnsaSe4OLyskm/xuA8fNU0LUdIK1rTaaw4G9H3FMNs7P5XbUFsYSM0dWtuSWpgb5V7AIiPfMNvl2ufsQH9vP58J1/0Hii7lrsMYGFTC692QTGxqvArWGLKrpAjHXdXM+BR5+MV0BQ56npeYRKAadZWf7jPdw5sFo3jYvdT9jGW/5/wAdKVDfo2JeQAAAAABJRU5ErkJggg==\n", "text/latex": [ "$$N_{Le} = 0.888446215139442$$" ], "text/plain": [ "Le = 0.888446215139442" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "vdict[T_a] = 300.\n", "subs_eq(eq_Le, [eq_alphaa, eq_Dva], vdict)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Evaluation of equations for long lists of variable sets\n", "Substitution of variables into equations takes a lot of time if they need to be evaluated for a large number of variables. We can use theano to speed this up:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] } ], "source": [ "#import theano\n", "from sympy.printing.theanocode import theano_function\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now create two long lists of values representing T_g and n_g respectively and show how long it takes to compute ideal gas law values." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "npoints = 10000\n", "xmin = 290.\n", "xmax = 310.\n", "Tvals = np.arange(xmin, xmax, (xmax - xmin)/npoints)\n", "xmin = 0.1\n", "xmax = 0.5\n", "nvals = np.arange(xmin, xmax, (xmax-xmin)/npoints)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 7.46 s, sys: 43.4 ms, total: 7.5 s\n", "Wall time: 7.68 s\n" ] } ], "source": [ "%%time\n", "# looping\n", "expr = eq_ideal_gas_law.rhs.subs(Variable.__defaults__)\n", "resvals0 = []\n", "for i in range(len(Tvals)):\n", " resvals0.append(expr.subs({T_g: Tvals[i], n_g: nvals[i]}))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 106 ms, sys: 11.7 ms, total: 118 ms\n", "Wall time: 623 ms\n" ] } ], "source": [ "%%time\n", "# Using theano\n", "f1 = theano_function([T_g, n_g], [eq_ideal_gas_law.rhs.subs(Variable.__defaults__)], dims={T_g:1, n_g:1})\n", "resvals1 = f1(Tvals,nvals)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(resvals0) == list(resvals1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Both approaches give identical results, but `theano_function` makes it a lot faster.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Numerical solution\n", "Some equations cannot be solved analytically for a given variable, e.g. eq_Nu_forced_all cannot be solved analytically for Re if Nu is given, so we can use numerical solvers instead:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "from sympy import nsolve" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIwAAAAPBAMAAADEyjp7AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAiXYyEM1Embsi72bdVKu+2mc6AAAACXBIWXMAAA7EAAAOxAGVKw4bAAACeklEQVQ4Ea2TvWtTURjGfze5TfPZJO3QDiK1i1AR08mxqSBCKSaguBTbbIIgRlGvdDE4aBGkt5P4Ac1QqtQPguggDuogTtaIg5tecRA/oK2mhVpjfN9z0v/AOzw553nO+Z2vN0D3jhwHsl/B2dUDrjedk+ZAD7FsXx39jD+bvartvdqXEMLlLVEA+0udRfcuN/IcZcznIM467hHOco3UpozE+LEyH2vSXpKFNIRk0BYFwGWila4qXeVEFafIGVgkWmaaPXX+yEiMH2qSfg6hB9gQBgPJVBRAZEN66SLOSjQg1eAeLPsnfTGPl1iTH4wfWdBB3JwAExJ7HGDEAIhWZWShQmgjHRBpMplnOXdFPPnsoayPOVRGMDZ0OgK5JREDIH2o+xxx2c16oUzkp85+k18b2lmSxlhGu23ffSHn8wVjwyeKUTEACuMk684KHZtzGcKKiay6azVew+yrvHSxfmyb3No+JrCh2y8EIwZAYYXOeS7ytjlXJrwq05KB24IpQYQWDKbtD5bIKMaEDoIxYgFy/4l1Ug+Hfm0dykNf6FlJEIvC2joU8flUTjA2fKcYIxaQLJLQ8uhqyFWmmrKFMrwUTO06zNQlMX7MJ9ocQTEaPsoIxlXBAqRkZDeIIa+faMB7Yv433U0rbzHGTzcE88Hzfp9GwxNfvAu3d6v0W0C4QWc1Mc+wb8svnCHkD+eZ4o6+mSxg/HiFpKyBlJ8JIR5IV8QA4BOj9dBhdxx5nNEct7K9P+jKuPd5qnXjSgWqn6jZP4OwTSiEimJEFCAldOozZL0SjAxIa7LVkgvu3Z4j4p0Xc6ntH/MuybS+1ncb4sz8rbVFAf/n+wfURvBtNAe4UwAAAABJRU5ErkJggg==\n", "text/latex": [ "$$690263.0346446$$" ], "text/plain": [ "690263.034644600" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "vdict = Variable.__defaults__.copy()\n", "vdict[Pr] = 0.71\n", "vdict[Re_c] = 3000.\n", "vdict[Nu] = 1000.\n", "expr = eq_Nu_forced_all.subs(vdict)\n", "nsolve(expr, 1000.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now applying to a long list of Nu-values:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "npoints = 100\n", "xmin = 1000.\n", "xmax = 1200.\n", "Nuvals = np.arange(xmin, xmax, (xmax - xmin)/npoints)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.56 s, sys: 3.9 ms, total: 1.56 s\n", "Wall time: 1.56 s\n" ] } ], "source": [ "%%time\n", "# Solving for a range of Nu values\n", "vdict = Variable.__defaults__.copy()\n", "vdict[Pr] = 0.71\n", "vdict[Re_c] = 3000.\n", "resvals = []\n", "for Nu1 in Nuvals:\n", " vdict[Nu] = Nu1\n", " resvals.append(nsolve(eq_Nu_forced_all.subs(vdict), 1000.))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now again use a theano function to make it faster. First we import optimize from scipy and preapre the theano_function:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "import scipy.optimize as sciopt\n", "vdict = Variable.__defaults__.copy()\n", "vdict[Pr] = 0.71\n", "vdict[Re_c] = 3000.\n", "expr = eq_Nu_forced_all.subs(vdict)\n", "expr1 = expr.rhs - expr.lhs\n", "fun_tf = theano_function([Re, Nu], [expr1], dims={Nu:1, Re:1})\n", "x0vals = np.full(Nuvals.shape, fill_value=2000.) # array of same shape as Nuvals, with initial guess" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 7.35 ms, sys: 0 ns, total: 7.35 ms\n", "Wall time: 7.24 ms\n" ] } ], "source": [ "%%time\n", "# Solving for a range of Nu values\n", "resvals1 = sciopt.fsolve(fun_tf, args=Nuvals, x0=x0vals)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAN8AAAAVBAMAAAA0pCbNAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAzWYQMplU74mrdiK7RN1/7zyFAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADSElEQVRIDeVUO0xTURj+6G1LSx804KCDocEB4+JhIjhIFxOdqGKFQeyNCTMYjAy+aiKGNiYgkWAcpCaSoA7chcTEREmMMUaC3VyM1EEXBwXig/Co/6O3wuBaB//k/Pf7vvP957/n9NwC1QrvPMCjWuFNFcGjehGhZjyqFv9Pw5p93lN8rKGRmzaUBB5dh9WTdcr00Ewv0GCO2OqxUgsJRUBuJqFJRHDhzsjtJ4OGe6S1pVKGlTn4VyDE2ovjaIuHBpUiiaCNq6VfZU87ApvqRmvBF9MkohTq6pT9g0CNbX3B+wUKIxeGm9ZkH4tn1saWknASC5gEJpT647Bi2H3XAcTTBQwpwkuE5zWJKIVuw/bsKtACXHMFd4fBsvCugGUIWaKlCWNUaWAeVhFjbBPPFWDKERRZJ00SRJRC9kl4qOFroN+Uub9Iu6bhNiS2ouQVOawN2o8jc5ENE4hrQ/E0GUwlBIUHyCkJInLhn+CGa8CSrVLozFqeB4LN0wmVDmeULB+dLuAc7bCgc1MrN4DxAyfZRR6Kr/TWhOp35Xo1qSiFDCWoofWdGubL3H2ETWhT8JNPBkys5Tg+03miKa5zgVIj8ALPCwB76BTp/jDqaEDUliSiFrrrghqGyNiZrygVsKio7jI/F60S0G3Cef+ozRTo+/CTeiFcpCSe6Jiijm/wDUoiHh3TQoIWDUhD2mFnRtiOlHaUDvHrp50t4GABZx+MJoR6xjA7SMjHVwTsSTEgVF+EZ1MSURK1EOiI8fzfjvQW/1g0/xR4awv5SA3jpEwZobUG3h/BGH994kFdsuyOxuBZkURbJ9EtlBPghnxp+m1aanvQ7UgbEkqGGgqhT5B2COzhm5M2HQTHa2PwraoHfQg54g4O0A4l0bmT6Bbq8txwgZbixbeHvhlwiS+fkE6DbrQY7zr9x9Bb1zrAGx/9SDH1eDOoc8TtXYVvQBJElMLK2txw+4fvTszBc5FPvZm/QyHBjHUB901rXGmkEZ58KIkWRz3PRobvKcIdtNmaRJRCd13UUcOobZ2vCGUQ6pkwdMUQSZ0uQAiGswkEhukfW+mxHoJz2S76HtjTVCptKUKg5zY0iSiF7vq+yeUTsHIPE67wr56/AbImTGJbOXpVAAAAAElFTkSuQmCC\n", "text/latex": [ "$$5.35695853236626 \\cdot 10^{-11}$$" ], "text/plain": [ "5.35695853236626e-11" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.mean(abs((resvals - resvals1)/resvals))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Using theano and scipy makes it 2 orders of magnitude faster and the results are different only by 10$^{-10}$%!**\n", "**Note, however, that scipy gets slowed down for large arrays, so it is more efficient to re-run it repreatedly with subsections of the arra:**" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "npoints = 1000\n", "xmin = 1000.\n", "xmax = 1200.\n", "Nuvals = np.arange(xmin, xmax, (xmax - xmin)/npoints)\n", "x0vals = np.full(Nuvals.shape, fill_value=2000.)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.6 s, sys: 4.13 ms, total: 1.61 s\n", "Wall time: 1.61 s\n" ] } ], "source": [ "%%time\n", "# Solving for a range of Nu values\n", "resvals1 = sciopt.fsolve(fun_tf, args=Nuvals, x0=x0vals)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now test that we can process Nuvals bit by bit and re-create it consistently:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True\n" ] } ], "source": [ "# Solving for a range of Nu values\n", "imax = len(Nuvals)\n", "i0 = 0\n", "idiff = 100\n", "i1 = i0\n", "resvals2 = []\n", "while i1 < imax - 1:\n", " i0 = i1 # note that resvals[0:2] + resvals[2:4] = resvals[0:4]\n", " i1 = min(i0+idiff, imax)\n", " resvals0 = Nuvals[i0:i1]\n", " resvals2 = np.append(resvals2,resvals0)\n", "print(list(resvals2) == list(Nuvals))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will run fsolve for portions of Nuvals bit by bit:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 56.5 ms, sys: 16 µs, total: 56.5 ms\n", "Wall time: 56 ms\n" ] } ], "source": [ "%%time\n", "# Solving for a range of Nu values\n", "imax = len(Nuvals)\n", "i0 = 0\n", "idiff = 100\n", "i1 = i0\n", "resvals2 = []\n", "while i1 < imax - 1:\n", " i0 = i1 # note that resvals[0:2] + resvals[2:4] = resvals[0:4]\n", " i1 = min(i0+idiff, imax)\n", " resvals0 = sciopt.fsolve(fun_tf, args=Nuvals[i0:i1], x0=x0vals[i0:i1])\n", " resvals2 = np.append(resvals2,resvals0)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOIAAAAPBAMAAAAL5A64AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMt2rmYlmIkR2uxDNVO+L8+I6AAAACXBIWXMAAA7EAAAOxAGVKw4bAAADaklEQVQ4Eb2UTWhcVRTHf/P15uNlPih0HTfuKh0lxY1tBpxqodo8TUaIGwfESqmQR9FtO4i7Cs0iKtnUobjQivpo4y52piBFUOkI6tIEkbYw1aYVa8okPs8597UMpG49MOe9e89/fv/7zr3vsWuq1QpIwm+er2GJ3PILkHtf0lutz6Q8/55Uluc6TBxPPZ9o7F8rrSfBKFaBVCOh7bwYbyGO47rWykO4SG7TJf8bnoPvuRDmG6y1We1k+2RhP5U4DpwGfonwr3EywihWgbThdtiJ1ng8CxmtXpwdwWH4w6VSg/N46+T6xW2qj/ARpS4/wyoTsx8kQn/fdERmnUzDUawCZx/maFrjEYB0TsITx19hOrR0JpSpUp3yqHCDar/wt0pehgtufaaBhUhq5DYcxVZO/o2HOZrWeMIpLCrNHK9GTNcsfa5T1TqFbbmutUvrOj7zKS85R9MYpdelqMsRinPMTfyn431eSWHOUS57I0tbu5c69BoU7oH/KNXL84ek/fFrAzJHl2qqV+FCREWecVOGQnGV19XRn18KVTQWok14smwL7aos9S9L/labxzgVkLpH/sXD9K6QHsArWxGlKK8OJhSKtHTiroyFYhV/UR3fDvOBTI6HaB0Pf8PNO0c7aOm6H8OR6FSDlK7gbKe3QXZI+dDeGyo+Jz8TCoVP+Ek6n1DOkUMc/d9ZqalyLNTR8Yp9N+0cmzqQdEda1rnfhcpQTqu3yasU/xEPZkLTWFcpH9t9GxLKTPiDOmauLR9X0Op3Gpf1Vncg4aW7NmEnh2JDBpr2iKQtO13ezoeUttN9vLsMYXLwlDS344RKkcjIfihFK6cDdZQ/7gjRGs8OpFXtGX9EDDQd0GeUV9AbVUfiKG+dt6mdy7RvyjNGpnHPiFqYh1a+mmrO/LZYDXYY6uqMJ6tuuKo6pgKKoaXJiCPuC1Dpkh6lRmTXuS6nIxT5Hic0R2/IZOgoVoFKnbQ4LjrsgyyO7gvAmlR7femlOJ5ozT3hUiaQ3ZfjeqnmtfUrt49LA96NOMabeNedBj6MKH7rXxGqUqwijl3d1ZXwgZe7Ea3x4OTADl72wNYzXI3jOy4xN1uDL5cfh3eaH8s3/aDcpabmOuQPTkWJ5ulbX7dpNTuCVIpVyH3xZ5ulo22ZHA/TGm989n+5/xc1gERnYT01lAAAAABJRU5ErkJggg==\n", "text/latex": [ "$$7.122603685219754e-10$$" ], "text/plain": [ "7.122603685219754e-10" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.mean(abs((resvals1 - resvals2)/resvals1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**It is strange that resvals1 and resvals2 are different at all, but anyway, it is clear that slicing the data in relatively small portions is important to keep `scipy.optimize.fsolve` time-efficient.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generate code from sympy expressions and execute\n", "Need to install gfortran system-wide first!" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "from sympy.utilities.autowrap import autowrap" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "from sympy import symbols\n", "x, y, z = symbols('x y z')\n", "expr = ((x - y + z)**(13)).expand()\n", "autowrap_func = autowrap(expr)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 8 µs, sys: 0 ns, total: 8 µs\n", "Wall time: 13.4 µs\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAACkAAAAPBAMAAACLu/vuAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM0yVO+riWZ2md0iu0S3uypJAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAg0lEQVQYGWNgQAeV7QZIQiyBYA6TA+MehKhY2Ucwp5WBoQghysABEb3GwOAvgBCGin5jYHjvgC7K+AsoOgFdlOc7A4M9pihQrX0Bg5ASCCgzQG3DbgID0DZ/DNsY7jIw9GO6DM0XnEBfyCswcDswBiMcxpX6aQUDdwIDY+U0A4QoCgsAmwIj8jixtwMAAAAASUVORK5CYII=\n", "text/latex": [ "$$-1.0$$" ], "text/plain": [ "-1.0" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "autowrap_func(1, 4, 2)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 70 ms, sys: 3.77 ms, total: 73.8 ms\n", "Wall time: 72.4 ms\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAABgAAAAPBAMAAAAMihLoAAAAJFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADHJj5lAAAAC3RSTlMAEM0yVO+riWZ2mSMU5s8AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAAsSURBVAgdY2DAClgCEcJiZRsRHAYOMjlCSiCgzECWAZxIlnKlblqB5B4GBgDVtwt2YFScIgAAAABJRU5ErkJggg==\n", "text/latex": [ "$$-1$$" ], "text/plain": [ "-1" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "expr.subs({x:1, y:4, z:2})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use of `autowrap` made the calculation 3 orders of magnitude faster than substitution of values into the original expression!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another way is to use `binary_function`:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "from sympy.utilities.autowrap import binary_function\n", "f = binary_function('f', expr)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 2.47 ms, sys: 0 ns, total: 2.47 ms\n", "Wall time: 2.65 ms\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAACkAAAAPBAMAAACLu/vuAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM0yVO+riWZ2md0iu0S3uypJAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAg0lEQVQYGWNgQAeV7QZIQiyBYA6TA+MehKhY2Ucwp5WBoQghysABEb3GwOAvgBCGin5jYHjvgC7K+AsoOgFdlOc7A4M9pihQrX0Bg5ASCCgzQG3DbgID0DZ/DNsY7jIw9GO6DM0XnEBfyCswcDswBiMcxpX6aQUDdwIDY+U0A4QoCgsAmwIj8jixtwMAAAAASUVORK5CYII=\n", "text/latex": [ "$$-1.0$$" ], "text/plain": [ "-1.0" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "f(x,y,z).evalf(2, subs={x:1, y:4, z:2})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, for the above example, `binary_function` was 2 orders of magnitude slower than `autowrap`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.4" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }