### Set up virtual environment and install pystencils

Source: https://github.com/pycodegen/pystencils/blob/master/CONTRIBUTING.md

Create a virtual environment and install the pystencils package in editable mode. This is recommended for local development.

```bash
mkvirtualenv pystencils
cd pystencils/
pip install -e .
```

--------------------------------

### Initialize pystencils fields

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/02_tutorial_basic_kernels.ipynb

Setup source and destination arrays and map them to pystencils fields.

```python
src_arr = np.zeros([20, 30])
dst_arr = np.zeros_like(src_arr)

dst, src = ps.fields(dst=dst_arr, src=src_arr)
```

--------------------------------

### Install pystencils

Source: https://github.com/pycodegen/pystencils/blob/master/README.md

Commands to install the pystencils package with various optional dependencies.

```bash
pip install pystencils[interactive]
```

```bash
pip install pystencils[interactive, gpu, doc]
```

--------------------------------

### GPU Data Handling Setup

Source: https://context7.com/pycodegen/pystencils/llms.txt

Initializes data handling for GPU execution, specifying the default target as GPU. Requires pystencils.

```python
# GPU data handling
dh_gpu = ps.create_data_handling(
    domain_size=(256, 256),
    default_target=ps.Target.GPU
)
src_gpu = dh_gpu.add_array("src")
dst_gpu = dh_gpu.add_array("dst")
```

--------------------------------

### Create Domain Kernel Example

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/kernel_compile_and_call.md

Demonstrates creating a domain kernel for a 2D stencil computation. This involves defining fields, an assignment, kernel configuration, generating the AST, compiling it, and executing the compiled kernel with NumPy arrays. Ensure necessary imports are present.

```python
>>> import pystencils as ps
>>> import numpy as np
>>> from pystencils.kernelcreation import create_domain_kernel
>>> from pystencils.node_collection import NodeCollection
>>> s, d = ps.fields('s, d: [2D]')
>>> assignment = ps.Assignment(d[0,0], s[0, 1] + s[0, -1] + s[1, 0] + s[-1, 0])
>>> kernel_config = ps.CreateKernelConfig(cpu_openmp=True)
>>> kernel_ast = create_domain_kernel(NodeCollection([assignment]), config=kernel_config)
>>> kernel = kernel_ast.compile()
>>> d_arr = np.zeros([5, 5])
>>> kernel(d=d_arr, s=np.ones([5, 5]))
>>> d_arr
array([[0., 0., 0., 0., 0.],
       [0., 4., 4., 4., 0.],
       [0., 4., 4., 4., 0.],
       [0., 4., 4., 4., 0.],
       [0., 0., 0., 0., 0.]])
```

--------------------------------

### Initialize PyStencils Printing

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_phasefield_dentritic_3D.ipynb

Initializes pretty printing for symbolic expressions in PyStencils. This setup is typically done once at the beginning of a session.

```python
from pystencils.session import *
sp.init_printing()
frac = sp.Rational
```

--------------------------------

### Pure Diffusion Example (Heat Equation)

Source: https://context7.com/pycodegen/pystencils/llms.txt

Applies discretization to a pure diffusion equation (heat equation) using ps.fd. Requires pystencils and sympy.

```python
# Pure diffusion example (heat equation)
heat_eq = ps.fd.transient(c) - ps.fd.diffusion(c, D)
heat_update = discretize(heat_eq)
```

--------------------------------

### Execute GPU Simulation

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_wave_equation.ipynb

Calls the GPU simulation function with a specified number of timesteps. Requires CuPy to be installed.

```python
if cupy:
    run_on_gpu(400)
```

--------------------------------

### Create Data Handling for CPU or GPU

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/03_tutorial_datahandling.ipynb

Initialize data handling with a `default_target` to specify whether arrays should be allocated on the CPU or GPU. This simplifies GPU simulation setup.

```python
if gpu is False:
    dh = ps.create_data_handling(domain_size=(30, 30), default_target=ps.Target.CPU)
else:
    dh = ps.create_data_handling(domain_size=(30, 30), default_target=ps.Target.GPU)
```

--------------------------------

### PyStencils Field Access Examples

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/field.md

Demonstrates creating a 2D vector field and accessing its components, including neighbors and central values. Shows how to shift accesses and change components using methods like get_shifted and at_index.

```python
vector_field_2d = fields("v(2): double[2D]")  # create a 2D vector field
>>> northern_neighbor_y_component = vector_field_2d[0, 1](1)
>>> northern_neighbor_y_component
v_N^1
>>> central_y_component = vector_field_2d(1)
>>> central_y_component
v_C^1
>>> central_y_component.get_shifted(1, 0)  # move the existing access
v_E^1
>>> central_y_component.at_index(0)  # change component
v_C^0
```

--------------------------------

### Create Animation with ffmpeg

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_wave_equation.ipynb

Generates an animation of the simulation results using matplotlib's animation module. Requires ffmpeg to be installed.

```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")
```

--------------------------------

### Implement Periodic Boundary Conditions

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/03_tutorial_datahandling.ipynb

