Skip to content
Snippets Groups Projects
Commit 61d1bae6 authored by Martin Bauer's avatar Martin Bauer
Browse files

Tests and documentation for derivative module

parent beec6d3e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 665 additions and 38 deletions
.gitlab-ci.yml
+ 1
− 1
View file @ 61d1bae6
......@@ -130,4 +130,4 @@ pages:
tags:
- docker
only:
- master@software/pystencils
- master@pycodegen/pystencils
doc/notebooks/04_tutorial_advection_diffusion.ipynb
+ 19
− 9
View file @ 61d1bae6
Source diff could not be displayed: it is too large. Options to address this: view the blob.
doc/notebooks/05_tutorial_phasefield_spinodal_decomposition.ipynb
+ 30
− 24
View file @ 61d1bae6
Source diff could not be displayed: it is too large. Options to address this: view the blob.
doc/notebooks/demo_derivatives.ipynb 0 → 100644
+ 579
− 0
View file @ 61d1bae6
%% Cell type:code id: tags:
``` python
from pystencils.session import *
```
%% Cell type:markdown id: tags:
# 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.
%% Cell type:code id: tags:
``` python
f = ps.fields("f: [2D]")
first_derivative_x = ps.fd.diff(f, 0)
first_derivative_x
```
%% Output
$\displaystyle {\partial_{0} {{f}_{(0,0)}}}$
D(f[0,0])
%% Cell type:markdown id: tags:
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:
%% Cell type:code id: tags:
``` python
discretize_2nd_order = ps.fd.Discretization2ndOrder(dx=sp.symbols("h"))
discretize_2nd_order(first_derivative_x)
```
%% Output
!
$\displaystyle \frac{- {{f}_{(-1,0)}} + {{f}_{(1,0)}}}{2 h}$
-f_W + f_E
──────────
2⋅h
%% Cell type:markdown id: tags:
Strictly speaking, derivative objects act on *field accesses*, not *fields*, that why there is a $(0,0)$ index at the field:
%% Cell type:code id: tags:
``` python
first_derivative_x
```
%% Output
$\displaystyle {\partial_{0} {{f}_{(0,0)}}}$
D(f[0,0])
%% Cell type:markdown id: tags:
Sometimes it might be useful to specify derivatives at an offset e.g.
%% Cell type:code id: tags:
``` python
derivative_offset = ps.fd.diff(f[0, 1], 0)
derivative_offset, discretize_2nd_order(derivative_offset)
```
%% Output
$\displaystyle \left( {\partial_{0} {{f}_{(0,1)}}}, \ \frac{- {{f}_{(-1,1)}} + {{f}_{(1,1)}}}{2 h}\right)$
⎛ -f_NW + f_NE⎞
⎜D(f[0,1]), ────────────⎟
⎝ 2⋅h ⎠
%% Cell type:markdown id: tags:
Another example with second order derivatives:
%% Cell type:code id: tags:
``` python
laplacian = ps.fd.diff(f, 0, 0) + ps.fd.diff(f, 1, 1)
laplacian
```
%% Output
$\displaystyle {\partial_{0} {\partial_{0} {{f}_{(0,0)}}}} + {\partial_{1} {\partial_{1} {{f}_{(0,0)}}}}$
D(D(f[0,0])) + D(D(f[0,0]))
%% Cell type:code id: tags:
``` python
discretize_2nd_order(laplacian)
```
%% Output
$\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}}$
f_W - 2⋅f_C + f_E f_S - 2⋅f_C + f_N
───────────────── + ─────────────────
2 2
h h
%% Cell type:markdown id: tags:
## Working with derivative terms
No automatic simplifications are done on derivative terms i.e. linearity relations or product rule are not applied automatically.
%% Cell type:code id: tags:
``` python
f, g = ps.fields("f, g :[2D]")
c = sp.symbols("c")
δ = ps.fd.diff
expr = δ( δ(f, 0) + δ(g, 0) + c + 5 , 0)
expr
```
%% Output
$\displaystyle {\partial_{0} (c + {\partial_{0} {{f}_{(0,0)}}} + {\partial_{0} {{g}_{(0,0)}}} + 5) }$
D(c + Diff(f_C, 0, -1) + Diff(g_C, 0, -1) + 5)
%% Cell type:markdown id: tags:
This nested term can not be discretized automatically.
%% Cell type:code id: tags:
``` python
try:
discretize_2nd_order(expr)
except ValueError as e:
print(e)
```
%% Output
Only derivatives with field or field accesses as arguments can be discretized
%% Cell type:markdown id: tags:
### Linearity
The following function expands all derivatives exploiting linearity:
%% Cell type:code id: tags:
``` python
ps.fd.expand_diff_linear(expr)
```
%% Output
$\displaystyle {\partial_{0} c} + {\partial_{0} {\partial_{0} {{f}_{(0,0)}}}} + {\partial_{0} {\partial_{0} {{g}_{(0,0)}}}}$
D(c) + D(D(f[0,0])) + D(D(g[0,0]))
%% Cell type:markdown id: tags:
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:
%% Cell type:code id: tags:
``` python
ps.fd.expand_diff_linear(expr, functions=(f[0,0], g[0, 0]))
```
%% Output
$\displaystyle {\partial_{0} {\partial_{0} {{f}_{(0,0)}}}} + {\partial_{0} {\partial_{0} {{g}_{(0,0)}}}}$
D(D(f[0,0])) + D(D(g[0,0]))
%% Cell type:code id: tags:
``` python
ps.fd.expand_diff_linear(expr, constants=[c])
```
%% Output
$\displaystyle {\partial_{0} {\partial_{0} {{f}_{(0,0)}}}} + {\partial_{0} {\partial_{0} {{g}_{(0,0)}}}}$
D(D(f[0,0])) + D(D(g[0,0]))
%% Cell type:markdown id: tags:
The expanded term can then be discretized:
%% Cell type:code id: tags:
``` python
discretize_2nd_order(ps.fd.expand_diff_linear(expr, constants=[c]))
```
%% Output
$\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}}$
f_W - 2⋅f_C + f_E g_W - 2⋅g_C + g_E
───────────────── + ─────────────────
2 2
h h
%% Cell type:markdown id: tags:
### Product rule
The next cells show how to apply product rule and its reverse:
%% Cell type:code id: tags:
``` python
expr = δ(f[0, 0] * g[0, 0], 0 )
expr
```
%% Output
$\displaystyle {\partial_{0} ({{f}_{(0,0)}} {{g}_{(0,0)}}) }$
D(f_C*g_C)
%% Cell type:code id: tags:
``` python
expanded_expr = ps.fd.expand_diff_products(expr)
expanded_expr
```
%% Output
$\displaystyle {{f}_{(0,0)}} {\partial_{0} {{g}_{(0,0)}}} + {{g}_{(0,0)}} {\partial_{0} {{f}_{(0,0)}}}$
f_C⋅D(g[0,0]) + g_C⋅D(f[0,0])
%% Cell type:code id: tags:
``` python
recombined_expr = ps.fd.combine_diff_products(expanded_expr)
recombined_expr
```
%% Output
$\displaystyle {\partial_{0} ({{f}_{(0,0)}} {{g}_{(0,0)}}) }$
D(f_C*g_C)
%% Cell type:code id: tags:
``` python
assert recombined_expr == expr
```
%% Cell type:markdown id: tags:
### 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.
%% Cell type:code id: tags:
``` python
k = sp.symbols("k")
expr = δ(k**3 + 2 * k, 0 )
expr
```
%% Output
$\displaystyle {\partial_{0} (k^{3} + 2 k) }$
D(k**3 + 2*k)
%% Cell type:code id: tags:
``` python
ps.fd.evaluate_diffs(expr, var=k)
```
%% Output
$\displaystyle 3 k^{2} + 2$
2
3⋅k + 2
doc/sphinx/tutorials.rst
+ 1
− 1
View file @ 61d1bae6
......@@ -15,6 +15,6 @@ It is a good idea to download them and run them directly to be able to play arou
/notebooks/05_tutorial_phasefield_spinodal_decomposition.ipynb
/notebooks/06_tutorial_phasefield_dentritic_growth.ipynb
/notebooks/demo_assignment_collection.ipynb
/notebooks/demo_derivatives.ipynb
/notebooks/demo_benchmark.ipynb
/notebooks/demo_wave_equation.ipynb
pystencils/fd/derivative.py
+ 5
− 1
View file @ 61d1bae6
import sympy as sp
from collections import namedtuple, defaultdict
from pystencils import Field
from pystencils.sympyextensions import normalize_product, prod
......@@ -214,6 +213,11 @@ def diff_terms(expr):
This function yields different results than 'expr.atoms(Diff)' when nested derivatives are in the expression,
since this function only returns the outer derivatives
Example:
>>> x, y = sp.symbols("x, y")
>>> diff_terms( diff(x, 0, 0) )
{Diff(Diff(x, 0, -1), 0, -1)}
"""
result = set()
......
pystencils/fd/spatial.py
+ 1
− 1
View file @ 61d1bae6
......@@ -11,6 +11,7 @@ from .derivation import FiniteDifferenceStencilDerivation
def fd_stencils_standard(indices, dx, fa):
order = len(indices)
assert all(i >= 0 for i in indices), "Can only discretize objects with (integer) subscripts"
if order == 1:
idx = indices[0]
return (fa.neighbor(idx, 1) - fa.neighbor(idx, -1)) / (2 * dx)
......@@ -122,7 +123,6 @@ def discretize_spatial(expr, dx, stencil=fd_stencils_standard):
def discretize_spatial_staggered(expr, dx, stencil=fd_stencils_standard):
def staggered_visitor(e, coordinate, sign):
if isinstance(e, Diff):
arg, *indices = diff_args(e)
......
pystencils_tests/test_derivative.py
+ 29
− 1
View file @ 61d1bae6
import sympy as sp
import pystencils as ps
from sympy.printing.latex import LatexPrinter
from pystencils.fd import *
from sympy.abc import a, b, t, x, y, z
def test_derivative_basic():
x, y, z, t = sp.symbols("x y z t")
d = diff
op1, op2, op3 = DiffOperator(), DiffOperator(target=x), DiffOperator(target=x, superscript=1)
......@@ -18,4 +19,31 @@ def test_derivative_basic():
assert diff_term == dx**2 + 2 * dx * dy + dy**2 + 1
assert DiffOperator.apply(diff_term, t) == d(t, x, x) + 2 * d(t, x, y) + d(t, y, y) + t
assert ps.fd.Diff(0) == 0
expr = ps.fd.diff(ps.fd.diff(x, 0, 0), 1, 1)
assert expr.get_arg_recursive() == x
assert expr.change_arg_recursive(y).get_arg_recursive() == y
def test_derivative_expand_collect():
original = Diff(x*y*z)
result = combine_diff_products(combine_diff_products(expand_diff_products(original))).expand()
assert original == result
original = -3 * y * z * Diff(x) + 2 * x * z * Diff(y)
result = expand_diff_products(combine_diff_products(original)).expand()
assert original == result
original = a + b * Diff(x ** 2 * y * z)
expanded = expand_diff_products(original)
collect_res = combine_diff_products(combine_diff_products(combine_diff_products(expanded)))
assert collect_res == original
def test_diff_expand_using_linearity():
eps = sp.symbols("epsilon")
funcs = [a, b]
test = Diff(eps * Diff(a+b))
result = expand_diff_linear(test, functions=funcs)
assert result == eps * Diff(Diff(a)) + eps * Diff(Diff(b))
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment