doc/notebooks/06_tutorial_phasefield_dentritic_growth.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from pystencils.session import *
+import pystencils as ps
+sp.init_printing()
+frac = sp.Rational
+```
+
+%% Cell type:markdown id: tags:
+
+# Tutorial 06: Phase-field simulation of dentritic solidification
+
+This is the second tutorial on phase field methods with pystencils. Make sure to read the previous tutorial first.
+
+In this tutorial we again implement a model described in **Programming Phase-Field Modelling** by S. Bulent Biner.
+This time we implement the model from chapter 4.7 that describes dentritic growth. So get ready for some beautiful snowflake pictures.
+
+We start again by adding all required arrays fields. This time we explicitly store the change of the phase variable φ in time, since the dynamics is calculated using an Allen-Cahn formulation where a term $\partial_t \phi$ occurs.
+
+%% Cell type:code id: tags:
+
+``` python
+dh = ps.create_data_handling(domain_size=(300, 300), periodicity=True,
+                             default_target=ps.Target.CPU)
+φ_field = dh.add_array('phi', latex_name='φ')
+φ_field_tmp = dh.add_array('phi_temp', latex_name='φ_temp')
+φ_delta_field = dh.add_array('phidelta', latex_name='φ_D')
+t_field = dh.add_array('T')
+t_field_tmp = dh.add_array('T_tmp')
+```
+
+%% Cell type:markdown id: tags:
+
+This model has a lot of parameters that are created here in a symbolic fashion.
+
+%% Cell type:code id: tags:
+
+``` python
+ε, m, δ, j, θzero, α, γ, Teq, κ, τ = sp.symbols("ε m δ j θ_0 α γ T_eq κ τ")
+εb = sp.Symbol("\\bar{\\epsilon}")
+
+φ = φ_field.center
+φ_tmp = φ_field_tmp.center
+T = t_field.center
+
+def f(φ, m):
+    return φ**4 / 4 - (frac(1, 2) - m/3) * φ**3 + (frac(1,4)-m/2)*φ**2
+
+free_energy_density = ε**2 / 2 * (ps.fd.Diff(φ,0)**2 + ps.fd.Diff(φ,1)**2 ) + f(φ, m)
+free_energy_density
+```
+
+%% Output
+
+    $\displaystyle \frac{{φ}_{(0,0)}^{4}}{4} - {φ}_{(0,0)}^{3} \left(\frac{1}{2} - \frac{m}{3}\right) + {φ}_{(0,0)}^{2} \left(\frac{1}{4} - \frac{m}{2}\right) + \frac{ε^{2} \left({\partial_{0} {φ}_{(0,0)}}^{2} + {\partial_{1} {φ}_{(0,0)}}^{2}\right)}{2}$
+       4                                  2 ⎛         2            2⎞
+    φ_C       3 ⎛1   m⎞      2 ⎛1   m⎞   ε ⋅⎝D(φ[0,0])  + D(φ[0,0]) ⎠
+    ──── - φ_C ⋅⎜─ - ─⎟ + φ_C ⋅⎜─ - ─⎟ + ────────────────────────────
+     4          ⎝2   3⎠        ⎝4   2⎠                2
+
+%% Cell type:markdown id: tags:
+
+The free energy is again composed of a bulk and interface part.
+
+%% Cell type:code id: tags:
+
+``` python
+plt.figure(figsize=(7,4))
+plt.sympy_function(f(φ, m=1), x_values=(-1.05, 1.5) )
+plt.xlabel("φ")
+plt.title("Bulk free energy");
+```
+
+%% Output
+
+
+
+%% Cell type:markdown id: tags:
+
+Compared to last tutorial, this bulk free energy has also two minima, but at different values.
+
+Another complexity of the model is its anisotropy. The gradient parameter $\epsilon$ depends on the interface normal.
+
+%% Cell type:code id: tags:
+
+``` python
+def σ(θ):
+    return 1 + δ * sp.cos(j * (θ - θzero))
+
+θ = sp.atan2(ps.fd.Diff(φ, 1), ps.fd.Diff(φ, 0))
+
+ε_val = εb * σ(θ)
+ε_val
+```
+
+%% Output
+
+    $\displaystyle \bar{\epsilon} \left(δ \cos{\left(j \left(- θ_{0} + \operatorname{atan_{2}}{\left({\partial_{1} {φ}_{(0,0)}},{\partial_{0} {φ}_{(0,0)}} \right)}\right) \right)} + 1\right)$
+    \bar{\epsilon}⋅(δ⋅cos(j⋅(-θ₀ + atan2(D(φ[0,0]), D(φ[0,0])))) + 1)
+
+%% Cell type:code id: tags:
+
+``` python
+def m_func(T):
+    return (α / sp.pi) * sp.atan(γ * (Teq - T))
+```
+
+%% Cell type:markdown id: tags:
+
+However, we just insert these parameters into the free energy formulation before doing the functional derivative, to make the dependence of $\epsilon$ on $\nabla \phi$ explicit.
+
+%% Cell type:code id: tags:
+
+``` python
+fe = free_energy_density.subs({
+    m: m_func(T),
+    ε: εb * σ(θ),
+})
+
+dF_dφ = ps.fd.functional_derivative(fe, φ)
+dF_dφ = ps.fd.expand_diff_full(dF_dφ, functions=[φ])
+dF_dφ
+```
+
+%% Output
+
+    $\displaystyle {φ}_{(0,0)}^{3} - \frac{{φ}_{(0,0)}^{2} α \operatorname{atan}{\left({T}_{(0,0)} γ - T_{eq} γ \right)}}{\pi} - \frac{3 {φ}_{(0,0)}^{2}}{2} + \frac{{φ}_{(0,0)} α \operatorname{atan}{\left({T}_{(0,0)} γ - T_{eq} γ \right)}}{\pi} + \frac{{φ}_{(0,0)}}{2} - \bar{\epsilon}^{2} δ^{2} \cos^{2}{\left(j θ_{0} - j \operatorname{atan_{2}}{\left({\partial_{1} {φ}_{(0,0)}},{\partial_{0} {φ}_{(0,0)}} \right)} \right)} {\partial_{0} {\partial_{0} {φ}_{(0,0)}}} - \bar{\epsilon}^{2} δ^{2} \cos^{2}{\left(j θ_{0} - j \operatorname{atan_{2}}{\left({\partial_{1} {φ}_{(0,0)}},{\partial_{0} {φ}_{(0,0)}} \right)} \right)} {\partial_{1} {\partial_{1} {φ}_{(0,0)}}} - 2 \bar{\epsilon}^{2} δ \cos{\left(j θ_{0} - j \operatorname{atan_{2}}{\left({\partial_{1} {φ}_{(0,0)}},{\partial_{0} {φ}_{(0,0)}} \right)} \right)} {\partial_{0} {\partial_{0} {φ}_{(0,0)}}} - 2 \bar{\epsilon}^{2} δ \cos{\left(j θ_{0} - j \operatorname{atan_{2}}{\left({\partial_{1} {φ}_{(0,0)}},{\partial_{0} {φ}_{(0,0)}} \right)} \right)} {\partial_{1} {\partial_{1} {φ}_{(0,0)}}} - \bar{\epsilon}^{2} {\partial_{0} {\partial_{0} {φ}_{(0,0)}}} - \bar{\epsilon}^{2} {\partial_{1} {\partial_{1} {φ}_{(0,0)}}}$
+              2                               2
+       3   φ_C ⋅α⋅atan(T_C⋅γ - T_eq⋅γ)   3⋅φ_C    φ_C⋅α⋅atan(T_C⋅γ - T_eq⋅γ)   φ_C                 2  2    2
+    φ_C  - ─────────────────────────── - ────── + ────────────────────────── + ─── - \bar{\epsilon} ⋅δ ⋅cos (j⋅θ₀ - j⋅atan2(
+                        π                  2                  π                 2
+    
+    
+                                                        2  2    2
+    D(φ[0,0]), D(φ[0,0])))⋅D(D(φ[0,0])) - \bar{\epsilon} ⋅δ ⋅cos (j⋅θ₀ - j⋅atan2(D(φ[0,0]), D(φ[0,0])))⋅D(D(φ[0,0])) - 2⋅\ba
+    
+    
+    
+               2                                                                            2
+    r{\epsilon} ⋅δ⋅cos(j⋅θ₀ - j⋅atan2(D(φ[0,0]), D(φ[0,0])))⋅D(D(φ[0,0])) - 2⋅\bar{\epsilon} ⋅δ⋅cos(j⋅θ₀ - j⋅atan2(D(φ[0,0])
+    
+    
+    
+                                               2                              2
+    , D(φ[0,0])))⋅D(D(φ[0,0])) - \bar{\epsilon} ⋅D(D(φ[0,0])) - \bar{\epsilon} ⋅D(D(φ[0,0]))
+    
+
+%% Cell type:markdown id: tags:
+
+Then we insert all the numeric parameters and discretize:
+
+%% Cell type:code id: tags:
+
+``` python
+discretize = ps.fd.Discretization2ndOrder(dx=0.03, dt=1e-5)
+parameters = {
+    τ: 0.0003,
+    κ: 1.8,
+    εb: 0.01,
+    δ: 0.02,
+    γ: 10,
+    j: 6,
+    α: 0.9,
+    Teq: 1.0,
+    θzero: 0.2,
+    sp.pi: sp.pi.evalf()
+}
+parameters
+```
+
+%% Output
+
+
+    $\displaystyle \left\{ \pi : 3.14159265358979, \  T_{eq} : 1.0, \  \bar{\epsilon} : 0.01, \  j : 6, \  α : 0.9, \  γ : 10, \  δ : 0.02, \  θ_{0} : 0.2, \  κ : 1.8, \  τ : 0.0003\right\}$
+    {π: 3.14159265358979, T_eq: 1.0, \bar{\epsilon}: 0.01, j: 6, α: 0.9, γ: 10, δ: 0.02, θ₀: 0.2, κ: 1.8, τ: 0.0003}
+
+%% Cell type:code id: tags:
+
+``` python
+dφ_dt = - dF_dφ / τ
+φ_eqs = ps.simp.sympy_cse_on_assignment_list([ps.Assignment(φ_delta_field.center,
+                                                            discretize(dφ_dt.subs(parameters)))])
+φ_eqs.append(ps.Assignment(φ_tmp, discretize(ps.fd.transient(φ) - φ_delta_field.center)))
+
+temperature_evolution = -ps.fd.transient(T) + ps.fd.diffusion(T, 1) + κ * φ_delta_field.center
+temperature_eqs = [
+    ps.Assignment(t_field_tmp.center, discretize(temperature_evolution.subs(parameters)))
+]
+```
+
+%% Cell type:markdown id: tags:
+
+When creating the kernels we pass as target (which may be 'Target.CPU' or 'Target.GPU') the default target of the target handling. This enables to switch to a GPU simulation just by changing the parameter of the data handling.
+
+The rest is similar to the previous tutorial.
+
+%% Cell type:code id: tags:
+
+``` python
+φ_kernel = ps.create_kernel(φ_eqs, cpu_openmp=True, target=dh.default_target).compile()
+temperature_kernel = ps.create_kernel(temperature_eqs, target=dh.default_target).compile()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+def timeloop(steps=200):
+    φ_sync = dh.synchronization_function(['phi'])
+    temperature_sync = dh.synchronization_function(['T'])
+    dh.all_to_gpu()  # this does nothing when running on CPU
+    for t in range(steps):
+        φ_sync()
+        dh.run_kernel(φ_kernel)
+        dh.swap(φ_field.name, φ_field_tmp.name)
+        temperature_sync()
+        dh.run_kernel(temperature_kernel)
+        dh.swap(t_field.name, t_field_tmp.name)
+    dh.all_to_cpu()
+    return dh.gather_array('phi')
+
+def init(nucleus_size=np.sqrt(5)):
+    for b in dh.iterate():
+        x, y = b.cell_index_arrays
+        x, y = x-b.shape[0]//2, y-b.shape[0]//2
+        bArr = (x**2 + y**2) < nucleus_size**2
+        b['phi'].fill(0)
+        b['phi'][(x**2 + y**2) < nucleus_size**2] = 1.0
+        b['T'].fill(0.0)
+
+def plot():
+    plt.subplot(1,3,1)
+    plt.scalar_field(dh.gather_array('phi'))
+    plt.title("φ")
+    plt.colorbar()
+    plt.subplot(1,3,2)
+    plt.title("T")
+    plt.scalar_field(dh.gather_array('T'))
+    plt.colorbar()
+    plt.subplot(1,3,3)
+    plt.title("∂φ")
+    plt.scalar_field(dh.gather_array('phidelta'))
+    plt.colorbar()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+ps.show_code(φ_kernel)
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+timeloop(10)
+init()
+plot()
+print(dh)
+```
+
+%% Output
+
+        Name|      Inner (min/max)|     WithGl (min/max)
+    ----------------------------------------------------
+           T|            (  0,  0)|            (  0,  0)
+       T_tmp|            (  0,  0)|            (  0,  0)
+         phi|            (  0,  1)|            (  0,  1)
+    phi_temp|            (  0,  0)|            (  0,  0)
+    phidelta|            (  0,  0)|            (  0,  0)
+    
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+result = None
+if 'is_test_run' not in globals():
+    ani = ps.plot.scalar_field_animation(timeloop, rescale=True, frames=600)
+    result = ps.jupyter.display_as_html_video(ani)
+result
+```
+
+%% Cell type:code id: tags:
+
+``` python
+assert np.isfinite(dh.max('phi'))
+assert np.isfinite(dh.max('T'))
+```