Set up a data handling with periodicity enabled for specific dimensions. Use `dh.synchronization_function()` to get a function that copies data to implement periodicity.

```python
dh = ps.create_data_handling((4, 4), periodicity=(True, False))
dh.add_array("src")

# get copy function
copy_fct = dh.synchronization_function(['src'])

dh.fill('src', 0.0, ghost_layers=True)
dh.fill('src', 3.0, ghost_layers=False)
before_sync = dh.gather_array('src', ghost_layers=True).copy()

copy_fct() # copy over to get periodicity in x direction

after_sync = dh.gather_array('src', ghost_layers=True).copy()
plt.subplot(1,2,1)
plt.scalar_field(before_sync);
plt.title("Before")
plt.subplot(1,2,2)
plt.scalar_field(after_sync);
plt.title("After");
```

--------------------------------

### Generate Example Phase Field Data

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_plotting_and_animation.ipynb

Generates sample data for phase-field simulations. The last coordinate represents the fraction of a phase (0 to 1), and the sum over this dimension should be 1.

```python
def example_phase_field():
    if scipy_avail:
        from scipy.ndimage.filters import gaussian_filter

    shape=(80, 80)
    result = np.zeros(shape + (4,))
    result[20:40, 20:40, 0] = 1
    
    if scipy_avail:
        gaussian_filter(result[..., 0], sigma=3, output=result[..., 0])
    result[50:70, 30:50, 1] = 1
    if scipy_avail:
        gaussian_filter(result[..., 1], sigma=3, output=result[..., 1])

    result[:, :10, 2] = 1
    result[:, :, 3] = 1 - np.sum(result, axis=2)
    return result
data = example_phase_field()
```

--------------------------------

### NestedScopes Symbol Visibility Example

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/ast.md

Demonstrates how NestedScopes manages symbol visibility, tracking defined and free parameters across nested scopes. Use this to understand symbol access and definition rules.

```python
>>> s = NestedScopes()
>>> s.access_symbol("a")
>>> s.is_defined("a")
False
>>> s.free_parameters
{'a'}
>>> s.define_symbol("b")
>>> s.is_defined("b")
True
>>> s.push()
>>> s.is_defined_locally("b")
False
>>> s.define_symbol("c")
>>> s.pop()
>>> s.is_defined("c")
False
```

--------------------------------

### Checking for CuPy Installation

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/03_tutorial_datahandling.ipynb

Detects if the cupy library is available for GPU-accelerated computations.

```python
try:
    import cupy
    gpu = True
except ImportError:
    gpu = False
    cupy = None
    print('No cupy installed')
```

--------------------------------

### Initialize data handling and fields

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_phasefield_dentritic_3D.ipynb

Sets up the simulation domain and registers the necessary fields for the phase-field and temperature variables.

```python
target = ps.Target.GPU
gpu = target == ps.Target.GPU
domain_size = (25, 25, 25) if 'is_test_run' in globals() else (300, 300, 300)

dh = ps.create_data_handling(domain_size=domain_size, periodicity=True, default_target=target)
φ_field = dh.add_array('phi', latex_name='φ')
φ_field_tmp = dh.add_array('phi_tmp', latex_name='φ_tmp')
φ_delta_field = dh.add_array('phidelta', latex_name='φ_D')
t_field = dh.add_array('T')
t_field_tmp = dh.add_array('T_tmp')
```

--------------------------------

### Get Kernel Parameters

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_wave_equation.ipynb

Retrieves the list of parameters that the compiled kernel expects.

```python
ast.get_parameters()
```

--------------------------------

### Initialize data handling and fields

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/06_tutorial_phasefield_dentritic_growth.ipynb

Sets up the data handling environment and initializes the required arrays for the phase field and temperature.

```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')
```

--------------------------------

### Initialize SimplificationStrategy

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_assignment_collection.ipynb

Create a strategy instance and register simplification rules such as expansion and common subexpression elimination.

```python
strategy = ps.simp.SimplificationStrategy()
strategy.add(ps.simp.apply_to_all_assignments(sp.expand))
strategy.add(ps.simp.sympy_cse)
```

--------------------------------

### Generated C kernel code

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/03_tutorial_datahandling.ipynb

Example of low-level C code generated by pystencils for array operations.

```c
FUNC_PREFIX void kernel(double * RESTRICT  _data_dst, double * RESTRICT const _data_src, int64_t const _size_dst_0, int64_t const _size_dst_1, int64_t const _stride_dst_0, int64_t const _stride_dst_1, int64_t const _stride_src_0, int64_t const _stride_src_1)
{
   for (int64_t ctr_0 = 1; ctr_0 < _size_dst_0 - 1; ctr_0 += 1)
   {
      double * RESTRICT  _data_dst_00 = _data_dst + _stride_dst_0*ctr_0;
      double * RESTRICT _data_src_01 = _data_src + _stride_src_0*ctr_0 + _stride_src_0;
      double * RESTRICT _data_src_00 = _data_src + _stride_src_0*ctr_0;
      double * RESTRICT _data_src_0m1 = _data_src + _stride_src_0*ctr_0 - _stride_src_0;
      for (int64_t ctr_1 = 1; ctr_1 < _size_dst_1 - 1; ctr_1 += 1)
      {
         _data_dst_00[_stride_dst_1*ctr_1] = 0.25*_data_src_00[_stride_src_1*ctr_1 + _stride_src_1] + 0.25*_data_src_00[_stride_src_1*ctr_1 - _stride_src_1] + 0.25*_data_src_01[_stride_src_1*ctr_1] + 0.25*_data_src_0m1[_stride_src_1*ctr_1];
      }
   }
}
```

