[1]:

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.

[2]:

f = ps.fields("f: [2D]")
first_derivative_x = ps.fd.diff(f, 0)
first_derivative_x

[2]:

$\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:

[3]:

discretize_2nd_order = ps.fd.Discretization2ndOrder(dx=sp.symbols("h"))
discretize_2nd_order(first_derivative_x)

!

[3]:

$\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:

[4]:

first_derivative_x

[4]:

$\displaystyle {\partial_{0} {{f}_{(0,0)}}}$

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

[5]:

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

[5]:

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

Another example with second order derivatives:

[6]:

laplacian = ps.fd.diff(f, 0, 0) + ps.fd.diff(f, 1, 1)
laplacian

[6]:

$\displaystyle {\partial_{0} {\partial_{0} {{f}_{(0,0)}}}} + {\partial_{1} {\partial_{1} {{f}_{(0,0)}}}}$
[7]:

discretize_2nd_order(laplacian)

[7]:

$\displaystyle \frac{{{f}_{(-1,0)}} - 2 {{f}_{(0,0)}} + {{f}_{(1,0)}}}{h^{2}} + \frac{{{f}_{(0,-1)}} - 2 {{f}_{(0,0)}} + {{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.

[8]:

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

expr = δ( δ(f, 0) +  δ(g, 0) + c + 5 , 0)
expr

[8]:

$\displaystyle {\partial_{0} (c + {\partial_{0} {{f}_{(0,0)}}} + {\partial_{0} {{g}_{(0,0)}}} + 5) }$

This nested term can not be discretized automatically.

[9]:

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:

[10]:

ps.fd.expand_diff_linear(expr)

[10]:

$\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:

[11]:

ps.fd.expand_diff_linear(expr, functions=(f[0,0], g[0, 0]))

[11]:

$\displaystyle {\partial_{0} {\partial_{0} {{f}_{(0,0)}}}} + {\partial_{0} {\partial_{0} {{g}_{(0,0)}}}}$
[12]:

ps.fd.expand_diff_linear(expr, constants=[c])

[12]:

$\displaystyle {\partial_{0} {\partial_{0} {{f}_{(0,0)}}}} + {\partial_{0} {\partial_{0} {{g}_{(0,0)}}}}$

The expanded term can then be discretized:

[13]:

discretize_2nd_order(ps.fd.expand_diff_linear(expr, constants=[c]))

[13]:

$\displaystyle \frac{{{f}_{(-1,0)}} - 2 {{f}_{(0,0)}} + {{f}_{(1,0)}}}{h^{2}} + \frac{{{g}_{(-1,0)}} - 2 {{g}_{(0,0)}} + {{g}_{(1,0)}}}{h^{2}}$

### Product rule¶

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

[14]:

expr = δ(f[0, 0] * g[0, 0], 0 )
expr

[14]:

$\displaystyle {\partial_{0} ({{f}_{(0,0)}} {{g}_{(0,0)}}) }$
[15]:

expanded_expr = ps.fd.expand_diff_products(expr)
expanded_expr

[15]:

$\displaystyle {{f}_{(0,0)}} {\partial_{0} {{g}_{(0,0)}}} + {{g}_{(0,0)}} {\partial_{0} {{f}_{(0,0)}}}$
[16]:

recombined_expr = ps.fd.combine_diff_products(expanded_expr)
recombined_expr

[16]:

$\displaystyle {\partial_{0} ({{f}_{(0,0)}} {{g}_{(0,0)}}) }$
[17]:

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.

[29]:

k = sp.symbols("k")
expr = δ(k**3 + 2 * k, 0 )
expr

[29]:

$\displaystyle {\partial_{0} (k^{3} + 2 k) }$
[31]:

ps.fd.evaluate_diffs(expr, var=k)

[31]:

$\displaystyle 3 k^{2} + 2$