Demo: Working with derivatives

import pystencils as ps
import sympy as sp

Demo: Working with derivatives#

Overview#

This notebook demonstrates how to formulate continuous differential operators in pystencils and automatically derive finite difference stencils from them.

Instead of using the built-in derivatives in sympy, pystencils comes with its own derivative objects. They represent spatial derivatives of pystencils fields.

f = ps.fields("f: [2D]")
first_derivative_x = ps.fd.diff(f, 0)
first_derivative_x
\[\displaystyle {\partial_{0} {f}_{(0,0)}}\]

This object is the derivative of the field \(f\) with respect to the first spatial coordinate \(x\). To get a finite difference approximation a discretization strategy is required:

discretize_2nd_order = ps.fd.Discretization2ndOrder(dx=sp.symbols("h"))
discretize_2nd_order(first_derivative_x)
\[\displaystyle \frac{{f}_{(1,0)} - {f}_{(-1,0)}}{2 h}\]

Strictly speaking, derivative objects act on field accesses, not fields, that why there is a \((0,0)\) index at the field:

first_derivative_x
\[\displaystyle {\partial_{0} {f}_{(0,0)}}\]

Sometimes it might be useful to specify derivatives at an offset e.g.

derivative_offset = ps.fd.diff(f[0, 1], 0)
derivative_offset, discretize_2nd_order(derivative_offset)

(Diff(f_N, 0, -1), (f_NE - f_NW)/(2*h))

Another example with second order derivatives:

laplacian = ps.fd.diff(f, 0, 0) + ps.fd.diff(f, 1, 1)
laplacian
\[\displaystyle {\partial_{0} {\partial_{0} {f}_{(0,0)}}} + {\partial_{1} {\partial_{1} {f}_{(0,0)}}}\]
discretize_2nd_order(laplacian)
\[\displaystyle \frac{- 2 {f}_{(0,0)} + {f}_{(1,0)} + {f}_{(-1,0)}}{h^{2}} + \frac{- 2 {f}_{(0,0)} + {f}_{(0,1)} + {f}_{(0,-1)}}{h^{2}}\]

Working with derivative terms#

No automatic simplifications are done on derivative terms i.e. linearity relations or product rule are not applied automatically.

f, g = ps.fields("f, g :[2D]")
c = sp.symbols("c")
δ = ps.fd.diff

expr = δ( δ(f, 0) +  δ(g, 0) + c + 5 , 0)
expr
\[\displaystyle {\partial_{0} (c + {\partial_{0} {f}_{(0,0)}} + {\partial_{0} {g}_{(0,0)}} + 5) }\]

This nested term can not be discretized automatically.

try:
    discretize_2nd_order(expr)
except ValueError as e:
    print(e)

Only derivatives with field or field accesses as arguments can be discretized

Linearity#

The following function expands all derivatives exploiting linearity:

ps.fd.expand_diff_linear(expr)
\[\displaystyle {\partial_{0} c} + {\partial_{0} {\partial_{0} {f}_{(0,0)}}} + {\partial_{0} {\partial_{0} {g}_{(0,0)}}}\]

The symbol \(c\) that was included is interpreted as a function by default. We can control the simplification behaviour by specifying all functions or all constants:

ps.fd.expand_diff_linear(expr, functions=(f[0,0], g[0, 0]))
\[\displaystyle {\partial_{0} {\partial_{0} {f}_{(0,0)}}} + {\partial_{0} {\partial_{0} {g}_{(0,0)}}}\]
ps.fd.expand_diff_linear(expr, constants=[c])
\[\displaystyle {\partial_{0} {\partial_{0} {f}_{(0,0)}}} + {\partial_{0} {\partial_{0} {g}_{(0,0)}}}\]

The expanded term can then be discretized:

discretize_2nd_order(ps.fd.expand_diff_linear(expr, constants=[c]))
\[\displaystyle \frac{- 2 {f}_{(0,0)} + {f}_{(1,0)} + {f}_{(-1,0)}}{h^{2}} + \frac{- 2 {g}_{(0,0)} + {g}_{(1,0)} + {g}_{(-1,0)}}{h^{2}}\]

Product rule#

The next cells show how to apply product rule and its reverse:

expr = δ(f[0, 0] * g[0, 0], 0 )
expr
\[\displaystyle {\partial_{0} ({f}_{(0,0)} {g}_{(0,0)}) }\]
expanded_expr = ps.fd.expand_diff_products(expr)
expanded_expr
\[\displaystyle {f}_{(0,0)} {\partial_{0} {g}_{(0,0)}} + {g}_{(0,0)} {\partial_{0} {f}_{(0,0)}}\]
recombined_expr = ps.fd.combine_diff_products(expanded_expr)
recombined_expr
\[\displaystyle {\partial_{0} ({f}_{(0,0)} {g}_{(0,0)}) }\]
assert recombined_expr == expr

Evaluate derivatives#

Arguments of derivative need not be to be fields, only when trying to discretize them. The next cells show how to transform them to sympy derivatives and evaluate them.

k = sp.symbols("k")
expr = δ(k**3 + 2 * k, 0 )
expr
\[\displaystyle {\partial_{0} (k^{3} + 2 k) }\]
ps.fd.evaluate_diffs(expr, var=k)
\[\displaystyle 3 k^{2} + 2\]
Contents