--------------------------------

### Configure and Run Custom GPU Kernels

Source: https://context7.com/pycodegen/pystencils/llms.txt

Shows how to check for GPU availability, configure custom block sizes, and execute kernels using cupy arrays.

```python
import pystencils as ps
import numpy as np

# Check GPU availability
try:
    import cupy
    gpu_available = True
except ImportError:
    gpu_available = False

if gpu_available:
    # Create fields with fixed size for optimal GPU code
    src_arr = np.random.rand(1024, 1024).astype(np.float64)
    dst_arr = np.zeros_like(src_arr)

    src, dst = ps.fields("src, dst: [2D]", src=src_arr, dst=dst_arr)
    update = ps.Assignment(dst[0, 0], (src[1,0] + src[-1,0] + src[0,1] + src[0,-1]) / 4)

    # Create GPU kernel with custom block size
    from pystencils.gpu import BlockIndexing

    config = ps.CreateKernelConfig(
        target=ps.Target.GPU,
        gpu_indexing=BlockIndexing,
        gpu_indexing_params={'blockSize': (32, 8, 1)}
    )

    gpu_ast = ps.create_kernel(update, config=config)
    gpu_kernel = gpu_ast.compile()

    # View generated CUDA code
    ps.show_code(gpu_ast)

    # Transfer arrays to GPU and execute
    src_gpu = cupy.asarray(src_arr)
    dst_gpu = cupy.asarray(dst_arr)

    gpu_kernel(src=src_gpu, dst=dst_gpu)

    # Transfer result back to CPU
    result = dst_gpu.get()
```

--------------------------------

### Process Information Properties

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/datahandling.md

Provides properties to check if the current process is the root process and to get the world rank.

```APIDOC
## is_root (property)

### Description
Returns `True` for exactly one process in the simulation, indicating it is the root process.

### Method
N/A (Property)

### Parameters
N/A

### Returns
bool - True if the current process is the root process, False otherwise.
```

```APIDOC
## world_rank (property)

### Description
Returns the rank (identifier) of the current process within the simulation's world communicator.

### Method
N/A (Property)

### Parameters
N/A

### Returns
int - The rank of the current process.
```

--------------------------------

### Get Inverse Direction Tuple

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/stencil.md

Returns the inverse (negative) of a given direction tuple. Use when you need the opposite direction.

```python
>>> inverse_direction((1, -1, 0))
(-1, 1, 0)
```

--------------------------------

### Periodic Boundary Conditions with Data Handling

Source: https://context7.com/pycodegen/pystencils/llms.txt

Configures data handling for periodic boundary conditions and demonstrates running a simulation loop with synchronization for periodicity. Requires pystencils.

```python
# Periodic boundary conditions
dh_periodic = ps.create_data_handling(
    domain_size=(100, 100),
    periodicity=(True, True)  # Periodic in both directions
)
src_p = dh_periodic.add_array("src")
sync = dh_periodic.synchronization_function(["src"])

# Run simulation with periodic boundaries
for step in range(100):
    sync()  # Copy ghost layers for periodicity
    dh_periodic.run_kernel(kernel)
    dh_periodic.swap("src", "dst")
```

--------------------------------

### Get Fields Accessed by Kernel

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_wave_equation.ipynb

Lists the fields that the compiled kernel accesses, along with their expected data types and dimensions.

```python
ast.fields_accessed
```

--------------------------------

### Define Evolution Equations

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/06_tutorial_phasefield_dentritic_growth.ipynb

Set up the symbolic equations for the phase field and temperature evolution using pystencils.

```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)))
]
```

--------------------------------

### Initialize PyStencils Printing and Rational Numbers

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/06_tutorial_phasefield_dentritic_growth.ipynb

Import necessary modules from pystencils and initialize pretty printing. Define a shorthand for rational numbers.

```python
from pystencils.session import *
import pystencils as ps
sp.init_printing()
frac = sp.Rational
```

--------------------------------

### Initialize Wave Equation Simulation Fields

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_wave_equation.ipynb

Sets up the initial condition for the wave equation simulation using a sine wave pattern. The previous and current time step fields are initialized with this pattern.

```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)
```

--------------------------------

### Transfer and Execute GPU Kernels

Source: https://context7.com/pycodegen/pystencils/llms.txt

Demonstrates the basic workflow for transferring data to a GPU, compiling a kernel, and executing it using the DataHandling class.

