from pystencils.session import *

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)

$\displaystyle \left( {\partial_{0} {f}_{(0,1)}}, \ \frac{{f}_{(1,1)} - {f}_{(-1,1)}}{2 h}\right)$

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$