ofuhrer / hpc4wc Goto Github PK

The apply_difussion() function in these two modules (day1/stencil2d.py, day5/stencil2d-gt4py-v0.py) do not work as expected when num_iter is even. This is because the swapping inside the function is only valid within the scope of the function. When num_iter is even, the swapping happens an odd number of times, so outside the function, the actual result is stored in in_field array and not in out_field. This means that when num_iter is even, out_field corresponds to num_iter - 1 iterations of the diffusion.

The following block code demonstrates this:

import numpy as np
from stencil2d import apply_diffusion

NX = NY = 128
NZ = 80
alpha = 1.0 / 32.0

in_field1 = np.zeros([NZ, NY, NX])
in_field1[:, NY//2 - NY//4 : NY//2 + NY//4, NX//2 - NX//4 : NX//2 + NX//4] = 1

in_field2 = np.copy(in_field1)

out_field1 = np.zeros_like(in_field1)
out_field2 = np.zeros_like(in_field1)

apply_diffusion(in_field1, out_field1, alpha=alpha, num_halo=2, num_iter=49)
apply_diffusion(in_field2, out_field2, alpha=alpha, num_halo=2, num_iter=50)

# out_field after 49 iterations equal to out_field after 50 iterations
np.allclose(out_field1, out_field2, rtol=1e-5, atol=1e-8, equal_nan=True)

The bug was found by @dymastery when working on the correctness checks in our project: #46.

A possible solution is to return out_field as we have done in our modified stencil2d and stencil2d_gt4py modules. This does not impact performance, as only a pointer is returned and no extra copies are made, and ensures that we keep track of the right field no matter the value of num_iter.

Computational Science and Engineering is in the (D-)MATH department, not MAVT.
Slides01-Admin.pdf, slide 9.