```python
dh_gpu.fill("src", 0.0)
dh_gpu.all_to_gpu()
config = ps.CreateKernelConfig(target=ps.Target.GPU)
gpu_kernel = ps.create_kernel(update, config=config).compile()
dh_gpu.run_kernel(gpu_kernel)
dh_gpu.all_to_cpu()

# Access results
result = dh_gpu.gather_array("dst", ghost_layers=False)
```

--------------------------------

### Get Offsets and Coefficients

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/stencil.md

Returns two lists: one with accessed offsets and another with their corresponding coefficients. The expression must not contain any nonlinear parts.

```python
>>> import pystencils as ps
>>> f = ps.fields("f(3) : double[2D]")
>>> coff = coefficients(2 * f[0, 1](1) + 3 * f[-1, 0](1))
```

--------------------------------

### Initialize SymPy

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/01_tutorial_getting_started.ipynb

Imports the SymPy library and enables LaTeX-formatted output.

```python
import sympy as sp
sp.init_printing()  # enable nice LaTeX output
```

--------------------------------

### Extract Outer Derivatives with diff_terms()

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/finite_differences.md

Use `diff_terms` to get a set of all outer derivatives present in an expression. This differs from `expr.atoms(Diff)` for nested derivatives.

```python
>>> x, y = sp.symbols("x, y")
>>> diff_terms( diff(x, 0, 0) )
{Diff(Diff(x, 0, -1), 0, -1)}
>>> diff_terms( diff(x, 0, 0) + y )
{Diff(Diff(x, 0, -1), 0, -1)}
```

--------------------------------

### Defining and Executing a Pystencils Kernel

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/03_tutorial_datahandling.ipynb

Demonstrates the variable-field-size workflow: defining fields, creating an update rule, compiling the kernel, and executing it on numpy arrays.

```python
# 1. field definitions
src_field, dst_field = ps.fields("src, dst:[2D]")

# 2. define update rule
update_rule = [ps.Assignment(lhs=dst_field[0, 0],
                             rhs=(src_field[1, 0] + src_field[-1, 0] + 
                                  src_field[0, 1] + src_field[0, -1]) / 4)]

# 3. compile update rule to function
kernel_function = ps.create_kernel(update_rule).compile()

# 4. create numpy arrays and call kernel
src_arr, dst_arr = np.random.rand(30, 30), np.zeros([30, 30])
```

--------------------------------

### Generate example scalar field

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_plotting_and_animation.ipynb

A helper function to generate a 2D scalar field for demonstration purposes. It uses NumPy for meshgrid and trigonometric functions.

```python
def example_scalar_field(t=0):
    x, y = np.meshgrid(np.linspace(0, 2 * np.pi, 100), np.linspace(0, 2 * np.pi, 100))
    z = np.cos(x + 0.1 * t) * np.sin(y + 0.1 * t) + 0.1 * x * y
    return z
```

--------------------------------

### CreateKernelConfig Parameters

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/kernel_compile_and_call.md

Configuration options for creating kernels, including settings for GPU indexing, default assignment simplifications, and optimization strategies.

```APIDOC
#### base_pointer_specification *: List[Iterable[str]] | List[Iterable[int]]* *= None*

Specification of how many and which intermediate pointers are created for a field access. For example [ (0), (2,3,)] creates on base pointer for coordinates 2 and 3 and writes the offset for coordinate zero directly in the field access. These specifications are defined dependent on the loop ordering. This function translates more readable version into the specification above. For more information see: [`pystencils.transformations.create_intermediate_base_pointer`](ast.md#pystencils.transformations.create_intermediate_base_pointer)

#### gpu_indexing *: str* *= 'block'

Either ‘block’ or ‘line’ , or custom indexing class, see [`pystencils.gpu.AbstractIndexing`](#pystencils.gpu.AbstractIndexing)

#### gpu_indexing_params *: MappingProxyType*

Dict with indexing parameters (constructor parameters of indexing class) e.g. for ‘block’ one can specify ‘{‘block_size’: (20, 20, 10) }’.

#### default_assignment_simplifications *: bool* *= False*

If `True` default simplifications are first performed on the Assignments. If problems occur during the simplification a warning will be thrown. Furthermore, it is essential to know that this is a two-stage process. The first stage of the process acts on the level of the [`pystencils.AssignmentCollection`](simplifications.md#pystencils.AssignmentCollection). In this part, `pystencil.simp.create_simplification_strategy` from pystencils.simplificationfactory will be used to apply optimisations like insertion of constants to remove pressure from the registers. Thus the first part of the optimisations can only be executed if an [`AssignmentCollection`](simplifications.md#pystencils.AssignmentCollection) is passed. The second part of the optimisation acts on the level of each Assignment individually. In this stage, all optimisations from `sympy.codegen.rewriting.optims_c99` are applied to each Assignment. Thus this stage can also be applied if a list of Assignments is passed.

#### cpu_prepend_optimizations *: List[Callable]*

List of extra optimizations to perform first on the AST.

#### use_auto_for_assignments *: bool* *= False*

If set to `True`, auto can be used in the generated code for data types. This makes the type system more robust.

#### index_fields *: List[[Field](field.md#pystencils.field.Field)]* *= None*

List of index fields, i.e. 1D fields with struct data type. If not `None`, `create_index_kernel` instead of [`create_domain_kernel`](#pystencils.kernelcreation.create_domain_kernel) is used.

#### coordinate_names *: Tuple[str, Any]* *= ('x', 'y', 'z')*

Name of the coordinate fields in the struct data type.

#### allow_double_writes *: bool* *= False*

If True, don’t check if every field is only written at a single location. This is required for example for kernels that are compiled with loop step sizes > 1, that handle multiple cells at once. Use with care!

#### skip_independence_check *: bool* *= False*

By default the assignment list is checked for read/write independence. This means fields are only written at locations where they are read. Doing so guarantees thread safety. In some cases e.g. for periodicity kernel, this can not be assured and does the check needs to be deactivated. Use with care!

#### *class* DataTypeFactory(dt)

Because of pickle, we need to have a nested class, instead of a lambda in 
__post_init__
```