doc/notebooks/demo_assignment_collection.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from pystencils.session import *
+```
+
+%% Cell type:markdown id: tags:
+
+# Demo: Assignment collections and simplification
+
+
+## Assignment collections
+
+The assignment collection class helps to formulate and simplify assignments for numerical kernels.
+
+An ``AssignmentCollection`` is an ordered collection of assignments, together with an optional ordered collection of subexpressions, that are required to evaluate the main assignments. There are various simplification rules available that operate on ``AssignmentCollection``s.
+
+%% Cell type:markdown id: tags:
+
+We start by defining some stencil update rule. Here we also use the *pystencils* ``Field``, note however that the assignment collection module works purely on the *sympy* level.
+
+%% Cell type:code id: tags:
+
+``` python
+a,b,c = sp.symbols("a b c")
+f = ps.fields("f(2) : [2D]")
+g = ps.fields("g(2) : [2D]")
+
+a1 = ps.Assignment(g[0,0](1), (a**2 +b) * f[0,1] + \
+                  (a**2 - c) * f[1,0] + \
+                  (a**2 - 2*c) * f[-1,0] + \
+                  (a**2) * f[0, -1])
+
+a2 = ps.Assignment(g[0,0](0), (c**2 +b) * f[0,1] + \
+                  (c**2 - c) * f[1,0] + \
+                  (c**2 - 2*c) * f[-1,0] + \
+                  (c**2 - a**2) * f[0, -1])
+
+
+ac = ps.AssignmentCollection([a1, a2], subexpressions=[])
+ac
+```
+
+%% Output
+
+    AssignmentCollection: g_C^0, g_C^1 <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
+
+%% Cell type:markdown id: tags:
+
+*sympy* operations can be applied on an assignment collection: In this example we first expand the collection, then look for common subexpressions.
+
+%% Cell type:code id: tags:
+
+``` python
+expand_all = ps.simp.apply_to_all_assignments(sp.expand)
+expandedEc = expand_all(ac)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+ac_cse = ps.simp.sympy_cse(expandedEc)
+ac_cse
+```
+
+%% Output
+
+    AssignmentCollection: g_C^0, g_C^1 <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
+
+%% Cell type:markdown id: tags:
+
+Symbols occuring in assignment collections are classified into 3 categories:
+- ``free_symbols``: symbols that occur in right-hand-sides but never on left-hand-sides
+- ``bound_symbols``: symbols that occur on left-hand-sides
+- ``defined_symbols``: symbols that occur on left-hand-sides of a main assignment
+
+%% Cell type:code id: tags:
+
+``` python
+ac_cse.free_symbols
+```
+
+%% Output
+
+
+    $\displaystyle \left\{{f}_{(1,0)}^{0}, {f}_{(0,1)}^{0}, {f}_{(0,-1)}^{0}, {f}_{(-1,0)}^{0}, a, b, c\right\}$
+    {f_E__0, f_N__0, f_S__0, f_W__0, a, b, c}
+
+%% Cell type:code id: tags:
+
+``` python
+ac_cse.bound_symbols
+```
+
+%% Output
+
+
+    $\displaystyle \left\{{g}_{(0,0)}^{0}, {g}_{(0,0)}^{1}, \xi_{0}, \xi_{1}, \xi_{2}, \xi_{3}\right\}$
+    {g_C__0, g_C__1, ξ₀, ξ₁, ξ₂, ξ₃}
+
+%% Cell type:code id: tags:
+
+``` python
+ac_cse.defined_symbols
+```
+
+%% Output
+
+
+    $\displaystyle \left\{{g}_{(0,0)}^{0}, {g}_{(0,0)}^{1}\right\}$
+    {g_C__0, g_C__1}
+
+%% Cell type:markdown id: tags:
+
+Assignment collections can be splitted up, and merged together. For splitting, a list of symbols that occur on the left-hand-side in the main assignments has to be passed. The returned assignment collection only contains these main assignments together with all necessary subexpressions.
+
+%% Cell type:code id: tags:
+
+``` python
+ac_f0 = ac_cse.new_filtered([g(0)])
+ac_f1 = ac_cse.new_filtered([g(1)])
+ac_f1
+```
+
+%% Output
+
+    AssignmentCollection: g_C^1, <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
+
+%% Cell type:markdown id: tags:
+
+Note here that $\xi_4$ is no longer part of the subexpressions, since it is not used in the main assignment of $f_C^1$.
+
+If we merge both collections together, we end up with the original collection.
+
+%% Cell type:code id: tags:
+
+``` python
+ac_f0.new_merged(ac_f1)
+```
+
+%% Output
+
+    AssignmentCollection: g_C^0, g_C^1 <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
+
+%% Cell type:markdown id: tags:
+
+There is also a method that inserts all subexpressions into the main assignments. This is the inverse operation of common subexpression elimination.
+
+%% Cell type:code id: tags:
+
+``` python
+assert sp.simplify(ac_f0.new_without_subexpressions().main_assignments[0].rhs - a2.rhs) == 0
+ac_f0.new_without_subexpressions()
+```
+
+%% Output
+
+    AssignmentCollection: g_C^0, <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
+
+%% Cell type:markdown id: tags:
+
+To evaluate an assignment collection, use the ``lambdify`` method. It is very similar to *sympy*s ``lambdify`` function.
+
+%% Cell type:code id: tags:
+
+``` python
+evalFct = ac_cse.lambdify([f[0,1], f[1,0]],  # new parameters of returned function
+                          fixed_symbols={a:1, b:2, c:3, f[0,-1]: 4, f[-1,0]: 5}) # fix values of other symbols
+evalFct(2,1)
+```
+
+%% Output
+
+
+    $\displaystyle \left\{ {g}_{(0,0)}^{0} : 75, \  {g}_{(0,0)}^{1} : -17\right\}$
+    {g_C__0: 75, g_C__1: -17}
+
+%% Cell type:markdown id: tags:
+
+lambdify is rather slow for evaluation. The intended way to evaluate an assignment collection is *pystencils* i.e. create a fast kernel, that applies the update at every site of a structured grid. The collection can be directly passed to the `create_kernel` function.
+
+%% Cell type:code id: tags:
+
+``` python
+func = ps.create_kernel(ac_cse).compile()
+```
+
+%% Cell type:markdown id: tags:
+
+## Simplification Strategies
+
+In above examples, we already applied simplification rules to assignment collections. Simplification rules are functions that take, as a single argument, an assignment collection and return an modified/simplified copy of it. The ``SimplificationStrategy`` class holds a list of simplification rules and can apply all of them in the specified order. Additionally it provides useful printing and reporting functions.
+
+We start by creating a simplification strategy, consisting of the expand and CSE simplifications we have already applied above:
+
+%% Cell type:code id: tags:
+
+``` python
+strategy = ps.simp.SimplificationStrategy()
+strategy.add(ps.simp.apply_to_all_assignments(sp.expand))
+strategy.add(ps.simp.sympy_cse)
+```
+
+%% Cell type:markdown id: tags:
+
+This strategy can be applied to any assignment collection:
+
+%% Cell type:code id: tags:
+
+``` python
+strategy(ac)
+```
+
+%% Output
+
+    AssignmentCollection: g_C^0, g_C^1 <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
+
+%% Cell type:markdown id: tags:
+
+The strategy can also print the simplification results at each stage.
+The report contains information about the number of operations after each simplification as well as the runtime of each simplification routine.
+
+%% Cell type:code id: tags:
+
+``` python
+strategy.create_simplification_report(ac)
+```
+
+%% Output
+
+    <pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x147de3e90>
+
+%% Cell type:markdown id: tags:
+
+The strategy can also print the full collection after each simplification...
+
+%% Cell type:code id: tags:
+
+``` python
+strategy.show_intermediate_results(ac)
+```
+
+%% Output
+
+    <pystencils.simp.simplificationstrategy.SimplificationStrategy.show_intermediate_results.<locals>.IntermediateResults at 0x147e09c90>
+
+%% Cell type:markdown id: tags:
+
+... or only specific assignments for better readability
+
+%% Cell type:code id: tags:
+
+``` python
+strategy.show_intermediate_results(ac, symbols=[g(1)])
+```
+
+%% Output
+
+    <pystencils.simp.simplificationstrategy.SimplificationStrategy.show_intermediate_results.<locals>.IntermediateResults at 0x1265a1b90>
+
+%% Cell type:code id: tags:
+
+``` python
+```

doc/notebooks/demo_benchmark.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from pystencils.session import *
+import timeit
+%load_ext Cython
+```
+
+%% Cell type:markdown id: tags:
+
+# Demo: Benchmark numpy, Cython, pystencils
+
+In this notebook we compare and benchmark different ways of implementing a simple stencil kernel in Python.
+Our simple example computes the average of the four neighbors in 2D and stores it in a second array. To prevent out-of-bounds accesses, we skip the cells at the border and compute values only in the range `[1:-1, 1:-1]`
+
+%% Cell type:markdown id: tags:
+
+## Implementations
+
+The first implementation is a pure Python implementation:
+
+%% Cell type:code id: tags:
+
+``` python
+def avg_pure_python(src, dst):
+    for x in range(1, src.shape[0] - 1):
+        for y in range(1, src.shape[1] - 1):
+            dst[x, y] = (src[x + 1, y] + src[x - 1, y] +
+                         src[x, y + 1] + src[x, y - 1]) / 4
+```
+
+%% Cell type:markdown id: tags:
+
+Obviously, this will be a rather slow version, since the loops are written directly in Python.
+
+Next, we use *numpy* functions to delegate the looping to numpy. The first version uses the `roll` function to shift the array by one element in each direction. This version has to allocate a new array for each accessed neighbor.
+
+%% Cell type:code id: tags:
+
+``` python
+def avg_numpy_roll(src, dst):
+    neighbors = [np.roll(src, axis=a, shift=s) for a in (0, 1) for s in (-1, 1)]
+    np.divide(sum(neighbors), 4, out=dst)
+```
+
+%% Cell type:markdown id: tags:
+
+Using views, we can get rid of the additional copies:
+
+%% Cell type:code id: tags:
+
+``` python
+def avg_numpy_slice(src, dst):
+    dst[1:-1, 1:-1] = src[2:, 1:-1] + src[:-2, 1:-1] + \
+                      src[1:-1, 2:] + src[1:-1, :-2]
+```
+
+%% Cell type:markdown id: tags:
+
+To further optimize the kernel we switch to Cython, to get a compiled C version.
+
+%% Cell type:code id: tags:
+
+``` python
+%%cython
+import cython
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def avg_cython(object[double, ndim=2] src, object[double, ndim=2] dst):
+    cdef int xs, ys, x, y
+    xs, ys = src.shape
+    for x in range(1, xs - 1):
+        for y in range(1, ys - 1):
+            dst[x, y] = (src[x + 1, y] + src[x - 1, y] +
+                         src[x, y + 1] + src[x, y - 1]) / 4
+```
+
+%% Cell type:markdown id: tags:
+
+If available, we also try the numba just-in-time compiler
+
+%% Cell type:code id: tags:
+
+``` python
+try:
+    from numba import jit
+
+    @jit(nopython=True)
+    def avg_numba(src, dst):
+        dst[1:-1, 1:-1] = src[2:, 1:-1] + src[:-2, 1:-1] + \
+                          src[1:-1, 2:] + src[1:-1, :-2]
+except ImportError:
+    pass
+```
+
+%% Cell type:markdown id: tags:
+
+And finally we also create a *pystencils* version of the same stencil code:
+
+%% Cell type:code id: tags:
+
+``` python
+src, dst = ps.fields("src, dst: [2D]")
+
+update = ps.Assignment(dst[0,0],
+                       (src[1, 0] + src[-1, 0] + src[0, 1] + src[0, -1]) / 4)
+avg_pystencils = ps.create_kernel(update).compile()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+all_implementations = {
+    'pure Python': avg_pure_python,
+    'numpy roll': avg_numpy_roll,
+    'numpy slice': avg_numpy_slice,
+    'pystencils': avg_pystencils,
+}
+if 'avg_cython' in globals():
+    all_implementations['Cython'] = avg_cython
+if 'avg_numba' in globals():
+    all_implementations['numba'] = avg_numba
+```
+
+%% Cell type:markdown id: tags:
+
+## Benchmark functions
+
+We implement a short function to get in- and output arrays of a given shape and to measure the runtime.
+
+%% Cell type:code id: tags:
+
+``` python
+def get_arrays(shape):
+    in_arr = np.random.rand(*shape)
+    out_arr = np.empty_like(in_arr)
+    return in_arr, out_arr
+
+def do_benchmark(func, shape):
+    in_arr, out_arr = get_arrays(shape)
+    func(src=in_arr, dst=out_arr) # warmup
+    timer = timeit.Timer('f(src=src, dst=dst)', globals={'f': func, 'src': in_arr, 'dst': out_arr})
+    calls, time_taken = timer.autorange()
+    return time_taken / calls
+```
+
+%% Cell type:markdown id: tags:
+
+## Comparison
+
+%% Cell type:code id: tags:
+
+``` python
+plot_order = ['pystencils', 'Cython', 'numba', 'numpy slice', 'numpy roll', 'pure Python']
+plot_order = [p for p in plot_order if p in all_implementations]
+
+def bar_plot(*shape):
+    names = plot_order
+    runtimes = tuple(do_benchmark(all_implementations[name], shape) for name in names)
+    speedups = tuple(runtime / min(runtimes) for runtime in runtimes)
+    y_pos = np.arange(len(names))
+    labels = tuple(f"{name} ({round(speedup, 1)} x)" for name, speedup in zip(names, speedups))
+
+    plt.text(0.5, 0.5, f"Size {shape}", horizontalalignment='center', fontsize=16,
+             verticalalignment='center', transform=plt.gca().transAxes)
+    plt.barh(y_pos, runtimes, log=True)
+
+    plt.yticks(y_pos, labels);
+    plt.xlabel('Runtime of single iteration')
+
+plt.figure(figsize=(8, 8))
+
+plt.subplot(3, 1, 1)
+bar_plot(32, 32)
+
+plt.subplot(3, 1, 2)
+bar_plot(128, 128)
+
+plt.subplot(3, 1, 3)
+bar_plot(2048, 2048)
+```
+
+%% Output
+
+
+
+%% Cell type:markdown id: tags:
+
+All runtimes are plotted logarithmically. Numbers next to the labels show how much slower the version is than the fastest one.

