This jupyter notebook can be found at: https://github.com/environmentalscience/essm/blob/master/docs/examples/api_features.ipynb
Below, we will import some generic python packages that are used in this notebook:
In [1]:
# Checking for essm version installed
import pkg_resources
pkg_resources.get_distribution("essm").version
Out[1]:
In [2]:
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:120% !important; }</style>"))
In [3]:
from IPython.display import display
from sympy import init_printing, latex
init_printing()
from sympy.printing import StrPrinter
StrPrinter._print_Quantity = lambda self, expr: str(expr.abbrev) # displays short units (m instead of meter)
In [4]:
import scipy as sc
# Import various functions from sympy
from sympy import Derivative, Eq, exp, log, solve, Symbol
from essm.equations import Equation
from essm.variables import Variable
from essm.variables.utils import ListTable, generate_metadata_table
In [5]:
from essm.variables.chamber import *
from essm.variables.leaf import *
from essm.variables.physics.thermodynamics import *
from essm.equations.leaf import *
from essm.equations.physics.thermodynamics import *
In [6]:
import matplotlib.pyplot as plt
from sympy import latex
from sympy import N
from numpy import arange
from essm.variables.units import derive_unit, SI, Quantity
from essm.variables.utils import markdown
def plot_expr2(xvar_min_max, yldata, yllabel=None, yrdata=None,
yrlabel='', clf=True, npoints=100, ylmin=None, ylmax=None,
yrmin=None, yrmax=None, xlabel=None,
colors=None,
loc_legend_left='best', loc_legend_right='right',
linestylesl=['-', '--', '-.', ':'],
linestylesr=['-', '--', '-.', ':'],
fontsize=None, fontsize_ticks=None, fontsize_labels=None,
fontsize_legend=None,
fig1=None, **args):
'''
Plot expressions as function of xvar from xmin to xmax.
**Examples:**
from essm.variables import Variable
from essm.variables.physics.thermodynamics import T_a
from essm.equations.physics.thermodynamics import eq_nua, eq_ka
vdict = Variable.__defaults__.copy()
expr = eq_nua.subs(vdict)
exprr = eq_ka.subs(vdict)
xvar = T_a
yldata = [(expr.rhs, 'full'), (expr.rhs/2, 'half')]
yrdata = exprr
plot_expr2((T_a, 273, 373), yldata, yllabel = (nu_a), yrdata=yrdata)
plot_expr2((T_a, 273, 373), yldata, yllabel = (nu_a),
yrdata=[(1/exprr.lhs, 1/exprr.rhs)],
loc_legend_right='lower right')
plot_expr2((T_a, 273, 373), expr)
plot_expr2((T_a, 273, 373), yldata, yllabel = (nu_a))
'''
(xvar, xmin, xmax) = xvar_min_max
if not colors:
if yrdata is not None:
colors = ['black', 'blue', 'red', 'green']
else:
colors = ['blue', 'black', 'red', 'green']
if fontsize:
fontsize_labels = fontsize
fontsize_legend = fontsize
fontsize_ticks = fontsize
if not fig1:
plt.close
plt.clf
fig = plt.figure(**args)
else:
fig = fig1
if hasattr(xvar, 'definition'):
unit1 = derive_unit(xvar)
if unit1 != 1:
strunit = ' (' + markdown(unit1) + ')'
else:
strunit = ''
if not xlabel:
xlabel = '$'+latex(xvar)+'$'+ strunit
else:
if not xlabel:
xlabel = xvar
if hasattr(yldata, 'lhs'):
yldata = (yldata.rhs, yldata.lhs)
if not yllabel:
if type(yldata) is tuple:
yllabel = yldata[1]
else:
try:
yllabel = yldata[0][1]
except Exception as e1:
print(e1)
print('yldata must be equation or list of (expr, name) tuples')
if type(yllabel) is not str:
unit1 = derive_unit(yllabel)
if unit1 != 1:
strunit = ' (' + markdown(unit1) + ')'
else:
strunit = ''
yllabel = '$'+latex(yllabel)+'$'+ strunit
if type (yldata) is not list and type(yldata) is not tuple:
# If only an expression given
yldata = [(yldata, '')]
if type(yldata[0]) is not tuple:
yldata = [yldata]
if yrdata is not None:
if yrlabel == '':
if hasattr(yrdata, 'lhs'):
yrlabel = yrdata.lhs
if type (yrdata) is not list and type(yrdata) is not tuple:
# If only an expression given
yrdata = [yrdata]
if type(yrlabel) is not str:
yrlabel = '$'+latex(yrlabel)+'$'+ ' (' + markdown(derive_unit(yrlabel)) + ')'
xstep = (xmax - xmin)/npoints
xvals = arange(xmin, xmax, xstep)
ax1 = fig.add_subplot(1, 1, 1)
if yrdata is not None:
color = colors[0]
else:
color = 'black'
if ylmin: ax1.set_ylim(ymin=float(ylmin))
if ylmax: ax1.set_ylim(ymax=float(ylmax))
ax1.set_xlabel(xlabel)
ax1.set_ylabel(yllabel, color=color)
ax1.tick_params(axis='y', labelcolor=color)
i = 0
for (expr1, y1var) in yldata:
linestyle = linestylesl[i]
if yrdata is None:
color = colors[i]
i= i + 1
try:
y1vals = [expr1.subs(xvar, dummy).n() for dummy in xvals]
ax1.plot(xvals, y1vals, color=color, linestyle=linestyle, label=y1var)
except Exception as e1:
print([expr1.subs(xvar, dummy) for dummy in xvals])
print(e1)
if i > 1 or yrdata is not None:
plt.legend(loc=loc_legend_left, fontsize=fontsize_legend)
if yrdata is not None:
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
color = colors[1]
ax2.set_ylabel(yrlabel, color=color)
i = 0
for item in yrdata:
if type(item) is tuple: # if item is tuple
(expr2, y2var) = item
else:
try:
(y2var, expr2) = (item.lhs, item.rhs)
except Exception as e1:
print(e1)
print('yrdata must be a list of equations or tuples (var, expr)')
return
linestyle = linestylesr[i]
i = i + 1
try:
y2vals = [expr2.subs(xvar, dummy).n() for dummy in xvals]
ax2.plot(xvals, y2vals, color=color, linestyle=linestyle, label=y2var)
except Exception as e1:
print(expr2)
print([expr2.subs(xvar, dummy).n() for dummy in xvals])
print(e1)
if not yrlabel:
if hasattr(yrdata[0], 'lhs'):
yrlabel = yrdata[0].lhs
if type(yrlabel) is not str:
yrlabel = '$'+latex(yrlabel)+'$'+ ' (' + markdown(derive_unit(yrlabel)) + ')'
ax2.tick_params(axis='y', labelcolor=color)
if yrmin: ax2.set_ylim(ymin=float(yrmin))
if yrmax: ax2.set_ylim(ymax=float(yrmax))
leg=ax2.legend(loc=loc_legend_right, fontsize=fontsize_legend)
ax2.add_artist(leg);
for item in ([ax2.xaxis.label, ax2.yaxis.label]):
item.set_fontsize(fontsize_labels)
ax2.tick_params(axis='both', which='major', labelsize=fontsize_ticks)
for item in ([ax1.xaxis.label, ax1.yaxis.label]):
item.set_fontsize(fontsize_labels)
ax1.tick_params(axis='both', which='major', labelsize=fontsize_ticks)
fig.tight_layout() # otherwise the right y-label is slightly clipped
return fig
In [7]:
vdict = Variable.__defaults__.copy()
expr = eq_nua.subs(vdict)
exprr = eq_ka.subs(vdict)
xvar = T_a
yldata = [(expr.rhs, 'full'), (expr.rhs/2, 'half')]
yrdata = exprr
fig=plot_expr2((T_a, 273, 373), yldata=expr, yrdata=exprr, yrmin=-0.0001, fontsize=14) # note that yrmin=0 would have no effect
fig=plot_expr2((T_a, 273, 373), yldata=expr, yrdata=exprr, colors=['red', 'blue'], linestylesr=['--'])
fig=plot_expr2((T_a, 273, 373), yldata, yllabel = (nu_a), yrdata=yrdata)
fig=plot_expr2((T_a, 273, 373), yldata, yllabel = (nu_a), yrdata=[(1/exprr.rhs, 1/exprr.lhs)],
loc_legend_right='lower right')
fig=plot_expr2((T_a, 273, 373), expr)
fig=plot_expr2((T_a, 273, 373), yldata, yllabel = (nu_a))
We can also manipulate the figure later, e.g. change the width:
In [8]:
%matplotlib inline
fig.set_figwidth(8)
fig
Out[8]:
In [9]:
from essm.variables import Variable
To define units, you can either import these units from the library, e.g.
from essm.variables.units import joule, kelvin, meter
or import the appropriate units from sympy, e.g.
from sympy.physics.units import joule, kelvin, meter
In [10]:
from sympy.physics.units import joule, kelvin, meter, mole, pascal, second
Then you can define a custom variable with its name, description, domain, latex_name, unit, and an optional default value, e.g.:
In [11]:
class R_mol(Variable):
"""Molar gas constant."""
unit = joule/(kelvin*mole)
latex_name = 'R_{mol}'
default = 8.314472
The variables defined above hold information about their docstring, units, latex representations and default values if any. Each can be accessed by e.g.:
In [12]:
print(R_mol.__doc__)
print(R_mol.definition.unit)
print(R_mol.definition.latex_name)
print(R_mol.definition.default)
We will now define a few additional variables.
In [13]:
class P_g(Variable):
"""Pressure of gas."""
unit = pascal
class V_g(Variable):
"""Volume of gas."""
unit = meter**3
class n_g(Variable):
"""Amount of gas."""
unit = mole
class n_w(Variable):
"""Amount of water."""
unit = mole
class T_g(Variable):
"""Temperature of gas."""
unit = kelvin
class P_wa(Variable):
"""Partial pressure of water vapour in air."""
unit = pascal
latex_name = 'P_{wa}'
In [14]:
class Delta_Pwa(Variable):
"""Slope of saturated vapour pressure, $\partial P_{wa} / \partial T_g$"""
expr = Derivative(P_wa,T_g)
latex_name = r'\Delta'
In [15]:
Delta_Pwa.definition.unit
Out[15]:
In [16]:
Delta_Pwa.definition.expr
Out[16]:
In [17]:
generate_metadata_table([Delta_Pwa])
Out[17]:
In [18]:
class x(Variable):
"""Positive real variable."""
assumptions = {'positive': True, 'real': True}
print(solve(x**2 - 1))
In [19]:
from essm.equations import Equation
We will now define an equation representing the ideal gas law, based on the variables defined above:
In [20]:
class eq_ideal_gas_law(Equation):
"""Ideal gas law."""
expr = Eq(P_g*V_g, n_g*R_mol*T_g)
Note that whenever an equation is defined, its units are checked for consistency in the background and if they are not consistent, an error message will be printed. To illustrate this, we will try to define the above equation again, but omit temperature on the right hand side:
In [21]:
try:
class eq_ideal_gas_law(Equation):
"""Ideal gas law."""
expr = Eq(P_g*V_g, n_g*R_mol)
except Exception as exc1:
print(exc1)
The equation can be displayed in typesetted form, and the documentation string can be accessed in a similar way as for Variable:
In [22]:
display(eq_ideal_gas_law)
print(eq_ideal_gas_law.__doc__)
In [23]:
soln = solve(eq_ideal_gas_law, P_g, dict=True); print(soln)
If we want to define a new equation based on a manipulation of eq_ideal_gas_law we can specify that the parent of the new equation is eq_ideal_gas_law.definition:
In [24]:
class eq_Pg(eq_ideal_gas_law.definition):
"""Calculate pressure of ideal gas."""
expr = Eq(P_g, soln[0][P_g])
eq_Pg
Out[24]:
We can also have nested inheritance, if we now define another equation based on eq_Pg:
In [25]:
class eq_Pwa_nw(eq_Pg.definition):
"""Calculate vapour pressure from amount of water in gas."""
expr = Eq(P_wa, eq_Pg.rhs.subs(n_g, n_w))
eq_Pwa_nw
Out[25]:
In [26]:
eq_Pwa_nw.definition.__bases__
Out[26]:
In [27]:
[parent.name for parent in eq_Pwa_nw.definition.__bases__]
Out[27]:
In [28]:
[parent.expr for parent in eq_Pwa_nw.definition.__bases__]
Out[28]:
We can also use a pre-defined function to extract the set of parents:
In [29]:
from essm.variables.utils import get_parents
get_parents(eq_Pwa_nw)
Out[29]:
We can use another function to get all parents recursively:
In [30]:
from essm.variables.utils import get_allparents
get_allparents(eq_Pwa_nw)
Out[30]:
In [31]:
class eq_Pg1(eq_ideal_gas_law.definition):
"""Calculate pressure of ideal gas."""
from sympy import solve
soln = solve(eq_ideal_gas_law, P_g, dict=True); print(soln)
expr = Eq(P_g, soln[0][P_g])
eq_Pg1
Out[31]:
In [32]:
%time
eq_Pg.subs({R_mol: 8.314, T_g: 300, n_g: 0.1, V_g: 1})
Out[32]:
In [33]:
%time
eq_Pg1.subs({R_mol: 8.314, T_g: 300, n_g: 0.1, V_g: 1})
Out[33]:
There is actually no difference!
Empirical equations not only contain variables but also numbers. As an example, we will try to define the Clausius-Clapeyron equation for saturation vapour pressure in the following example, after defining a few additional variables used in this equation. $$ P_{wa} = 611 e^\frac{-M_w \lambda_E (1/T_g - 1/273)}{R_{mol}}$$
In [34]:
from sympy.physics.units import joule, kilogram
class lambda_E(Variable):
"""Latent heat of evaporation."""
unit = joule/kilogram
latex_name = '\\lambda_E'
default = 2.45e6
class M_w(Variable):
"""Molar mass of water."""
unit = kilogram/mole
default = 0.018
In [35]:
from sympy import exp
try:
class eq_Pwa_CC(Equation):
"""Clausius-Clapeyron P_wa as function of T_g.
\cite[Eq. B3]{hartmann_global_1994}
"""
expr = Eq(P_wa, 611.*exp(-M_w*lambda_E*(1/T_g - 1/273.)/R_mol))
except Exception as exc1:
print(exc1)
The unit mismatch reported in the error message stems from the fact that the numbers in the empirical equation actually need units. Since the term in the exponent has to be non-dimensional, the units of 611 must be the same as those of P_wa, i.e. pascal. The units of the subtraction term in the exponent must match, meaning that 273 needs units of kelvin. To avoid the error message, we can define the empirical numbers as internal variables to the equation we want to define:
In [36]:
class eq_Pwa_CC(Equation):
"""Clausius-Clapeyron P_wa as function of T_g.
Eq. B3 in :cite{hartmann_global_1994}
"""
class p_CC1(Variable):
"""Internal parameter of eq_Pwl."""
unit = pascal
latex_name = '611'
default = 611.
class p_CC2(Variable):
"""Internal parameter of eq_Pwl."""
unit = kelvin
latex_name = '273'
default = 273.
expr = Eq(P_wa, p_CC1*exp(-M_w*lambda_E*(1/T_g - 1/p_CC2)/R_mol))
In the above, we defined the latex representation of the empirical constants as their actual values, so the equation displays in the familiar way:
In [37]:
eq_Pwa_CC
Out[37]:
All default values of variables defined along with the variable definitions are stored in a dictionary that can be accessed as Variable.__defaults__. We can substitute the values from this dictionary into our empirical equation to plot saturation vapour pressure as a function of temperature:
In [38]:
expr = eq_Pwa_CC.subs(Variable.__defaults__)
print(expr)
xvar = T_g
p = plot_expr2((xvar, 273, 373), expr)
In [39]:
from sympy import Piecewise
expr = Eq(P_wa, Piecewise((0, T_a < 0), (eq_Pwa_CC.rhs, T_a >= 0)))
expr
Out[39]:
In [40]:
try:
class eq1(Equation):
"""Test"""
expr = Eq(P_wa, Piecewise((0, T_a < 0), (eq_Pwa_CC.rhs, T_a >= 0)))
display(eq1)
except Exception as e1:
print(e1)
If the above returns a dimension error, then unit checking for Piecewise has not been implemented yet.
Above, we defined Delta_Pwa as a variable that represents the partial derivative of P_wa with respect to T_g:
class Delta_Pwa(Variable):
"""Slope of saturated vapour pressure, $\partial P_{ws} / \partial T_g"""
expr = P_wa(T_g).diff(T_g)
#unit = pascal/kelvin
latex_name = r'\Delta'This definition can be accessed by typing Delta_Pwa.definition.expr. Example:
In [41]:
print(Delta_Pwa.definition.expr)
display(Eq(Delta_Pwa, Delta_Pwa.definition.expr))
We also defined the Clausius-Clapeyron approximation to $P_{wa}(T_g)$ as eq_Pwa_CC.
In [42]:
display(eq_Pwa_CC)
print(eq_Pwa_CC.__doc__)
If we want to substitute this approximation into Delta_Pwa.definition.expr, we need to use replace instead of subs and evaluate the derivative using doit():
In [43]:
expr = Eq(Delta_Pwa, Delta_Pwa.definition.expr.replace(P_wa, eq_Pwa_CC.rhs).doit())
display(expr)
p = plot_expr2((T_g, 273, 373), expr.subs(Variable.__defaults__))
If we only had the slope of the curve, we could take the integral to get the absolute value:
In [44]:
from sympy import Integral
class T_a1(Variable):
"""Air temperature"""
unit = kelvin
latex_name = r'T_{a1}'
class T_a2(Variable):
"""Air temperature"""
unit = kelvin
latex_name = r'T_{a2}'
class P_wa1(Variable):
"""P_wa at T1"""
unit = pascal
latex_name = r'P_{wa1}'
class eq_Pwa_Delta(Equation):
"""P_wa deduced from the integral of Delta"""
expr = Eq(P_wa, P_wa1 + Integral(Delta_Pwa, (T_g, T_a1, T_a2)))
display(eq_Pwa_Delta)
In [45]:
expr_Delta = eq_Pwa_CC.rhs.diff(T_g)
expr = Eq(P_wa, eq_Pwa_Delta.rhs.replace(Delta_Pwa, expr_Delta).doit())
vdict = Variable.__defaults__.copy()
vdict[T_a1] = 273.
vdict[P_wa1] = eq_Pwa_CC.rhs.subs(T_g, T_a1).subs(vdict)
display(expr.subs(vdict))
p = plot_expr2((T_a2, 273, 373), expr.subs(vdict))
In [46]:
from sympy.physics.units import convert_to, kilo, mega, joule, kilogram, meter, second, inch, hour
from essm.variables.units import SI_EXTENDED_DIMENSIONS, SI_EXTENDED_UNITS
value1 = 0.3
unit1 = inch/hour
print(value1*unit1)
unit2 = Variable.get_dimensional_expr(unit1).subs(SI_EXTENDED_DIMENSIONS)
print(convert_to(value1*unit1, unit2))
The below example creates a file called test_variable-definitions.py with all the variable definitions and
relevant imports used in this notebook.
Then, we re-import all the variables from the newly created file and create another file called
test_equation-definitions.py with all the equation definitions and relevant imports.
To re-use the equations in another notebook, just execute from test_equation-definitions.py import *
In [47]:
from essm._generator import EquationWriter, VariableWriter
In [48]:
StrPrinter._print_Quantity = lambda self, expr: str(expr.name) # displays long units (meter instead of m)
writer = VariableWriter(docstring='Variables defined in api_features.ipynb and dependencies.')
for variable in Variable.__registry__.keys():
writer.var(variable)
writer.write('test_variable_definitions.py')
StrPrinter._print_Quantity = lambda self, expr: str(expr.abbrev) # displays short units (m instead of meter)
In [49]:
from test_variable_definitions import *
Since we re-imported names that already had definitions associated with them, we got a warning for each of them before it was overwritten. This tells us that the same variable used before may now have a different meaning, which could introduce inconsistency in the notebook. Here, however, we know that they are all the same. Now that we have re-imported all variables, we can use the EquationWriter to generate
another python file with all the equations and variables they depend on, as shown below. The re-import was necessary so that the import statements do not point to __main__
for locally defined variables.
In [50]:
StrPrinter._print_Quantity = lambda self, expr: str(expr.name) # displays long units (meter instead of m)
writer = EquationWriter(docstring='Equations defined in api_features.ipynb and dependencies.')
eqs_with_deps = []
for eq in Equation.__registry__.keys():
parents = tuple(get_parents(eq))
if parents == ():
writer.eq(eq)
else:
eqs_with_deps.append(eq) # Equations with dependencies must be at the end
for eq in eqs_with_deps:
writer.eq(eq)
writer.write('test_equation_definitions.py')
StrPrinter._print_Quantity = lambda self, expr: str(expr.abbrev) # displays short units (m instead of meter)
These definitions are re-imported in examples_numerics.ipynb, where you can also find examples for numerical computations using essm equations and variables.