--------------------------------

### Animate Vector Field Magnitude

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_plotting_and_animation.ipynb

Creates an animation visualizing the magnitude of a vector field. Requires `ffmpeg` to be installed. The `run_func` should return updated field data.

```python
if shutil.which("ffmpeg") is not None:
    animation = plt.vector_field_magnitude_animation(run_func, frames=60)
    ps.jupyter.display_as_html_video(animation)
else:
    print("No ffmpeg installed")
```

--------------------------------

### Execute Simulation

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/06_tutorial_phasefield_dentritic_growth.ipynb

Run the simulation loop, initialize fields, and visualize the results.

```python
timeloop(10)
init()
plot()
print(dh)
```

--------------------------------

### Animate surface plot

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_plotting_and_animation.ipynb

Generates an animated 3D surface plot of an evolving scalar field. This visualization requires `ffmpeg` to be installed and is displayed as an HTML video.

```python
if shutil.which("ffmpeg") is not None:
    animation = plt.surface_plot_animation(run_func, frames=60)
    ps.jupyter.display_as_html_video(animation)
else:
    print("No ffmpeg installed")
```

--------------------------------

### Initialize Fields for 2D Wave Equation

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_wave_equation.ipynb

Creates numpy arrays and converts them into PyStencils fields for simulation. Fields are given LaTeX names for better representation.

```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,)
```

--------------------------------

### Animate scalar field (video)

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_plotting_and_animation.ipynb

Creates an animated video of an evolving scalar field. Requires `ffmpeg` to be installed. The animation is displayed as an HTML video element.

```python
t = 0
def run_func():
    global t
    t += 1
    return example_scalar_field(t)
```

```python
if shutil.which("ffmpeg") is not None:
    plt.figure()
    animation = plt.scalar_field_animation(run_func, frames=60)
    ps.jupyter.display_as_html_video(animation)
else:
    print("No ffmpeg installed")
```

--------------------------------

### Create and Compile Kernel

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/04_tutorial_advection_diffusion.ipynb

Generates an Abstract Syntax Tree (AST) for the discretized equation and compiles it into an executable kernel.

```python
ast = ps.create_kernel([ps.Assignment(c_next.center(), discretization.subs(sp.Symbol("D"), 1))])
kernel = ast.compile()
```

--------------------------------

### Discretize the model

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/06_tutorial_phasefield_dentritic_growth.ipynb

Sets up the discretization parameters for the simulation.

```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,
```

--------------------------------

### Numba JIT stencil implementation

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_benchmark.ipynb

A just-in-time compiled version using Numba's `@jit` decorator with `nopython=True` for maximum performance. This implementation is only available if Numba is installed.

```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
```

--------------------------------

### Serial Data Handling with create_data_handling

Source: https://context7.com/pycodegen/pystencils/llms.txt

Sets up serial data handling for a 2D domain, adds arrays, initializes them, and runs a compiled kernel. Requires pystencils and numpy.

```python
import pystencils as ps
import numpy as np

# Create serial data handling for a 2D domain
dh = ps.create_data_handling(domain_size=(100, 100))

# Add arrays (automatically creates corresponding fields)
src = dh.add_array("src", values_per_cell=1)
dst = dh.add_array("dst", values_per_cell=1)

# Initialize arrays
dh.fill("src", 0.0)
dh.cpu_arrays["src"][50, 50] = 1.0  # Point source

# Create and run kernel
update = ps.Assignment(dst[0, 0], (src[1,0] + src[-1,0] + src[0,1] + src[0,-1]) / 4)
kernel = ps.create_kernel(update).compile()
dh.run_kernel(kernel)
```

--------------------------------

### Cython stencil implementation

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_benchmark.ipynb

A compiled version using Cython, disabling bounds checking and wraparound for potential performance gains. This requires Cython to be installed and the code to be compiled.

```cython
%%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
```

--------------------------------

### Define Benchmark Parameters

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_small_block_benchmark.ipynb