doc/notebooks/demo_derivatives.ipynb

+%% 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_E - f_W
+    ─────────
+       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_NE - f_NW⎞
+    ⎜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{- 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}}$
+    -2⋅f_C + f_E + f_W   -2⋅f_C + f_N + f_S
+    ────────────────── + ──────────────────
+             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{- 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}}$
+    -2⋅f_C + f_E + f_W   -2⋅g_C + g_E + g_W
+    ────────────────── + ──────────────────
+             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/notebooks/demo_wave_equation.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import psutil
+
+from pystencils.session import *
+
+import shutil
+```
+
+%% Cell type:markdown id: tags:
+
+# Demo: Finite differences - 2D wave equation
+
+In this tutorial we show how to use the finite difference module of *pystencils* to solve a 2D wave equations. The time derivative is discretized by a simple forward Euler method.
+
+%% Cell type:markdown id: tags:
+
+$$ \frac{\partial^2 u}{\partial t^2} = \mbox{div} \left( q(x,y) \nabla u \right)$$
+
+%% Cell type:markdown id: tags:
+
+We begin by creating three *numpy* arrays for the current, the previous and the next timestep. Actually we will see later that two fields are enough, but let's keep it simple. From these *numpy* arrays we create *pystencils* fields to formulate our update rule.
+
+%% Cell type:code id: tags:
+
+``` python
+size = (60, 70) #  domain size
+u_arrays = [np.zeros(size), np.zeros(size), np.zeros(size)]
+
+u_fields = [ps.Field.create_from_numpy_array("u%s" % (name,), arr)
+            for name, arr in zip(["0", "1", "2"], u_arrays)]
+
+# Nicer display for fields
+for i, field in enumerate(u_fields):
+    field.latex_name = "u^{(%d)}" % (i,)
+```
+
+%% Cell type:markdown id: tags:
+
+*pystencils* contains already simple rules to discretize the a diffusion term. The time discretization is done manually.
+
+%% Cell type:code id: tags:
+
+``` python
+discretize = ps.fd.Discretization2ndOrder()
+
+def central2nd_time_derivative(fields):
+    f_next, f_current, f_last = fields
+    return (f_next[0, 0] - 2 * f_current[0, 0] + f_last[0, 0]) / discretize.dt**2
+
+rhs = ps.fd.diffusion(u_fields[1], 1)
+
+wave_eq = sp.Eq(central2nd_time_derivative(u_fields), discretize(rhs))
+
+wave_eq = sp.simplify(wave_eq)
+wave_eq
+```
+
+%% Output
+
+
+    $\displaystyle \frac{{u^{(0)}}_{(0,0)} - 2 {u^{(1)}}_{(0,0)} + {u^{(2)}}_{(0,0)}}{dt^{2}} = \frac{- 4 {u^{(1)}}_{(0,0)} + {u^{(1)}}_{(1,0)} + {u^{(1)}}_{(0,1)} + {u^{(1)}}_{(0,-1)} + {u^{(1)}}_{(-1,0)}}{dx^{2}}$
+    u_0_C - 2⋅u_1_C + u_2_C   -4⋅u_1_C + u_1_E + u_1_N + u_1_S + u_1_W
+    ─────────────────────── = ────────────────────────────────────────
+                2                                 2
+              dt                                dx
+
+%% Cell type:markdown id: tags:
+
+The explicit Euler scheme is now obtained by solving above equation with respect to $u_C^{next}$.
+
+%% Cell type:code id: tags:
+
+``` python
+u_next_C = u_fields[-1][0, 0]
+update_rule = ps.Assignment(u_next_C, sp.solve(wave_eq, u_next_C)[0])
+update_rule
+```
+
+%% Output
+
+
+    $\displaystyle {u^{(2)}}_{(0,0)} \leftarrow \frac{- 4 {u^{(1)}}_{(0,0)} dt^{2} + {u^{(1)}}_{(1,0)} dt^{2} + {u^{(1)}}_{(0,1)} dt^{2} + {u^{(1)}}_{(0,-1)} dt^{2} + {u^{(1)}}_{(-1,0)} dt^{2} + dx^{2} \left(- {u^{(0)}}_{(0,0)} + 2 {u^{(1)}}_{(0,0)}\right)}{dx^{2}}$
+                         2           2           2           2           2     2
+             - 4⋅u_1_C⋅dt  + u_1_E⋅dt  + u_1_N⋅dt  + u_1_S⋅dt  + u_1_W⋅dt  + dx ⋅(-u_0_C + 2⋅u_1_C)
+    u_2_C := ──────────────────────────────────────────────────────────────────────────────────────
+                                                        2
+                                                      dx
+
+%% Cell type:markdown id: tags:
+
+Before creating the kernel, we substitute numeric values for $dx$ and $dt$.
+Then a kernel is created just like in the last tutorial.
+
+%% Cell type:code id: tags:
+
+``` python
+update_rule = update_rule.subs({discretize.dx: 0.1, discretize.dt: 0.05})
+ast = ps.create_kernel(update_rule)
+kernel = ast.compile()
+
+ps.show_code(ast)
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+ast.get_parameters()
+```
+
+%% Output
+
+    [_data_u0, _data_u1, _data_u2]
+
+%% Cell type:code id: tags:
+
+``` python
+ast.fields_accessed
+```
+
+%% Output
+
+    {u0: double[60,70], u1: double[60,70], u2: double[60,70]}
+
+%% Cell type:markdown id: tags:
+
+To run simulation a suitable initial condition and boundary treatment is required. We chose an initial condition which is zero at the borders of the domain. The outermost layer is not changed by the update kernel, so we have an implicit homogenous Dirichlet boundary condition.
+
+%% Cell type:code id: tags:
+
+``` python
+X,Y = np.meshgrid( np.linspace(0, 1, size[1]), np.linspace(0,1, size[0]))
+Z = np.sin(2*X*np.pi) * np.sin(2*Y*np.pi)
+
+# Initialize the previous and current values with the initial function
+np.copyto(u_arrays[0], Z)
+np.copyto(u_arrays[1], Z)
+# The values for the next timesteps do not matter, since they are overwritten
+u_arrays[2][:, :] = 0
+```
+
+%% Cell type:markdown id: tags:
+
+One timestep now consists of applying the kernel once, then shifting the arrays.
+
+%% Cell type:code id: tags:
+
+``` python
+def run(timesteps=1):
+    for t in range(timesteps):
+        kernel(u0=u_arrays[0], u1=u_arrays[1], u2=u_arrays[2])
+        u_arrays[0], u_arrays[1], u_arrays[2] = u_arrays[1], u_arrays[2], u_arrays[0]
+    return u_arrays[2]
+```
+
+%% Cell type:markdown id: tags:
+
+Lets create an animation of the solution:
+
+%% Cell type:code id: tags:
+
+``` python
+if shutil.which("ffmpeg") is not None:
+    ani = plt.surface_plot_animation(run, zlim=(-1, 1))
+    ps.jupyter.display_as_html_video(ani)
+else:
+    print("No ffmpeg installed")
+```
+
+%% Output
+
+    No ffmpeg installed
+
+%% Cell type:code id: tags:
+
+``` python
+assert np.isfinite(np.max(u_arrays[2]))
+```
+
+%% Cell type:markdown id: tags:
+
+__Runing on GPU__
+
+We can also run the same kernel on the GPU, by using the ``cupy`` package.
+
+%% Cell type:code id: tags:
+
+``` python
+try:
+    import cupy
+except ImportError:
+    cupy=None
+    print('No cupy installed')
+
+
+res = None
+if cupy:
+    gpu_ast = ps.create_kernel(update_rule, target=ps.Target.GPU)
+    gpu_kernel = gpu_ast.compile()
+    res = ps.show_code(gpu_ast)
+res
+```
+
+%% Output
+
+
+
+%% Cell type:markdown id: tags:
+
+The run function has to be changed now slightly, since the data has to be transfered to the GPU first, then the kernel can be executed, and in the end the data has to be transfered back
+
+%% Cell type:code id: tags:
+
+``` python
+if cupy:
+    def run_on_gpu(timesteps=1):
+        # Transfer arrays to GPU
+        gpuArrs = [cupy.asarray(cpu_array) for cpu_array in u_arrays]
+
+        for t in range(timesteps):
+            gpu_kernel(u0=gpuArrs[0], u1=gpuArrs[1], u2=gpuArrs[2])
+            gpuArrs[0], gpuArrs[1], gpuArrs[2] = gpuArrs[1], gpuArrs[2], gpuArrs[0]
+
+        # Transfer arrays to CPU
+        for gpuArr, cpuArr in zip(gpuArrs, u_arrays):
+            cpuArr[:] = gpuArr.get()
+assert np.isfinite(np.max(u_arrays[2]))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+if cupy:
+    run_on_gpu(400)
+```

doc/sphinx/api.rst

+API Reference
+=============
+
+.. toctree::
+   :maxdepth: 3
+
+   kernel_compile_and_call.rst
+   enums.rst
+   simplifications.rst
+   datahandling.rst
+   configuration.rst
+   field.rst
+   stencil.rst
+   finite_differences.rst
+   plot.rst
+   ast.rst

doc/sphinx/ast.rst

+*********************************************
+For developers: AST Nodes and Transformations
+*********************************************
+
+
+AST Nodes
+=========
+
+.. automodule:: pystencils.astnodes
+   :members:
+
+
+
+Transformations
+===============
+
+
+.. automodule:: pystencils.transformations
+   :members:

doc/sphinx/configuration.rst

+*************
+Configuration
+*************
+
+.. automodule:: pystencils.cpu.cpujit
\ No newline at end of file

doc/sphinx/datahandling.rst

+************
+DataHandling
+************
+
+.. autoclass:: pystencils.datahandling.DataHandling
+   :members:
\ No newline at end of file

doc/sphinx/enums.rst

+************
+Enumerations
+************
+
+.. automodule:: pystencils.enums
+   :members:

doc/sphinx/field.rst

+*****
+Field
+*****
+
+.. automodule:: pystencils.field
+   :members:
\ No newline at end of file

doc/sphinx/finite_differences.rst

+******************
+Finite Differences
+******************
+
+.. automodule:: pystencils.fd
+   :members:
+

doc/sphinx/kernel_compile_and_call.rst

+*****************************************
+Creating and calling kernels from Python
+*****************************************
+
+
+Creating kernels
+----------------
+
+.. autofunction:: pystencils.create_kernel
+
+.. autoclass:: pystencils.CreateKernelConfig
+    :members:
+
+.. autofunction:: pystencils.kernelcreation.create_domain_kernel
+
+.. autofunction:: pystencils.kernelcreation.create_indexed_kernel
+
+.. autofunction:: pystencils.kernelcreation.create_staggered_kernel
+
+
+Code printing
+-------------
+
+.. autofunction:: pystencils.show_code
+
+
+GPU Indexing
+-------------
+
+.. autoclass:: pystencils.gpu.AbstractIndexing
+   :members:
+
+.. autoclass:: pystencils.gpu.BlockIndexing
+   :members:
+
+.. autoclass:: pystencils.gpu.LineIndexing
+   :members:

doc/sphinx/plot.rst

+**********************
+Plotting and Animation
+**********************
+
+.. automodule:: pystencils.plot
+   :members:
+
+

doc/sphinx/simplifications.rst

+***************************************
+Assignment Collection & Simplifications
+***************************************
+
+
+AssignmentCollection
+====================
+
+.. autoclass:: pystencils.AssignmentCollection
+   :members:
+
+
+SimplificationStrategy
+======================
+
+.. autoclass:: pystencils.simp.SimplificationStrategy
+    :members:
+
+Simplifications
+===============
+
+.. automodule:: pystencils.simp.simplifications
+    :members:
+
+Subexpression insertion
+=======================
+
+The subexpression insertions have the goal to insert subexpressions which will not reduce the number of FLOPs.
+For example a constant value kept as subexpression will lead to a new variable in the code which will occupy
+a register slot. On the other side a single variable could just be inserted in all assignments.
+
+.. automodule:: pystencils.simp.subexpression_insertion
+    :members:
+
+
+
+
+

doc/sphinx/stencil.rst

+*******
+Stencil
+*******
+
+.. automodule:: pystencils.stencil
+   :members:

doc/sphinx/tutorials.rst

+Tutorials
+=========
+
+These tutorials are a good place to start if you are new to pystencils.
+All tutorials and demos listed here are based on Jupyter notebooks that you can find in the pystencils repository.
+It is a good idea to download them and run them directly to be able to play around with the code.
+
+.. toctree::
+   :maxdepth: 1
+
+   /notebooks/01_tutorial_getting_started.ipynb
+   /notebooks/02_tutorial_basic_kernels.ipynb
+   /notebooks/03_tutorial_datahandling.ipynb
+   /notebooks/04_tutorial_advection_diffusion.ipynb
+   /notebooks/05_tutorial_phasefield_spinodal_decomposition.ipynb
+   /notebooks/06_tutorial_phasefield_dentritic_growth.ipynb
+   /notebooks/demo_assignment_collection.ipynb
+   /notebooks/demo_plotting_and_animation.ipynb
+   /notebooks/demo_derivatives.ipynb
+   /notebooks/demo_benchmark.ipynb
+   /notebooks/demo_wave_equation.ipynb

field.py

-from itertools import chain
-import numpy as np
-import sympy as sp
-from sympy.core.cache import cacheit
-from sympy.tensor import IndexedBase
-from pystencils.typedsymbol import TypedSymbol
-
-
-class Field:
-    """
-    With fields one can formulate stencil-like update rules on structured grids.
-    This Field class knows about the dimension, memory layout (strides) and optionally about the size of an array.
-
-    Creating Fields:
-
-        To create a field use one of the static create* members. There are two options:
-
-        1. create a kernel with fixed loop sizes i.e. the shape of the array is already known. This is usually the
-           case if just-in-time compilation directly from Python is done. (see :func:`Field.createFromNumpyArray`)
-        2. create a more general kernel that works for variable array sizes. This can be used to create kernels
-           beforehand for a library. (see :func:`Field.createGeneric`)
-
-    Dimensions:
-        A field has spatial and index dimensions, where the spatial dimensions come first.
-        The interpretation is that the field has multiple cells in (usually) two or three dimensional space which are
-        looped over. Additionally  N values are stored per cell. In this case spatialDimensions is two or three,
-        and indexDimensions equals N. If you want to store a matrix on each point in a two dimensional grid, there
-        are four dimensions, two spatial and two index dimensions: ``len(arr.shape) == spatialDims + indexDims``
-
-    Indexing:
-        When accessing (indexing) a field the result is a FieldAccess which is derived from sympy Symbol.
-        First specify the spatial offsets in [], then in case indexDimension>0 the indices in ()
-        e.g. ``f[-1,0,0](7)``
-
-    Example without index dimensions:
-        >>> a = np.zeros([10, 10])
-        >>> f = Field.createFromNumpyArray("f", a, indexDimensions=0)
-        >>> jacobi = ( f[-1,0] + f[1,0] + f[0,-1] + f[0,1] ) / 4
-
-    Example with index dimensions: LBM D2Q9 stream pull
-        >>> stencil = np.array([[0,0], [0,1], [0,-1]])
-        >>> src = Field.createGeneric("src", spatialDimensions=2, indexDimensions=1)
-        >>> dst = Field.createGeneric("dst", spatialDimensions=2, indexDimensions=1)
-        >>> for i, offset in enumerate(stencil):
-        ...     sp.Eq(dst[0,0](i), src[-offset](i))
-        Eq(dst_C^0, src_C^0)
-        Eq(dst_C^1, src_S^1)
-        Eq(dst_C^2, src_N^2)
-
-    """
-    @staticmethod
-    def createFromNumpyArray(fieldName, npArray, indexDimensions=0):
-        """
-        Creates a field based on the layout, data type, and shape of a given numpy array.
-        Kernels created for these kind of fields can only be called with arrays of the same layout, shape and type.
-        :param fieldName: symbolic name for the field
-        :param npArray: numpy array
-        :param indexDimensions: see documentation of Field
-        """
-        spatialDimensions = len(npArray.shape) - indexDimensions
-        if spatialDimensions < 1:
-            raise ValueError("Too many index dimensions. At least one spatial dimension required")
-
-        fullLayout = getLayoutFromNumpyArray(npArray)
-        spatialLayout = tuple([i for i in fullLayout if i < spatialDimensions])
-        assert len(spatialLayout) == spatialDimensions
-
-        strides = tuple([s // np.dtype(npArray.dtype).itemsize for s in npArray.strides])
-        shape = tuple([int(s) for s in npArray.shape])
-
-        return Field(fieldName, npArray.dtype, spatialLayout, shape, strides)
-
-    @staticmethod
-    def createGeneric(fieldName, spatialDimensions, dtype=np.float64, indexDimensions=0, layout=None):
-        """
-        Creates a generic field where the field size is not fixed i.e. can be called with arrays of different sizes
-        :param fieldName: symbolic name for the field
-        :param dtype: numpy data type of the array the kernel is called with later
-        :param spatialDimensions: see documentation of Field
-        :param indexDimensions: see documentation of Field
-        :param layout: tuple specifying the loop ordering of the spatial dimensions e.g. (2, 1, 0 ) means that
-        the outer loop loops over dimension 2, the second outer over dimension 1, and the inner loop
-        over dimension 0
-        """
-        if not layout:
-            layout = tuple(reversed(range(spatialDimensions)))
-        if len(layout) != spatialDimensions:
-            raise ValueError("Layout")
-        shapeSymbol = IndexedBase(TypedSymbol(Field.SHAPE_PREFIX + fieldName, Field.SHAPE_DTYPE), shape=(1,))
-        strideSymbol = IndexedBase(TypedSymbol(Field.STRIDE_PREFIX + fieldName, Field.STRIDE_DTYPE), shape=(1,))
-        totalDimensions = spatialDimensions + indexDimensions
-        shape = tuple([shapeSymbol[i] for i in range(totalDimensions)])
-        strides = tuple([strideSymbol[i] for i in range(totalDimensions)])
-        return Field(fieldName, dtype, layout, shape, strides)
-
-    def __init__(self, fieldName, dtype, layout, shape, strides):
-        """Do not use directly. Use static create* methods"""
-        self._fieldName = fieldName
-        self._dtype = numpyDataTypeToC(dtype)
-        self._layout = layout
-        self._shape = shape
-        self._strides = strides
-
-    @property
-    def spatialDimensions(self):
-        return len(self._layout)
-
-    @property
-    def indexDimensions(self):
-        return len(self._shape) - len(self._layout)
-
-    @property
-    def layout(self):
-        return self._layout
-
-    @property
-    def name(self):
-        return self._fieldName
-
-    @property
-    def shape(self):
-        return self._shape
-
-    @property
-    def spatialShape(self):
-        return self._shape[:self.spatialDimensions]
-
-    @property
-    def indexShape(self):
-        return self._shape[self.spatialDimensions:]
-
-    @property
-    def spatialStrides(self):
-        return self._strides[:self.spatialDimensions]
-
-    @property
-    def indexStrides(self):
-        return self._strides[self.spatialDimensions:]
-
-    @property
-    def strides(self):
-        return self._strides
-
-    @property
-    def dtype(self):
-        return self._dtype
-
-    def __repr__(self):
-        return self._fieldName
-
-    def neighbor(self, coordId, offset):
-        offsetList = [0] * self.spatialDimensions
-        offsetList[coordId] = offset
-        return Field.Access(self, tuple(offsetList))
-
-    def __getitem__(self, offset):
-        if type(offset) is np.ndarray:
-            offset = tuple(offset)
-        if type(offset) is str:
-            offset = tuple(directionStringToOffset(offset, self.spatialDimensions))
-        if type(offset) is not tuple:
-            offset = (offset,)
-        if len(offset) != self.spatialDimensions:
-            raise ValueError("Wrong number of spatial indices: "
-                             "Got %d, expected %d" % (len(offset), self.spatialDimensions))
-        return Field.Access(self, offset)
-
-    def __call__(self, *args, **kwargs):
-        center = tuple([0]*self.spatialDimensions)
-        return Field.Access(self, center)(*args, **kwargs)
-
-    def __hash__(self):
-        return hash((self._layout, self._shape, self._strides, self._dtype, self._fieldName))
-
-    def __eq__(self, other):
-        selfTuple = (self.shape, self.strides, self.name, self.dtype)
-        otherTuple = (other.shape, other.strides, other.name, other.dtype)
-        return selfTuple == otherTuple
-
-    PREFIX = "f"
-    STRIDE_PREFIX = PREFIX + "stride_"
-    SHAPE_PREFIX = PREFIX + "shape_"
-    STRIDE_DTYPE = "const int *"
-    SHAPE_DTYPE = "const int *"
-    DATA_PREFIX = PREFIX + "d_"
-
-    class Access(sp.Symbol):
-        def __new__(cls, name, *args, **kwargs):
-            obj = Field.Access.__xnew_cached_(cls, name, *args, **kwargs)
-            return obj
-
-        def __new_stage2__(self, field, offsets=(0, 0, 0), idx=None):
-            fieldName = field.name
-            offsetsAndIndex = chain(offsets, idx) if idx is not None else offsets
-            constantOffsets = not any([isinstance(o, sp.Basic) for o in offsetsAndIndex])
-
-            if not idx:
-                idx = tuple([0] * field.indexDimensions)
-
-            if constantOffsets:
-                offsetName = offsetToDirectionString(offsets)
-
-                if field.indexDimensions == 0:
-                    obj = super(Field.Access, self).__xnew__(self, fieldName + "_" + offsetName)
-                elif field.indexDimensions == 1:
-                    obj = super(Field.Access, self).__xnew__(self, fieldName + "_" + offsetName + "^" + str(idx[0]))
-                else:
-                    idxStr = ",".join([str(e) for e in idx])
-                    obj = super(Field.Access, self).__xnew__(self, fieldName + "_" + offsetName + "^" + idxStr)
-
-            else:
-                offsetName = "%0.10X" % (abs(hash(tuple(offsetsAndIndex))))
-                obj = super(Field.Access, self).__xnew__(self, fieldName + "_" + offsetName)
-
-            obj._field = field
-            obj._offsets = []
-            for o in offsets:
-                if isinstance(o, sp.Basic):
-                    obj._offsets.append(o)
-                else:
-                    obj._offsets.append(int(o))
-            obj._offsetName = offsetName
-            obj._index = idx
-
-            return obj
-
-        __xnew__ = staticmethod(__new_stage2__)
-        __xnew_cached_ = staticmethod(cacheit(__new_stage2__))
-
-        def __call__(self, *idx):
-            if self._index != tuple([0]*self.field.indexDimensions):
-                print(self._index, tuple([0]*self.field.indexDimensions))
-                raise ValueError("Indexing an already indexed Field.Access")
-
-            idx = tuple(idx)
-            if len(idx) != self.field.indexDimensions and idx != (0,):
-                raise ValueError("Wrong number of indices: "
-                                 "Got %d, expected %d" % (len(idx), self.field.indexDimensions))
-            return Field.Access(self.field, self._offsets, idx)
-
-        @property
-        def field(self):
-            return self._field
-
-        @property
-        def offsets(self):
-            return self._offsets
-
-        @property
-        def requiredGhostLayers(self):
-            return int(np.max(np.abs(self._offsets)))
-
-        @property
-        def nrOfCoordinates(self):
-            return len(self._offsets)
-
-        @property
-        def offsetName(self):
-            return self._offsetName
-
-        @property
-        def index(self):
-            return self._index
-
-        def _hashable_content(self):
-            superClassContents = list(super(Field.Access, self)._hashable_content())
-            t = tuple([*superClassContents, hash(self._field), self._index] + self._offsets)
-            return t
-
-
-def extractCommonSubexpressions(equations):
-    """
-    Uses sympy to find common subexpressions in equations and returns
-    them in a topologically sorted order, ready for evaluation.
-    Usually called before list of equations is passed to :func:`createKernel`
-    """
-    replacements, newEq = sp.cse(equations)
-    replacementEqs = [sp.Eq(*r) for r in replacements]
-    equations = replacementEqs + newEq
-    topologicallySortedPairs = sp.cse_main.reps_toposort([[e.lhs, e.rhs] for e in equations])
-    equations = [sp.Eq(*a) for a in topologicallySortedPairs]
-    return equations
-
-
-def getLayoutFromNumpyArray(arr):
-    """
-    Returns a list indicating the memory layout (linearization order) of the numpy array.
-    Example:
-    >>> getLayoutFromNumpyArray(np.zeros([3,3,3]))
-    [0, 1, 2]
-
-    In this example the loop over the zeroth coordinate should be the outermost loop,
-    followed by the first and second. Elements arr[x,y,0] and arr[x,y,1] are adjacent in memory.
-    Normally constructed numpy arrays have this order, however by stride tricks or other frameworks, arrays
-    with different memory layout can be created.
-    """
-    coordinates = list(range(len(arr.shape)))
-    return [x for (y, x) in sorted(zip(arr.strides, coordinates), key=lambda pair: pair[0], reverse=True)]
-
-
-def numpyDataTypeToC(dtype):
-    """Mapping numpy data types to C data types"""
-    if dtype == np.float64:
-        return "double"
-    elif dtype == np.float32:
-        return "float"
-    elif dtype == np.int32:
-        return "int"
-    raise NotImplementedError()
-
-
-def offsetComponentToDirectionString(coordinateId, value):
-    """
-    Translates numerical offset to string notation.
-    x offsets are labeled with east 'E' and 'W',
-    y offsets with north 'N' and 'S' and
-    z offsets with top 'T' and bottom 'B'
-    If the absolute value of the offset is bigger than 1, this number is prefixed.
-    :param coordinateId: integer 0, 1 or 2 standing for x,y and z
-    :param value: integer offset
-
-    Example:
-    >>> offsetComponentToDirectionString(0, 1)
-    'E'
-    >>> offsetComponentToDirectionString(1, 2)
-    '2N'
-    """
-    nameComponents = (('W', 'E'),  # west, east
-                      ('S', 'N'),  # south, north
-                      ('B', 'T'),  # bottom, top
-                      )
-    if value == 0:
-        result = ""
-    elif value < 0:
-        result = nameComponents[coordinateId][0]
-    else:
-        result = nameComponents[coordinateId][1]
-    if abs(value) > 1:
-        result = "%d%s" % (abs(value), result)
-    return result
-
-
-def offsetToDirectionString(offsetTuple):
-    """
-    Translates numerical offset to string notation.
-    For details see :func:`offsetComponentToDirectionString`
-    :param offsetTuple: 3-tuple with x,y,z offset
-
-    Example:
-    >>> offsetToDirectionString([1, -1, 0])
-    'SE'
-    >>> offsetToDirectionString(([-3, 0, -2]))
-    '2B3W'
-    """
-    names = ["", "", ""]
-    for i in range(len(offsetTuple)):
-        names[i] = offsetComponentToDirectionString(i, offsetTuple[i])
-    name = "".join(reversed(names))
-    if name == "":
-        name = "C"
-    return name
-
-
-def directionStringToOffset(directionStr, dim=3):
-    """
-    Reverse mapping of :func:`offsetToDirectionString`
-    :param directionStr: string representation of offset
-    :param dim: dimension of offset, i.e the length of the returned list
-
-    >>> directionStringToOffset('NW', dim=3)
-    array([-1,  1,  0])
-    >>> directionStringToOffset('NW', dim=2)
-    array([-1,  1])
-    >>> directionStringToOffset(offsetToDirectionString([3,-2,1]))
-    array([ 3, -2,  1])
-    """
-    offsetMap = {
-        'C': np.array([0, 0, 0]),
-
-        'W': np.array([-1, 0, 0]),
-        'E': np.array([1, 0, 0]),
-
-        'S': np.array([0, -1, 0]),
-        'N': np.array([0, 1, 0]),
-
-        'B': np.array([0, 0, -1]),
-        'T': np.array([0, 0, 1]),
-    }
-    offset = np.array([0, 0, 0])
-
-    while len(directionStr) > 0:
-        factor = 1
-        firstNonDigit = 0
-        while directionStr[firstNonDigit].isdigit():
-            firstNonDigit += 1
-        if firstNonDigit > 0:
-            factor = int(directionStr[:firstNonDigit])
-            directionStr = directionStr[firstNonDigit:]
-        curOffset = offsetMap[directionStr[0]]
-        offset += factor * curOffset
-        directionStr = directionStr[1:]
-    return offset[:dim]