Sets up parameters for benchmarking, including the number of inner and outer repetitions and a list of domain sizes to test.

```python
inner_repeats = 100
outer_repeats = 5
sizes = [2**i for i in range(1, 8)]
sizes
```

--------------------------------

### Access Read-Only Array using Gather (Serial)

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/03_tutorial_datahandling.ipynb

Use `dh.gather_array()` to get a read-only copy of an array. Modifications to this copy do not affect the original data in serial simulations.

```python
read_only_copy = dh.gather_array('src', ps.make_slice[:, :], ghost_layers=False)
```

--------------------------------

### Configure Benchmark Functions

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_small_block_benchmark.ipynb

Creates a dictionary mapping benchmark names to partial functions, allowing easy selection and configuration of different benchmark scenarios.

```python
name_to_func = {
    'pure_extract': partial(benchmark_pure, extract_first=True),
    'pure_no_extract': partial(benchmark_pure, extract_first=False),
    'dh_serial': partial(benchmark_datahandling, parallel=False),
    'dh_parallel': partial(benchmark_datahandling, parallel=True),
}
```

--------------------------------

### Get Stencil Coefficients as Nested Lists

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/stencil.md

Returns stencil coefficients in the form of nested lists. The expression must not have any nonlinear part. Can optionally format output as a matrix.

```python
>>> import pystencils as ps
>>> f = ps.fields("f: double[2D]")
>>> coefficient_list(2 * f[0, 1] + 3 * f[-1, 0])
[[0, 0, 0], [3, 0, 0], [0, 2, 0]]
```

```python
>>> coefficient_list(2 * f[0, 1] + 3 * f[-1, 0], matrix_form=True)
Matrix([
[0, 2, 0],
[3, 0, 0],
[0, 0, 0]])
```

--------------------------------

### Initialize parallel block grid and kernel

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/03_tutorial_datahandling.ipynb

Sets up a parallel block grid with periodic boundaries and defines a simple stencil update rule.

```python
blocks = createUniformBlockGrid(blocks=(2,4,1), cellsPerBlock=(20, 10, 1), 
                                oneBlockPerProcess=False, periodic=(1, 1, 0))
dh = ParallelDataHandling(blocks, dim=2)
src_field = dh.add_array('src')
dst_field = dh.add_array('dst')

update_rule = [ps.Assignment(lhs=dst_field[0, 0 ],
                             rhs=(src_field[1, 0] + src_field[-1, 0] + 
                                  src_field[0, 1] + src_field[0, -1]) / 4)]
```

--------------------------------

### Import PyStencils Session and Utilities

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_small_block_benchmark.ipynb

Imports necessary components from pystencils.session, time, statistics, and functools for benchmarking and performance measurement.

```python
from pystencils.session import *
from time import perf_counter
from statistics import median
from functools import partial
```

--------------------------------

### Generate Example Vector Field Data

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_plotting_and_animation.ipynb

Creates a 2D vector field with shape and time-dependent components. The last dimension must have 2 entries for x and y components.

```python
def example_vector_field(t=0, shape=(40, 40)):
    result = np.empty(shape + (2,))
    x, y = np.meshgrid(np.linspace(0, 2 * np.pi, shape[0]), np.linspace(0, 2 * np.pi, shape[1]))
    result[..., 0] = np.cos(x + 0.1 * t) * np.sin(y + 0.1 * t) + 0.01 * x * y
    result[..., 1] = np.cos(0.001 * y)
    return result
```

--------------------------------

### Animate Vector Field Quivers

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_plotting_and_animation.ipynb

Creates an animation of the vector field using quivers. Requires `ffmpeg` to be installed for video generation. The `run_func` should return updated field data.

```python
t = 0
def run_func():
    global t
    t += 1
    return example_vector_field(t)
```

```python
if shutil.which("ffmpeg") is not None:
    plt.figure()
    animation = plt.vector_field_animation(run_func, frames=60)
    ps.jupyter.display_as_html_video(animation)
else:
    print("No ffmpeg installed")
```

--------------------------------

### Fill Array Manually

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_jupyter_extensions.ipynb

Manually fills the 'c' array with values using nested loops. Note that this example might be less efficient than using `dh.fill` for large arrays.

```python
for x in range(129):
    for y in range(258):
        dh.cpu_arrays['c'][x, y] = 1.0
```

--------------------------------

### Field Creation Methods

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/field.md

Methods for initializing different types of fields within pystencils.

```APIDOC
## Static Methods for Field Creation

### create_generic
Creates a generic field where the field size is not fixed.
- **field_name** (str) - Required - Symbolic name for the field
- **spatial_dimensions** (int) - Required - Number of spatial dimensions
- **dtype** (type) - Optional - Numpy data type
- **index_dimensions** (int) - Optional - Number of index dimensions
- **layout** (str/tuple) - Optional - Loop ordering
- **index_shape** (tuple) - Optional - Shape of index dimensions
- **field_type** (FieldType) - Optional - Type of field (GENERIC, INDEXED, BUFFER, STAGGERED)

### create_from_numpy_array
Creates a field based on the layout, data type, and shape of a given numpy array.
- **field_name** (str) - Required - Symbolic name
- **array** (ndarray) - Required - Numpy array to base field on
- **index_dimensions** (int) - Optional - Number of index dimensions
- **field_type** (FieldType) - Optional - Kind of field

### create_fixed_size
Creates a field with fixed sizes.
- **field_name** (str) - Required - Symbolic name
- **shape** (tuple) - Required - Overall shape of the array
- **index_dimensions** (int) - Optional - Number of index dimensions
- **dtype** (type) - Optional - Numpy data type
- **layout** (str) - Optional - Full layout of array
- **strides** (tuple) - Optional - Strides in bytes
- **field_type** (FieldType) - Optional - Kind of field
```

--------------------------------

### Create Data Handling and Add Array

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_jupyter_extensions.ipynb

Initializes data handling for a 2D domain with periodicity and adds a new array 'c', filling it with initial values.

```python
dh = ps.create_data_handling(domain_size=(256, 256), periodicity=True)
c_field = dh.add_array('c')
dh.fill("c", 0.0, ghost_layers=True)
```

--------------------------------

### Import Pytest and Skip if CuPy Not Available

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_phasefield_dentritic_3D.ipynb

This snippet is used to conditionally import and use CuPy within Pytest. It ensures that tests relying on CuPy will be skipped if the CuPy library is not installed.

```python
import pytest
pytest.importorskip('cupy')
```

--------------------------------

### Get Offset Name for Field Access in PyStencils

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/field.md

The offset_name property returns a string representation of the spatial offset for a field access, using directional abbreviations like 'NE' for North-East.

```python
f = fields("f: double[2D]")
>>> f[1, 1].offset_name  # north-east
'NE'
```

--------------------------------

### Initialize and Inspect Simulation

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_phasefield_dentritic_3D.ipynb

Basic initialization and inspection of the data handling object.

```python
init()
plot()
print(dh)
```

--------------------------------

### Initialize Data Handling and Fields

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/05_tutorial_phasefield_spinodal_decomposition.ipynb

Creates a DataHandling instance for managing numpy arrays and symbolic sympy fields for concentration (c) and chemical potential (μ) on a 2D periodic domain.

```python
dh = ps.create_data_handling(domain_size=(256, 256), periodicity=True)
μ_field = dh.add_array('mu', latex_name='μ')
c_field = dh.add_array('c')
```

--------------------------------

### Handle waLBerla Availability Error

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_small_block_benchmark.ipynb

This error message indicates that the 'walberla' module is not available, preventing the creation of parallel data handling objects. Ensure waLBerla is installed and accessible if parallel execution is required.

```text
ValueError: Cannot create parallel data handling because walberla module is not available
```

--------------------------------

### Import necessary libraries

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_plotting_and_animation.ipynb

Imports the PyStencils session and shutil for utility functions. Ensure scipy is available for certain functionalities.

```python
from pystencils.session import *

import shutil
```

```python
try:
    from scipy.ndimage.filters import gaussian_filter
    scipy_avail = True
except ImportError:
    scipy_avail = False
```

--------------------------------

### Implement Simulation Loop and Utilities

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/06_tutorial_phasefield_dentritic_growth.ipynb

Define the time loop, initialization, and plotting functions for the simulation.

```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()
```

--------------------------------

### Expand Derivatives Linearly with expand_diff_linear()

Source: https://github.com/pycodegen/pystencils/blob/master/doc/sphinx/finite_differences.md

Applies the linearity property of differentiation to expand `Diff` nodes. For example, `Diff(c*a + b)` becomes `c * Diff(a) + Diff(b)`. Specify `functions` and `constants` to control symbol behavior.

```python
# Example usage (conceptual, actual code not provided in source)
# expand_diff_linear(Diff(c*a + b, x)) might yield c*Diff(a, x) + Diff(b, x)
```

--------------------------------

### Compile Simulation Kernels

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/06_tutorial_phasefield_dentritic_growth.ipynb

Generate and compile the kernels for the phase field and temperature evolution.

```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()
```

--------------------------------

### Assemble phase-field evolution equation

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_phasefield_dentritic_3D.ipynb

Combines the functional derivatives into the final evolution equation for the phase field.

```python
dF_dφ = μ_b + sp.Piecewise((μ_if, nLen**2 > 1e-10), (0, True))
```

--------------------------------

### Create Data Handling for CPU Simulations

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/03_tutorial_datahandling.ipynb

Initialize a data handling object for CPU simulations, specifying the domain size. This object manages the creation and mapping of symbolic fields to NumPy arrays.

```python
dh = ps.create_data_handling(domain_size=(30, 30))
```

--------------------------------

### Import necessary libraries

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_benchmark.ipynb

Imports the required modules for PyStencils, timeit, and Cython. Ensure Cython is loaded as an extension.

```python
from pystencils.session import *
import timeit
%load_ext Cython
```

--------------------------------

### Perform Symbolic Operations

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/01_tutorial_getting_started.ipynb

Demonstrates basic arithmetic and manipulation of symbolic expressions.

```python
expr = x**2 * ( y + x + 5) + x**2
expr
```

```python
expr.expand()
```

```python
expr.factor()
```

```python
expr.subs(y, sp.cos(x))
```

--------------------------------

### Collect all implementations

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/demo_benchmark.ipynb

A dictionary mapping implementation names to their respective functions. This is used for easy iteration and selection during benchmarking.

```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
```

--------------------------------

### Clone the pystencils repository

Source: https://github.com/pycodegen/pystencils/blob/master/CONTRIBUTING.md

Clone your forked pystencils repository locally to begin development. Ensure you replace 'your-name' with your GitLab username.

```bash
git clone https://i10git.cs.fau.de/your-name/pystencils
```

--------------------------------

### Create Kernel with Default Ghost Layers

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/02_tutorial_basic_kernels.ipynb

Generate a kernel with automatic ghost layer detection based on neighbor accesses. PyStencils restricts the iteration region to ensure safe access.

```python
import pystencils as ps

kernel = ps.create_kernel(ps.Assignment(dst[0,0], src[2, 0] + src[-1, 0]))
ps.show_code(kernel)
```

--------------------------------

### Initialize Array using Data Handling (Serial)

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/03_tutorial_datahandling.ipynb

Use `dh.fill()` for initializing arrays, which is recommended over direct access for better compatibility with distributed systems.

```python
src_arr = dh.cpu_arrays['src']
src_arr.fill(0.0)
```

```python
dh.fill('src', 0.0)
```

--------------------------------

### Create and Compile Kernel

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_jupyter_extensions.ipynb

Defines an assignment operation and creates a kernel configuration. The kernel is then generated and compiled.

```python
ur = ps.Assignment(c_field[0, 0], c_field[1, 0])
config = ps.CreateKernelConfig(target=dh.default_target, cpu_openmp=False, skip_independence_check=True)
ast = ps.create_kernel(ur, config=config)
kernel = ast.compile()
```

--------------------------------

### Create and Compile Kernels with ps.create_kernel()

Source: https://context7.com/pycodegen/pystencils/llms.txt

Use `ps.create_kernel()` to transform symbolic assignments into an executable AST. Configuration options allow for CPU parallelization (OpenMP), SIMD vectorization, and GPU compilation. Generated C code can be viewed using `ps.show_code()`.

```python
import pystencils as ps
import numpy as np

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)

# Basic kernel creation and compilation
ast = ps.create_kernel(update)
kernel = ast.compile()

# Execute kernel with numpy arrays
src_arr = np.random.rand(100, 100)
dst_arr = np.zeros_like(src_arr)
kernel(src=src_arr, dst=dst_arr)

# Create kernel with OpenMP parallelization
config = ps.CreateKernelConfig(cpu_openmp=4)  # Use 4 threads
ast_parallel = ps.create_kernel(update, config=config)
kernel_parallel = ast_parallel.compile()

# Create kernel with SIMD vectorization
config_vec = ps.CreateKernelConfig(
    cpu_vectorize_info={'instruction_set': 'avx2', 'assume_aligned': True}
)
kernel_vectorized = ps.create_kernel(update, config=config_vec).compile()

# Create GPU kernel (requires cupy)
config_gpu = ps.CreateKernelConfig(target=ps.Target.GPU)
gpu_ast = ps.create_kernel(update, config=config_gpu)
gpu_kernel = gpu_ast.compile()

# View generated C code
ps.show_code(ast)
```

--------------------------------

### Conditional Kernel with @ps.kernel

Source: https://context7.com/pycodegen/pystencils/llms.txt

Demonstrates defining a kernel with a conditional expression using Python's ternary if-else syntax within the @ps.kernel decorator. Ensure necessary imports are present.

```python
# Conditional expressions with ternary if-else
@ps.kernel
def threshold_kernel(s):
    threshold = sp.Symbol("threshold")
    dst[0, 0] @= src[0, 0] if src[0, 0] > threshold else 0

kernel_thresh = ps.create_kernel(threshold_kernel).compile()
```

--------------------------------

### Execute Benchmarks and Collect Results

Source: https://github.com/pycodegen/pystencils/blob/master/tests/test_small_block_benchmark.ipynb

Iterates through defined domain sizes and benchmark functions, running each benchmark multiple times and storing the results in a dictionary.

```python
result = {'block_size': [],
          'name': [],
          'time': []}

for bs in sizes:
    print("Computing size ", bs)
    for name, func in name_to_func.items():
        for i in range(outer_repeats):
            time = func((bs, bs))
            result['block_size'].append(bs)
            result['name'].append(name)
            result['time'].append(time)
```

--------------------------------

### Define Simulation Parameters

Source: https://github.com/pycodegen/pystencils/blob/master/doc/notebooks/06_tutorial_phasefield_dentritic_growth.ipynb

Define the physical constants and parameters for the simulation.

```python
    Teq: 1.0,
    θzero: 0.2,
    sp.pi: sp.pi.evalf()
}
parameters
```