Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1440 additions and 111 deletions
pystencils_tests/test_datahandling.py → tests/test_datahandling.py
View file @ 5600b6b6
......@@ -6,8 +6,8 @@ import numpy as np
import pystencils as ps
from pystencils import create_data_handling, create_kernel
from pystencils.datahandling.pycuda import PyCudaArrayHandler
from pystencils.datahandling.pyopencl import PyOpenClArrayHandler
from pystencils.gpu.gpu_array_handler import GPUArrayHandler
from pystencils.enums import Target
try:
import pytest
......@@ -15,6 +15,12 @@ except ImportError:
import unittest.mock
pytest = unittest.mock.MagicMock()
try:
import cupy.cuda.runtime
device_numbers = range(cupy.cuda.runtime.getDeviceCount())
except ImportError:
device_numbers = []
SCRIPT_FOLDER = Path(__file__).parent.absolute()
INPUT_FOLDER = SCRIPT_FOLDER / "test_data"
......@@ -85,11 +91,7 @@ def access_and_gather(dh, domain_size):
def synchronization(dh, test_gpu=False):
field_name = 'comm_field_test'
if test_gpu:
try:
from pycuda import driver
import pycuda.autoinit
except ImportError:
return
pytest.importorskip("cupy")
field_name += 'Gpu'
dh.add_array(field_name, ghost_layers=1, dtype=np.int8, cpu=True, gpu=test_gpu)
......@@ -116,7 +118,7 @@ def synchronization(dh, test_gpu=False):
def kernel_execution_jacobi(dh, target):
test_gpu = target == 'gpu' or target == 'opencl'
test_gpu = target == Target.GPU
dh.add_array('f', gpu=test_gpu)
dh.add_array('tmp', gpu=test_gpu)
......@@ -132,7 +134,7 @@ def kernel_execution_jacobi(dh, target):
def jacobi():
dh.fields.tmp.center @= sum(dh.fields.f.neighbors(stencil)) / len(stencil)
kernel = create_kernel(jacobi, target).compile()
kernel = create_kernel(jacobi, config=ps.CreateKernelConfig(target=target)).compile()
for b in dh.iterate(ghost_layers=1):
b['f'].fill(42)
dh.run_kernel(kernel)
......@@ -211,26 +213,22 @@ def test_kernel():
for domain_shape in [(4, 5), (3, 4, 5)]:
dh = create_data_handling(domain_size=domain_shape, periodicity=True)
assert all(dh.periodicity)
kernel_execution_jacobi(dh, 'cpu')
kernel_execution_jacobi(dh, Target.CPU)
reduction(dh)
try:
import pycuda
import cupy
dh = create_data_handling(domain_size=domain_shape, periodicity=True)
kernel_execution_jacobi(dh, 'gpu')
kernel_execution_jacobi(dh, Target.GPU)
except ImportError:
pass
@pytest.mark.parametrize('target', ('cpu', 'gpu', 'opencl'))
@pytest.mark.parametrize('target', (Target.CPU, Target.GPU))
def test_kernel_param(target):
for domain_shape in [(4, 5), (3, 4, 5)]:
if target == 'gpu':
pytest.importorskip('pycuda')
if target == 'opencl':
pytest.importorskip('pyopencl')
from pystencils.opencl.opencljit import init_globally
init_globally()
if target == Target.GPU:
pytest.importorskip('cupy')
dh = create_data_handling(domain_size=domain_shape, periodicity=True, default_target=target)
kernel_execution_jacobi(dh, target)
......@@ -259,6 +257,20 @@ def test_add_arrays():
assert y == dh.fields['y']
@pytest.mark.parametrize('shape', [(17, 12), (7, 11, 18)])
@pytest.mark.parametrize('layout', ['zyxf', 'fzyx'])
def test_add_arrays_with_layout(shape, layout):
pytest.importorskip('cupy')
dh = create_data_handling(domain_size=shape, default_layout=layout, default_target=ps.Target.GPU)
f1 = dh.add_array("f1", values_per_cell=19)
dh.fill(f1.name, 1.0)
assert dh.cpu_arrays[f1.name].shape == dh.gpu_arrays[f1.name].shape
assert dh.cpu_arrays[f1.name].strides == dh.gpu_arrays[f1.name].strides
assert dh.cpu_arrays[f1.name].dtype == dh.gpu_arrays[f1.name].dtype
def test_get_kwarg():
domain_shape = (10, 10)
field_description = 'src, dst'
......@@ -269,7 +281,7 @@ def test_get_kwarg():
dh.fill("dst", 0.0, ghost_layers=True)
with pytest.raises(ValueError):
dh.add_array('src')
dh.add_array('src', values_per_cell=1)
ur = ps.Assignment(src.center, dst.center)
kernel = ps.create_kernel(ur).compile()
......@@ -280,22 +292,20 @@ def test_get_kwarg():
def test_add_custom_data():
pytest.importorskip('pycuda')
import pycuda.gpuarray as gpuarray
import pycuda.autoinit # noqa
pytest.importorskip('cupy')
import cupy as cp
def cpu_data_create_func():
return np.ones((2, 2), dtype=np.float64)
def gpu_data_create_func():
return gpuarray.zeros((2, 2), dtype=np.float64)
return cp.zeros((2, 2), dtype=np.float64)
def cpu_to_gpu_transfer_func(gpuarr, cpuarray):
gpuarr.set(cpuarray)
def gpu_to_cpu_transfer_func(gpuarr, cpuarray):
gpuarr.get(cpuarray)
cpuarray[:] = gpuarr.get()
dh = create_data_handling(domain_size=(10, 10))
dh.add_custom_data('custom_data',
......@@ -361,20 +371,11 @@ def test_load_data():
assert np.all(dh.cpu_arrays['dst2']) == 0
@pytest.mark.parametrize('target', ('gpu', 'opencl'))
def test_array_handler(target):
@pytest.mark.parametrize("device_number", device_numbers)
def test_array_handler(device_number):
size = (2, 2)
if target == 'gpu':
pytest.importorskip('pycuda')
array_handler = PyCudaArrayHandler()
if target == 'opencl':
pytest.importorskip('pyopencl')
import pyopencl as cl
from pystencils.opencl.opencljit import init_globally
init_globally()
ctx = cl.create_some_context(0)
queue = cl.CommandQueue(ctx)
array_handler = PyOpenClArrayHandler(queue)
pytest.importorskip('cupy')
array_handler = GPUArrayHandler(device_number)
zero_array = array_handler.zeros(size)
cpu_array = np.empty(size)
......@@ -388,8 +389,23 @@ def test_array_handler(target):
empty = array_handler.empty(size)
assert empty.strides == (16, 8)
empty = array_handler.empty(shape=size, layout=(1, 0))
empty = array_handler.empty(shape=size, order="F")
assert empty.strides == (8, 16)
random_array = array_handler.randn(size)
cpu_array = np.empty((20, 40), dtype=np.float64)
gpu_array = array_handler.to_gpu(cpu_array)
assert cpu_array.base is None
assert gpu_array.base is None
assert gpu_array.strides == cpu_array.strides
cpu_array2 = np.empty((20, 40), dtype=np.float64)
cpu_array2 = cpu_array2.swapaxes(0, 1)
gpu_array2 = array_handler.to_gpu(cpu_array2)
assert cpu_array2.base is not None
assert gpu_array2.base is not None
assert gpu_array2.strides == cpu_array2.strides
pystencils_tests/test_datahandling_parallel.py → tests/test_datahandling_parallel.py
View file @ 5600b6b6
import numpy as np
import waLBerla as wlb
import pystencils
from pystencils import make_slice
from tempfile import TemporaryDirectory
from pathlib import Path
from pystencils.boundaries import BoundaryHandling, Neumann
......@@ -10,7 +11,7 @@ from pystencils.slicing import slice_from_direction
from pystencils.datahandling.parallel_datahandling import ParallelDataHandling
from pystencils.datahandling import create_data_handling
from pystencils_tests.test_datahandling import (
from tests.test_datahandling import (
access_and_gather, kernel_execution_jacobi, reduction, synchronization, vtk_output)
SCRIPT_FOLDER = Path(__file__).parent.absolute()
......@@ -32,14 +33,12 @@ def test_access_and_gather():
dh = ParallelDataHandling(blocks, default_ghost_layers=2)
access_and_gather(dh, cells)
synchronization(dh, test_gpu=False)
if hasattr(wlb, 'cuda'):
if hasattr(wlb, 'gpu'):
synchronization(dh, test_gpu=True)
def test_gpu():
if not hasattr(wlb, 'cuda'):
print("Skip GPU tests because walberla was built without CUDA")
return
pytest.importorskip('waLBerla.gpu')
block_size = (4, 7, 1)
num_blocks = (3, 2, 1)
......@@ -57,24 +56,22 @@ def test_gpu():
np.testing.assert_equal(b['v'], 42)
def test_kernel():
for gpu in (True, False):
if gpu and not hasattr(wlb, 'cuda'):
print("Skipping CUDA tests because walberla was built without GPU support")
continue
@pytest.mark.parametrize('target', (pystencils.Target.CPU, pystencils.Target.GPU))
def test_kernel(target):
if target == pystencils.Target.GPU:
pytest.importorskip('waLBerla.gpu')
# 3D
blocks = wlb.createUniformBlockGrid(blocks=(3, 2, 4), cellsPerBlock=(3, 2, 5), oneBlockPerProcess=False)
dh = ParallelDataHandling(blocks)
kernel_execution_jacobi(dh, 'gpu')
reduction(dh)
# 3D
blocks = wlb.createUniformBlockGrid(blocks=(3, 2, 4), cellsPerBlock=(3, 2, 5), oneBlockPerProcess=False)
dh = ParallelDataHandling(blocks, default_target=target)
kernel_execution_jacobi(dh, target)
reduction(dh)
# 2D
blocks = wlb.createUniformBlockGrid(blocks=(3, 2, 1), cellsPerBlock=(3, 2, 1), oneBlockPerProcess=False)
dh = ParallelDataHandling(blocks, dim=2)
kernel_execution_jacobi(dh, 'gpu')
reduction(dh)
# 2D
blocks = wlb.createUniformBlockGrid(blocks=(3, 2, 1), cellsPerBlock=(3, 2, 1), oneBlockPerProcess=False)
dh = ParallelDataHandling(blocks, dim=2, default_target=target)
kernel_execution_jacobi(dh, target)
reduction(dh)
def test_vtk_output():
......@@ -88,7 +85,7 @@ def test_block_iteration():
num_blocks = (2, 2, 2)
blocks = wlb.createUniformBlockGrid(blocks=num_blocks, cellsPerBlock=block_size, oneBlockPerProcess=False)
dh = ParallelDataHandling(blocks, default_ghost_layers=2)
dh.add_array('v', values_per_cell=1, dtype=np.int64, ghost_layers=2, gpu=True)
dh.add_array('v', values_per_cell=1, dtype=np.int64, ghost_layers=2)
for b in dh.iterate():
b['v'].fill(1)
......@@ -111,10 +108,12 @@ def test_block_iteration():
def test_getter_setter():
pytest.importorskip('waLBerla.gpu')
block_size = (2, 2, 2)
num_blocks = (2, 2, 2)
blocks = wlb.createUniformBlockGrid(blocks=num_blocks, cellsPerBlock=block_size, oneBlockPerProcess=False)
dh = ParallelDataHandling(blocks, default_ghost_layers=2)
dh = ParallelDataHandling(blocks, default_ghost_layers=2, default_target=pystencils.Target.GPU)
dh.add_array('v', values_per_cell=1, dtype=np.int64, ghost_layers=2, gpu=True)
assert dh.shape == (4, 4, 4)
......@@ -132,22 +131,27 @@ def test_getter_setter():
def test_parallel_datahandling_boundary_conditions():
pytest.importorskip('waLBerla.cuda')
dh = create_data_handling(domain_size=(7, 7), periodicity=True, parallel=True, default_target="gpu")
src = dh.add_array('src')
src2 = dh.add_array('src2')
dh.fill("src", 0.0, ghost_layers=True)
dh.fill("src", 1.0, ghost_layers=False)
src_cpu = dh.add_array('src_cpu', gpu=False)
dh.fill("src_cpu", 0.0, ghost_layers=True)
dh.fill("src_cpu", 1.0, ghost_layers=False)
pytest.importorskip('waLBerla.gpu')
dh = create_data_handling(domain_size=(7, 7), periodicity=True, parallel=True,
default_target=pystencils.Target.GPU)
src = dh.add_array('src', values_per_cell=1)
dh.fill(src.name, 0.0, ghost_layers=True)
dh.fill(src.name, 1.0, ghost_layers=False)
src2 = dh.add_array('src2', values_per_cell=1)
src_cpu = dh.add_array('src_cpu', values_per_cell=1, gpu=False)
dh.fill(src_cpu.name, 0.0, ghost_layers=True)
dh.fill(src_cpu.name, 1.0, ghost_layers=False)
boundary_stencil = [(1, 0), (-1, 0), (0, 1), (0, -1)]
boundary_handling_cpu = BoundaryHandling(dh, src_cpu.name, boundary_stencil,
name="boundary_handling_cpu", target='cpu')
name="boundary_handling_cpu", target=pystencils.Target.CPU)
boundary_handling = BoundaryHandling(dh, src.name, boundary_stencil,
name="boundary_handling_gpu", target='gpu')
name="boundary_handling_gpu", target=pystencils.Target.GPU)
neumann = Neumann()
for d in ('N', 'S', 'W', 'E'):
......@@ -163,10 +167,11 @@ def test_parallel_datahandling_boundary_conditions():
boundary_handling()
dh.all_to_cpu()
for block in dh.iterate():
np.testing.assert_almost_equal(block["src_cpu"], block["src"])
np.testing.assert_almost_equal(block[src_cpu.name], block[src.name])
assert dh.custom_data_names == ('boundary_handling_cpuIndexArrays', 'boundary_handling_gpuIndexArrays')
dh.swap("src", "src2", gpu=True)
dh.swap(src.name, src2.name, gpu=True)
def test_save_data():
domain_shape = (2, 2)
......
tests/test_dot_printer.ipynb 0 → 100644
View file @ 5600b6b6
%% Cell type:code id: tags:
``` python
from pystencils.session import *
from pystencils.astnodes import Block, Conditional, SympyAssignment
```
%% Cell type:code id: tags:
``` python
src, dst = ps.fields("src, dst: double[2D]", layout='c')
true_block = Block([SympyAssignment(dst[0, 0], src[-1, 0])])
false_block = Block([SympyAssignment(dst[0, 0], src[1, 0])])
ur = [true_block, Conditional(dst.center() > 0.0, true_block, false_block)]
ast = ps.create_kernel(ur)
```
%% Cell type:code id: tags:
``` python
ps.show_code(ast)
```
%% Output
%% Cell type:code id: tags:
``` python
```
tests/test_dot_printer.py 0 → 100644
View file @ 5600b6b6
import pystencils as ps
from pystencils.astnodes import Block, Conditional, SympyAssignment
def test_dot_print():
src, dst = ps.fields("src, dst: double[2D]", layout='c')
true_block = Block([SympyAssignment(dst[0, 0], src[-1, 0])])
false_block = Block([SympyAssignment(dst[0, 0], src[1, 0])])
ur = [true_block, Conditional(dst.center() > 0.0, true_block, false_block)]
ast = ps.create_kernel(ur)
pystencils_tests/test_fast_approximation.py → tests/test_fast_approximation.py
View file @ 5600b6b6
......@@ -7,13 +7,13 @@ from pystencils.fast_approximation import (
def test_fast_sqrt():
pytest.importorskip('pycuda')
pytest.importorskip('cupy')
f, g = ps.fields("f, g: double[2D]")
expr = sp.sqrt(f[0, 0] + f[1, 0])
assert len(insert_fast_sqrts(expr).atoms(fast_sqrt)) == 1
assert len(insert_fast_sqrts([expr])[0].atoms(fast_sqrt)) == 1
ast_gpu = ps.create_kernel(ps.Assignment(g[0, 0], insert_fast_sqrts(expr)), target='gpu')
ast_gpu = ps.create_kernel(ps.Assignment(g[0, 0], insert_fast_sqrts(expr)), target=ps.Target.GPU)
ast_gpu.compile()
code_str = ps.get_code_str(ast_gpu)
assert '__fsqrt_rn' in code_str
......@@ -23,14 +23,14 @@ def test_fast_sqrt():
ac = ps.AssignmentCollection([expr], [])
assert len(insert_fast_sqrts(ac).main_assignments[0].atoms(fast_inv_sqrt)) == 1
ast_gpu = ps.create_kernel(insert_fast_sqrts(ac), target='gpu')
ast_gpu = ps.create_kernel(insert_fast_sqrts(ac), target=ps.Target.GPU)
ast_gpu.compile()
code_str = ps.get_code_str(ast_gpu)
assert '__frsqrt_rn' in code_str
def test_fast_divisions():
pytest.importorskip('pycuda')
pytest.importorskip('cupy')
f, g = ps.fields("f, g: double[2D]")
expr = f[0, 0] / f[1, 0]
assert len(insert_fast_divisions(expr).atoms(fast_division)) == 1
......@@ -38,7 +38,7 @@ def test_fast_divisions():
expr = 1 / f[0, 0] * 2 / f[0, 1]
assert len(insert_fast_divisions(expr).atoms(fast_division)) == 1
ast = ps.create_kernel(ps.Assignment(g[0, 0], insert_fast_divisions(expr)), target='gpu')
ast = ps.create_kernel(ps.Assignment(g[0, 0], insert_fast_divisions(expr)), target=ps.Target.GPU)
ast.compile()
code_str = ps.get_code_str(ast)
assert '__fdividef' in code_str
pystencils_tests/test_fd_derivation.ipynb → tests/test_fd_derivation.ipynb
View file @ 5600b6b6
%% Cell type:code id: tags:
``` python
from pystencils.session import *
from pystencils.fd.derivation import *
```
%% Cell type:markdown id: tags:
# 2D standard stencils
%% Cell type:code id: tags:
``` python
stencil = [(-1, 0), (1, 0), (0, -1), (0, 1), (0, 0)]
standard_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil)
f = ps.fields("f: [2D]")
standard_2d_00_res = standard_2d_00.get_stencil()
res = standard_2d_00_res.apply(f.center)
expected = f[-1, 0] - 2 * f[0, 0] + f[1, 0]
assert res == expected
```
%% Cell type:code id: tags:
``` python
assert standard_2d_00_res.accuracy == 2
assert not standard_2d_00_res.is_isotropic
standard_2d_00_res
```
%% Output
Finite difference stencil of accuracy 2, isotropic error: False
%% Cell type:code id: tags:
``` python
standard_2d_00.get_stencil().as_array()
```
%% Output
$\displaystyle \left[\begin{matrix}0 & 0 & 0\\1 & -2 & 1\\0 & 0 & 0\end{matrix}\right]$
⎡0 0 0⎤
⎢ ⎥
⎢1 -2 1⎥
⎢ ⎥
⎣0 0 0⎦
[[0, 0, 0], [1, -2, 1], [0, 0, 0]]
%% Cell type:markdown id: tags:
# 2D isotropic stencils
## second x-derivative
%% Cell type:code id: tags:
``` python
stencil = [(i, j) for i in (-1, 0, 1) for j in (-1, 0, 1)]
isotropic_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil)
isotropic_2d_00_res = isotropic_2d_00.get_stencil(isotropic=True)
assert isotropic_2d_00_res.is_isotropic
assert isotropic_2d_00_res.accuracy == 2
isotropic_2d_00_res
```
%% Output
Finite difference stencil of accuracy 2, isotropic error: True
%% Cell type:code id: tags:
``` python
isotropic_2d_00_res.as_array()
```
%% Output
$\displaystyle \left[\begin{matrix}\frac{1}{12} & - \frac{1}{6} & \frac{1}{12}\\\frac{5}{6} & - \frac{5}{3} & \frac{5}{6}\\\frac{1}{12} & - \frac{1}{6} & \frac{1}{12}\end{matrix}\right]$
⎡1/12 -1/6 1/12⎤
⎢ ⎥
⎢5/6 -5/3 5/6 ⎥
⎢ ⎥
⎣1/12 -1/6 1/12⎦
[[1/12, -1/6, 1/12], [5/6, -5/3, 5/6], [1/12, -1/6, 1/12]]
%% Cell type:code id: tags:
``` python
plt.figure(figsize=(2,2))
isotropic_2d_00_res.visualize()
```
%% Output
%% Cell type:code id: tags:
``` python
expected_result = sp.Array([[1, -2, 1], [10, -20, 10], [1, -2, 1]]) / 12
assert expected_result == isotropic_2d_00_res.as_array()
```
%% Cell type:markdown id: tags:
## Isotropic laplacian
%% Cell type:code id: tags:
``` python
isotropic_2d_11 = FiniteDifferenceStencilDerivation((1,1), stencil)
isotropic_2d_11_res = isotropic_2d_11.get_stencil(isotropic=True)
iso_laplacian = isotropic_2d_00_res.as_array() + isotropic_2d_11_res.as_array()
iso_laplacian
```
%% Output
$\displaystyle \left[\begin{matrix}\frac{1}{6} & \frac{2}{3} & \frac{1}{6}\\\frac{2}{3} & - \frac{10}{3} & \frac{2}{3}\\\frac{1}{6} & \frac{2}{3} & \frac{1}{6}\end{matrix}\right]$
⎡1/6 2/3 1/6⎤
⎢ ⎥
⎢2/3 -10/3 2/3⎥
⎢ ⎥
⎣1/6 2/3 1/6⎦
[[1/6, 2/3, 1/6], [2/3, -10/3, 2/3], [1/6, 2/3, 1/6]]
%% Cell type:code id: tags:
``` python
expected_result = sp.Array([[1, 4, 1], [4, -20, 4], [1, 4, 1]]) / 6
assert iso_laplacian == expected_result
```
%% Cell type:markdown id: tags:
# stencils for staggered fields
%% Cell type:code id: tags:
``` python
half = sp.Rational(1, 2)
fd_points_ex = (
(half, 0),
(-half, 0),
(half, 1),
(half, -1),
(-half, 1),
(-half, -1)
)
assert set(fd_points_ex) == set(FiniteDifferenceStaggeredStencilDerivation("E", 2).stencil)
fd_points_ey = (
(0, half),
(0, -half),
(-1,-half),
(-1, half),
(1, -half),
(1, half)
)
assert set(fd_points_ey) == set(FiniteDifferenceStaggeredStencilDerivation("N",2).stencil)
fd_points_c = (
(half, half),
(-half, -half),
(half, -half),
(-half, half)
)
assert set(fd_points_c) == set(FiniteDifferenceStaggeredStencilDerivation("NE",2).stencil)
assert len(FiniteDifferenceStaggeredStencilDerivation("E",3).points) == 10
assert len(FiniteDifferenceStaggeredStencilDerivation("NE",3).points) == 12
assert len(FiniteDifferenceStaggeredStencilDerivation("TNE",3).points) == 8
```
%% Cell type:code id: tags:
``` python
c = ps.fields("c: [2D]")
c3 = ps.fields("c3: [3D]")
assert FiniteDifferenceStaggeredStencilDerivation("E", 2, (0,)).apply(c.center) == c[1, 0] - c[0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("W", 2, (0,)).apply(c.center) == c[0, 0] - c[-1, 0]
assert FiniteDifferenceStaggeredStencilDerivation("N", 2, (1,)).apply(c.center) == c[0, 1] - c[0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("S", 2, (1,)).apply(c.center) == c[0, 0] - c[0, -1]
assert FiniteDifferenceStaggeredStencilDerivation("E", 3, (0,)).apply(c3.center) == c3[1, 0, 0] - c3[0, 0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("W", 3, (0,)).apply(c3.center) == c3[0, 0, 0] - c3[-1, 0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("N", 3, (1,)).apply(c3.center) == c3[0, 1, 0] - c3[0, 0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("S", 3, (1,)).apply(c3.center) == c3[0, 0, 0] - c3[0, -1, 0]
assert FiniteDifferenceStaggeredStencilDerivation("T", 3, (2,)).apply(c3.center) == c3[0, 0, 1] - c3[0, 0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("B", 3, (2,)).apply(c3.center) == c3[0, 0, 0] - c3[0, 0, -1]
assert FiniteDifferenceStaggeredStencilDerivation("S", 2, (0,)).apply(c.center) == \
(c[1, 0] + c[1, -1] - c[-1, 0] - c[-1, -1])/4
assert FiniteDifferenceStaggeredStencilDerivation("NE", 2, (0,)).apply(c.center) + \
FiniteDifferenceStaggeredStencilDerivation("NE", 2, (1,)).apply(c.center) == c[1, 1] - c[0, 0]
assert FiniteDifferenceStaggeredStencilDerivation("NE", 3, (0,)).apply(c3.center) + \
FiniteDifferenceStaggeredStencilDerivation("NE", 3, (1,)).apply(c3.center) == c3[1, 1, 0] - c3[0, 0, 0]
```
%% Cell type:code id: tags:
``` python
d = FiniteDifferenceStaggeredStencilDerivation("NE", 2, (0, 1))
assert d.apply(c.center) == c[0,0] + c[1,1] - c[1,0] - c[0,1]
d.visualize()
```
%% Output
%% Cell type:code id: tags:
``` python
v3 = ps.fields("v(3): [3D]")
for i in range(*v3.index_shape):
assert FiniteDifferenceStaggeredStencilDerivation("E", 3, (0,)).apply(v3.center_vector[i]) == \
v3[1,0,0](i) - v3[0,0,0](i)
```
......
pystencils_tests/test_field.py → tests/test_field.py
View file @ 5600b6b6
......@@ -4,22 +4,34 @@ import sympy as sp
import pystencils as ps
from pystencils import TypedSymbol
from pystencils.data_types import create_type
from pystencils.field import Field, FieldType, layout_string_to_tuple
from pystencils.typing import create_type
from pystencils.field import Field, FieldType, layout_string_to_tuple, spatial_layout_string_to_tuple
def test_field_basic():
f = Field.create_generic('f', spatial_dimensions=2)
f = Field.create_generic("f", spatial_dimensions=2)
assert FieldType.is_generic(f)
assert f['E'] == f[1, 0]
assert f['N'] == f[0, 1]
assert '_' in f.center._latex('dummy')
assert f.index_to_physical(index_coordinates=sp.Matrix([0, 0]), staggered=False)[0] == 0
assert f.index_to_physical(index_coordinates=sp.Matrix([0, 0]), staggered=False)[1] == 0
assert f.physical_to_index(physical_coordinates=sp.Matrix([0, 0]), staggered=False)[0] == 0
assert f.physical_to_index(physical_coordinates=sp.Matrix([0, 0]), staggered=False)[1] == 0
assert f["E"] == f[1, 0]
assert f["N"] == f[0, 1]
assert "_" in f.center._latex("dummy")
assert (
f.index_to_physical(index_coordinates=sp.Matrix([0, 0]), staggered=False)[0]
== 0
)
assert (
f.index_to_physical(index_coordinates=sp.Matrix([0, 0]), staggered=False)[1]
== 0
)
assert (
f.physical_to_index(physical_coordinates=sp.Matrix([0, 0]), staggered=False)[0]
== 0
)
assert (
f.physical_to_index(physical_coordinates=sp.Matrix([0, 0]), staggered=False)[1]
== 0
)
f1 = f.new_field_with_different_name("f1")
assert f1.ndim == f.ndim
......@@ -28,7 +40,7 @@ def test_field_basic():
fixed = ps.fields("f(5, 5) : double[20, 20]")
assert fixed.neighbor_vector((1, 1)).shape == (5, 5)
f = Field.create_fixed_size('f', (10, 10), strides=(80, 8), dtype=np.float64)
f = Field.create_fixed_size("f", (10, 10), strides=(80, 8), dtype=np.float64)
assert f.spatial_strides == (10, 1)
assert f.index_strides == ()
assert f.center_vector == sp.Matrix([f.center])
......@@ -37,20 +49,21 @@ def test_field_basic():
assert f1.ndim == f.ndim
assert f1.values_per_cell() == f.values_per_cell()
f = Field.create_fixed_size('f', (8, 8, 2, 2), index_dimensions=2)
assert f.center_vector == sp.Matrix([[f(0, 0), f(0, 1)],
[f(1, 0), f(1, 1)]])
f = Field.create_fixed_size("f", (8, 8, 2, 2), index_dimensions=2)
assert f.center_vector == sp.Matrix([[f(0, 0), f(0, 1)], [f(1, 0), f(1, 1)]])
field_access = f[1, 1]
assert field_access.nr_of_coordinates == 2
assert field_access.offset_name == 'NE'
assert field_access.offset_name == "NE"
neighbor = field_access.neighbor(coord_id=0, offset=-2)
assert neighbor.offsets == (-1, 1)
assert '_' in neighbor._latex('dummy')
assert "_" in neighbor._latex("dummy")
f = Field.create_fixed_size('f', (8, 8, 2, 2, 2), index_dimensions=3)
assert f.center_vector == sp.Matrix([[[f(i, j, k) for k in range(2)] for j in range(2)] for i in range(2)])
f = Field.create_fixed_size("f", (8, 8, 2, 2, 2), index_dimensions=3)
assert f.center_vector == sp.Array(
[[[f(i, j, k) for k in range(2)] for j in range(2)] for i in range(2)]
)
f = Field.create_generic('f', spatial_dimensions=5, index_dimensions=2)
f = Field.create_generic("f", spatial_dimensions=5, index_dimensions=2)
field_access = f[1, -1, 2, -3, 0](1, 0)
assert field_access.offsets == (1, -1, 2, -3, 0)
assert field_access.index == (1, 0)
......@@ -60,61 +73,71 @@ def test_error_handling():
struct_dtype = np.dtype([('a', np.int32), ('b', np.float64), ('c', np.uint32)])
Field.create_generic('f', spatial_dimensions=2, index_dimensions=0, dtype=struct_dtype)
with pytest.raises(ValueError) as e:
Field.create_generic('f', spatial_dimensions=2, index_dimensions=1, dtype=struct_dtype)
assert 'index dimension' in str(e.value)
Field.create_generic(
"f", spatial_dimensions=2, index_dimensions=1, dtype=struct_dtype
)
assert "index dimension" in str(e.value)
arr = np.array([[[(1,)*3, (2,)*3, (3,)*3]]*2], dtype=struct_dtype)
Field.create_from_numpy_array('f', arr, index_dimensions=0)
arr = np.array([[[(1,) * 3, (2,) * 3, (3,) * 3]] * 2], dtype=struct_dtype)
Field.create_from_numpy_array("f", arr, index_dimensions=0)
with pytest.raises(ValueError) as e:
Field.create_from_numpy_array('f', arr, index_dimensions=1)
assert 'Structured arrays' in str(e.value)
Field.create_from_numpy_array("f", arr, index_dimensions=1)
assert "Structured arrays" in str(e.value)
arr = np.zeros([3, 3, 3])
Field.create_from_numpy_array('f', arr, index_dimensions=2)
Field.create_from_numpy_array("f", arr, index_dimensions=2)
with pytest.raises(ValueError) as e:
Field.create_from_numpy_array('f', arr, index_dimensions=3)
assert 'Too many' in str(e.value)
Field.create_from_numpy_array("f", arr, index_dimensions=3)
assert "Too many" in str(e.value)
Field.create_fixed_size('f', (3, 2, 4), index_dimensions=0, dtype=struct_dtype, layout='reverse_numpy')
Field.create_fixed_size(
"f", (3, 2, 4), index_dimensions=0, dtype=struct_dtype, layout="reverse_numpy"
)
with pytest.raises(ValueError) as e:
Field.create_fixed_size('f', (3, 2, 4), index_dimensions=1, dtype=struct_dtype, layout='reverse_numpy')
assert 'Structured arrays' in str(e.value)
f = Field.create_fixed_size('f', (10, 10))
Field.create_fixed_size(
"f",
(3, 2, 4),
index_dimensions=1,
dtype=struct_dtype,
layout="reverse_numpy",
)
assert "Structured arrays" in str(e.value)
f = Field.create_fixed_size("f", (10, 10))
with pytest.raises(ValueError) as e:
f[1]
assert 'Wrong number of spatial indices' in str(e.value)
assert "Wrong number of spatial indices" in str(e.value)
f = Field.create_generic('f', spatial_dimensions=2, index_shape=(3,))
f = Field.create_generic("f", spatial_dimensions=2, index_shape=(3,))
with pytest.raises(ValueError) as e:
f(3)
assert 'out of bounds' in str(e.value)
assert "out of bounds" in str(e.value)
f = Field.create_fixed_size('f', (10, 10, 3, 4), index_dimensions=2)
f = Field.create_fixed_size("f", (10, 10, 3, 4), index_dimensions=2)
with pytest.raises(ValueError) as e:
f(3, 0)
assert 'out of bounds' in str(e.value)
assert "out of bounds" in str(e.value)
with pytest.raises(ValueError) as e:
f(1, 0)(1, 0)
assert 'Indexing an already indexed' in str(e.value)
assert "Indexing an already indexed" in str(e.value)
with pytest.raises(ValueError) as e:
f(1)
assert 'Wrong number of indices' in str(e.value)
assert "Wrong number of indices" in str(e.value)
with pytest.raises(ValueError) as e:
Field.create_generic('f', spatial_dimensions=2, layout='wrong')
assert 'Unknown layout descriptor' in str(e.value)
Field.create_generic("f", spatial_dimensions=2, layout="wrong")
assert "Unknown layout descriptor" in str(e.value)
assert layout_string_to_tuple('fzyx', dim=4) == (3, 2, 1, 0)
assert layout_string_to_tuple("fzyx", dim=4) == (3, 2, 1, 0)
with pytest.raises(ValueError) as e:
layout_string_to_tuple('wrong', dim=4)
assert 'Unknown layout descriptor' in str(e.value)
layout_string_to_tuple("wrong", dim=4)
assert "Unknown layout descriptor" in str(e.value)
def test_decorator_scoping():
dst = ps.fields('dst : double[2D]')
dst = ps.fields("dst : double[2D]")
def f1():
a = sp.Symbol("a")
......@@ -134,7 +157,7 @@ def test_decorator_scoping():
def test_string_creation():
x, y, z = ps.fields(' x(4), y(3,5) z : double[ 3, 47]')
x, y, z = ps.fields(" x(4), y(3,5) z : double[ 3, 47]")
assert x.index_shape == (4,)
assert y.index_shape == (3, 5)
assert z.spatial_shape == (3, 47)
......@@ -142,19 +165,85 @@ def test_string_creation():
def test_itemsize():
x = ps.fields('x: float32[1d]')
y = ps.fields('y: float64[2d]')
i = ps.fields('i: int16[1d]')
x = ps.fields("x: float32[1d]")
y = ps.fields("y: float64[2d]")
i = ps.fields("i: int16[1d]")
assert x.itemsize == 4
assert y.itemsize == 8
assert i.itemsize == 2
def test_spatial_memory_layout_descriptors():
assert (
spatial_layout_string_to_tuple("AoS", 3)
== spatial_layout_string_to_tuple("aos", 3)
== spatial_layout_string_to_tuple("ZYXF", 3)
== spatial_layout_string_to_tuple("zyxf", 3)
== (2, 1, 0)
)
assert (
spatial_layout_string_to_tuple("SoA", 3)
== spatial_layout_string_to_tuple("soa", 3)
== spatial_layout_string_to_tuple("FZYX", 3)
== spatial_layout_string_to_tuple("fzyx", 3)
== spatial_layout_string_to_tuple("f", 3)
== spatial_layout_string_to_tuple("F", 3)
== (2, 1, 0)
)
assert (
spatial_layout_string_to_tuple("c", 3)
== spatial_layout_string_to_tuple("C", 3)
== (0, 1, 2)
)
assert spatial_layout_string_to_tuple("C", 5) == (0, 1, 2, 3, 4)
with pytest.raises(ValueError):
spatial_layout_string_to_tuple("aos", -1)
with pytest.raises(ValueError):
spatial_layout_string_to_tuple("aos", 4)
def test_memory_layout_descriptors():
assert (
layout_string_to_tuple("AoS", 4)
== layout_string_to_tuple("aos", 4)
== layout_string_to_tuple("ZYXF", 4)
== layout_string_to_tuple("zyxf", 4)
== (2, 1, 0, 3)
)
assert (
layout_string_to_tuple("SoA", 4)
== layout_string_to_tuple("soa", 4)
== layout_string_to_tuple("FZYX", 4)
== layout_string_to_tuple("fzyx", 4)
== layout_string_to_tuple("f", 4)
== layout_string_to_tuple("F", 4)
== (3, 2, 1, 0)
)
assert (
layout_string_to_tuple("c", 4)
== layout_string_to_tuple("C", 4)
== (0, 1, 2, 3)
)
assert layout_string_to_tuple("C", 5) == (0, 1, 2, 3, 4)
with pytest.raises(ValueError):
layout_string_to_tuple("aos", -1)
with pytest.raises(ValueError):
layout_string_to_tuple("aos", 5)
def test_staggered():
# D2Q5
j1, j2, j3 = ps.fields('j1(2), j2(2,2), j3(2,2,2) : double[2D]', field_type=FieldType.STAGGERED)
j1, j2, j3 = ps.fields(
"j1(2), j2(2,2), j3(2,2,2) : double[2D]", field_type=FieldType.STAGGERED
)
assert j1[0, 1](1) == j1.staggered_access((0, sp.Rational(1, 2)))
assert j1[0, 1](1) == j1.staggered_access(np.array((0, sp.Rational(1, 2))))
......@@ -163,7 +252,7 @@ def test_staggered():
assert j1[0, 1](1) == j1.staggered_access("N")
assert j1[0, 0](1) == j1.staggered_access("S")
assert j1.staggered_vector_access("N") == sp.Matrix([j1.staggered_access("N")])
assert j1.staggered_stencil_name == 'D2Q5'
assert j1.staggered_stencil_name == "D2Q5"
assert j1.physical_coordinates[0] == TypedSymbol("ctr_0", create_type("int"), nonnegative=True)
assert j1.physical_coordinates[1] == TypedSymbol("ctr_1", create_type("int"), nonnegative=True)
......@@ -176,28 +265,40 @@ def test_staggered():
assert j2[0, 1](1, 1) == j2.staggered_access((0, sp.Rational(1, 2)), 1)
assert j2[0, 1](1, 1) == j2.staggered_access("N", 1)
assert j2.staggered_vector_access("N") == sp.Matrix([j2.staggered_access("N", 0), j2.staggered_access("N", 1)])
assert j2.staggered_vector_access("N") == sp.Matrix(
[j2.staggered_access("N", 0), j2.staggered_access("N", 1)]
)
assert j3[0, 1](1, 1, 1) == j3.staggered_access((0, sp.Rational(1, 2)), (1, 1))
assert j3[0, 1](1, 1, 1) == j3.staggered_access("N", (1, 1))
assert j3.staggered_vector_access("N") == sp.Matrix([[j3.staggered_access("N", (i, j))
for j in range(2)] for i in range(2)])
assert j3.staggered_vector_access("N") == sp.Matrix(
[[j3.staggered_access("N", (i, j)) for j in range(2)] for i in range(2)]
)
# D2Q9
k1, k2 = ps.fields('k1(4), k2(2) : double[2D]', field_type=FieldType.STAGGERED)
k1, k2 = ps.fields("k1(4), k2(2) : double[2D]", field_type=FieldType.STAGGERED)
assert k1[1, 1](2) == k1.staggered_access("NE")
assert k1[0, 0](2) == k1.staggered_access("SW")
assert k1[0, 0](3) == k1.staggered_access("NW")
a = k1.staggered_access("NE")
assert a._staggered_offset(a.offsets, a.index[0]) == [sp.Rational(1, 2), sp.Rational(1, 2)]
assert a._staggered_offset(a.offsets, a.index[0]) == [
sp.Rational(1, 2),
sp.Rational(1, 2),
]
a = k1.staggered_access("SW")
assert a._staggered_offset(a.offsets, a.index[0]) == [sp.Rational(-1, 2), sp.Rational(-1, 2)]
assert a._staggered_offset(a.offsets, a.index[0]) == [
sp.Rational(-1, 2),
sp.Rational(-1, 2),
]
a = k1.staggered_access("NW")
assert a._staggered_offset(a.offsets, a.index[0]) == [sp.Rational(-1, 2), sp.Rational(1, 2)]
assert a._staggered_offset(a.offsets, a.index[0]) == [
sp.Rational(-1, 2),
sp.Rational(1, 2),
]
# sign reversed when using as flux field
r = ps.fields('r(2) : double[2D]', field_type=FieldType.STAGGERED_FLUX)
r = ps.fields("r(2) : double[2D]", field_type=FieldType.STAGGERED_FLUX)
assert r[0, 0](0) == r.staggered_access("W")
assert -r[1, 0](0) == r.staggered_access("E")
pystencils_tests/test_field_access_poly.py → tests/test_field_access_poly.py
View file @ 5600b6b6
......@@ -7,7 +7,7 @@
"""
import pytest
from pystencils.session import *
from sympy import poly
......@@ -22,8 +22,14 @@ def test_field_access_poly():
def test_field_access_piecewise():
dh = ps.create_data_handling((20, 20))
ρ = dh.add_array('rho')
pw = sp.Piecewise((0, 1 < sp.Max(-0.5, ρ.center+0.5)), (1, True))
a = sp.simplify(pw)
print(a)
try:
a = sp.Piecewise((0, 1 < sp.Max(-0.5, sp.Symbol("test") + 0.5)), (1, True))
a.simplify()
except Exception as e:
pytest.skip(f"Bug in SymPy 1.10: {e}")
else:
dh = ps.create_data_handling((20, 20))
ρ = dh.add_array('rho')
pw = sp.Piecewise((0, 1 < sp.Max(-0.5, ρ.center+0.5)), (1, True))
a = sp.simplify(pw)
print(a)
pystencils_tests/test_field_equality.ipynb → tests/test_field_equality.ipynb
View file @ 5600b6b6
%% Cell type:code id: tags:
``` python
from pystencils.session import *
from pystencils.data_types import cast_func
```
%% Cell type:markdown id: tags:
## Test field equality behaviour
%% Cell type:markdown id: tags:
Fields create with same parameters are equal
%% Cell type:code id: tags:
``` python
f1 = ps.Field.create_generic('f', spatial_dimensions=2, index_dimensions=0)
f2 = ps.Field.create_generic('f', spatial_dimensions=2, index_dimensions=0)
assert f1 == f2
```
%% Cell type:code id: tags:
``` python
print("Field ids equal in accesses: ", id(f1.center._field) == id(f2.center._field))
print("Field accesses equal: ", f1.center == f2.center)
```
%% Output
Field ids equal in accesses: True
Field accesses equal: True
%% Cell type:code id: tags:
``` python
f1 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0)
f2 = ps.Field.create_generic('f', spatial_dimensions=2, index_dimensions=0)
assert f1 != f2
```
%% Cell type:code id: tags:
``` python
f1 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0)
f2 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0, dtype=np.float32)
assert f1 != f2
```
%% Cell type:markdown id: tags:
Properties of fields:
- `field_type`: enum distinguishing normal from index or buffer fields
- `_dtype`: data type of field elements
- `_layout`: tuple indicating the memory linearization order
- `shape`: size of field for each dimension
- `strides`: number of elements to jump over to increase coordinate of this dimension by one
- `latex_name`: optional display name when field is printed as latex
Equality compare of fields:
- field has `__eq__` and ``__hash__`` overridden
- all parameter but `latex_name` are considered for equality
%% Cell type:markdown id: tags:
## Test field access equality behaviour
%% Cell type:code id: tags:
``` python
f1 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0)
f2 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0)
assert f1.center == f2.center
```
%% Cell type:code id: tags:
``` python
f1 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0)
f2 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0, dtype=np.float32)
assert f1.center != f2.center
```
%% Cell type:code id: tags:
``` python
def print_field_accesses_debug(expr):
from pystencils import Field
fas = list(expr.atoms(Field.Access))
fields = list({e.field for e in fas})
print("Field Accesses:")
for fa in fas:
print(f" - {fa}, hash {hash(fa)}, offsets {fa._offsets}, index {fa._index}, {fa._hashable_content()}")
print("")
for i in range(len(fas)):
for j in range(len(fas)):
if not i < j:
continue
print( f" -> {i},{j} {fas[i]} == {fas[j]}: {fas[i] == {fas[j]}}")
print("Fields")
for f in fields:
print(f" - {f}, {id(f)}, shape {f.shape}, strides {f.strides}, {f._dtype}, {f.field_type}, layout {f.layout}")
print("")
for i in range(len(fields)):
for j in range(len(fields)):
if not i < j:
continue
print(f" - {fields[i]} == {fields[j]}: {fields[i] == fields[j]}, ids equal {id(fields[i])==id(fields[j])}, hash equal {hash(fields[i])==hash(fields[j])}")
```
%% Cell type:code id: tags:
``` python
print_field_accesses_debug(f1.center * f2.center)
```
%% Output
Field Accesses:
- f[0], hash -3276894289571194847, offsets (0,), index (), (('f_C', ('commutative', True)), ((0,), (_size_f_0,), (_stride_f_0,), 3146377891102027609, <FieldType.GENERIC: 0>, 'f', None), 0)
- f[0], hash -1516451775709390846, offsets (0,), index (), (('f_C', ('commutative', True)), ((0,), (_size_f_0,), (_stride_f_0,), -1421177580377734245, <FieldType.GENERIC: 0>, 'f', None), 0)
- f[0], hash -8859424145258271267, offsets (0,), index (), ((('f_C', ('commutative', True), ('complex', True), ('extended_real', True), ('finite', True), ('hermitian', True), ('imaginary', False), ('infinite', False), ('real', True)), 2305067722319023373), ((0,), (_size_f_0,), (_stride_f_0,), <FieldType.GENERIC: 0>, 'f', None, double), 0)
- f[0], hash -6454673863007224785, offsets (0,), index (), ((('f_C', ('commutative', True), ('complex', True), ('extended_real', True), ('finite', True), ('hermitian', True), ('imaginary', False), ('infinite', False), ('real', True)), 4093629613697528859), ((0,), (_size_f_0,), (_stride_f_0,), <FieldType.GENERIC: 0>, 'f', None, float), 0)
-> 0,1 f[0] == f[0]: False
Fields
- f, 140548694371968, shape (_size_f_0,), strides (_stride_f_0,), double, FieldType.GENERIC, layout (0,)
- f, 140548693963104, shape (_size_f_0,), strides (_stride_f_0,), float, FieldType.GENERIC, layout (0,)
- f, 4881406800, shape (_size_f_0,), strides (_stride_f_0,), double, FieldType.GENERIC, layout (0,)
- f, 4881445024, shape (_size_f_0,), strides (_stride_f_0,), float, FieldType.GENERIC, layout (0,)
- f == f: False, ids equal False, hash equal False
......
pystencils_tests/test_floor_ceil_int_optimization.py → tests/test_floor_ceil_int_optimization.py
View file @ 5600b6b6
......@@ -11,7 +11,7 @@
import sympy as sp
import pystencils
from pystencils.data_types import create_type
from pystencils.typing import create_type
def test_floor_ceil_int_optimization():
......
pystencils_tests/test_fvm.py → tests/test_fvm.py
View file @ 5600b6b6
......@@ -3,6 +3,9 @@ import pystencils as ps
import numpy as np
import pytest
from itertools import product
from pystencils.rng import random_symbol
from pystencils.astnodes import SympyAssignment
from pystencils.node_collection import NodeCollection
def advection_diffusion(dim: int):
......@@ -12,7 +15,7 @@ def advection_diffusion(dim: int):
elif dim == 3:
L = (16, 16, 16)
dh = ps.create_data_handling(domain_size=L, periodicity=True, default_target='cpu')
dh = ps.create_data_handling(domain_size=L, periodicity=True, default_target=ps.Target.CPU)
n_field = dh.add_array('n', values_per_cell=1)
j_field = dh.add_array('j', values_per_cell=3 ** dim // 2, field_type=ps.FieldType.STAGGERED_FLUX)
......@@ -118,12 +121,307 @@ def test_advection_diffusion_2d(velocity):
@pytest.mark.parametrize("velocity", list(product([0, -0.047, 0.041], [0, -0.031, 0.023], [0, -0.017, 0.011])))
@pytest.mark.longrun
def test_advection_diffusion_3d(velocity):
if 3 not in advection_diffusion.runners:
advection_diffusion.runners[3] = advection_diffusion(3)
advection_diffusion.runners[3](velocity)
def advection_diffusion_fluctuations(dim: int):
# parameters
if dim == 2:
L = (32, 32)
stencil_factor = np.sqrt(1 / (1 + np.sqrt(2)))
elif dim == 3:
L = (16, 16, 16)
stencil_factor = np.sqrt(1 / (1 + 2 * np.sqrt(2) + 4.0 / 3.0 * np.sqrt(3)))
dh = ps.create_data_handling(domain_size=L, periodicity=True, default_target=ps.Target.CPU)
n_field = dh.add_array('n', values_per_cell=1)
j_field = dh.add_array('j', values_per_cell=3 ** dim // 2, field_type=ps.FieldType.STAGGERED_FLUX)
velocity_field = dh.add_array('v', values_per_cell=dim)
D = 0.00666
time = 10000
def grad(f):
return sp.Matrix([ps.fd.diff(f, i) for i in range(dim)])
flux_eq = - D * grad(n_field)
fvm_eq = ps.fd.FVM1stOrder(n_field, flux=flux_eq)
vof_adv = ps.fd.VOF(j_field, velocity_field, n_field)
# merge calculation of advection and diffusion terms
flux = []
for adv, div in zip(vof_adv, fvm_eq.discrete_flux(j_field)):
assert adv.lhs == div.lhs
flux.append(ps.Assignment(adv.lhs, adv.rhs + div.rhs))
flux = ps.AssignmentCollection(flux)
rng_symbol_gen = random_symbol(flux.subexpressions, dim=dh.dim)
for i in range(len(flux.main_assignments)):
n = j_field.staggered_stencil[i]
assert flux.main_assignments[i].lhs == j_field.staggered_access(n)
# calculate mean density
dens = (n_field.neighbor_vector(n) + n_field.center_vector)[0] / 2
# multyply by smoothed haviside function so that fluctuation will not get bigger that the density
dens *= sp.Max(0, sp.Min(1.0, n_field.neighbor_vector(n)[0]) * sp.Min(1.0, n_field.center_vector[0]))
# lenght of the vector
length = sp.sqrt(len(j_field.staggered_stencil[i]))
# amplitude of the random fluctuations
fluct = sp.sqrt(2 * dens * D) * sp.sqrt(1 / length) * stencil_factor
# add fluctuations
fluct *= 2 * (next(rng_symbol_gen) - 0.5) * sp.sqrt(3)
flux.main_assignments[i] = ps.Assignment(flux.main_assignments[i].lhs, flux.main_assignments[i].rhs + fluct)
# Add the folding to the flux, so that the random numbers persist through the ghostlayers.
fold = {ps.astnodes.LoopOverCoordinate.get_loop_counter_symbol(i):
ps.astnodes.LoopOverCoordinate.get_loop_counter_symbol(i) % L[i] for i in range(len(L))}
flux.subs(fold)
flux_kernel = ps.create_staggered_kernel(flux).compile()
pde_kernel = ps.create_kernel(fvm_eq.discrete_continuity(j_field)).compile()
sync_conc = dh.synchronization_function([n_field.name])
# analytical density distribution calculation
def P(rho, density_init):
res = []
for r in rho:
res.append(np.power(density_init, r) * np.exp(-density_init) / np.math.gamma(r + 1))
return np.array(res)
def run(density_init: float, velocity: np.ndarray, time: int):
dh.fill(n_field.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
dh.fill(j_field.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
# set initial values for velocity and density
for i in range(dim):
dh.fill(velocity_field.name, velocity[i], i, ghost_layers=True, inner_ghost_layers=True)
dh.fill(n_field.name, density_init)
measurement_intervall = 10
warm_up = 1000
data = []
sync_conc()
for i in range(warm_up):
dh.run_kernel(flux_kernel, seed=42, time_step=i)
dh.run_kernel(pde_kernel)
sync_conc()
for i in range(time):
dh.run_kernel(flux_kernel, seed=42, time_step=i + warm_up)
dh.run_kernel(pde_kernel)
sync_conc()
if(i % measurement_intervall == 0):
data = np.append(data, dh.gather_array(n_field.name).ravel(), 0)
# test mass conservation
np.testing.assert_almost_equal(dh.gather_array(n_field.name).mean(), density_init)
n_bins = 50
density_value, bins = np.histogram(data, density=True, bins=n_bins)
bins_mean = bins[:-1] + (bins[1:] - bins[:-1]) / 2
analytical_value = P(bins_mean, density_init)
print(density_value - analytical_value)
np.testing.assert_allclose(density_value, analytical_value, atol=2e-3)
return lambda density_init, v: run(density_init, np.array(v), time)
advection_diffusion_fluctuations.runners = {}
@pytest.mark.parametrize("velocity", list(product([0, 0.00041], [0, -0.00031])))
@pytest.mark.parametrize("density", [27.0, 56.5])
@pytest.mark.longrun
def test_advection_diffusion_fluctuation_2d(density, velocity):
if 2 not in advection_diffusion_fluctuations.runners:
advection_diffusion_fluctuations.runners[2] = advection_diffusion_fluctuations(2)
advection_diffusion_fluctuations.runners[2](density, velocity)
@pytest.mark.parametrize("velocity", [(0.0, 0.0, 0.0), (0.00043, -0.00017, 0.00028)])
@pytest.mark.parametrize("density", [27.0, 56.5])
@pytest.mark.longrun
def test_advection_diffusion_fluctuation_3d(density, velocity):
if 3 not in advection_diffusion_fluctuations.runners:
advection_diffusion_fluctuations.runners[3] = advection_diffusion_fluctuations(3)
advection_diffusion_fluctuations.runners[3](density, velocity)
def diffusion_reaction(fluctuations: bool):
# parameters
L = (32, 32)
stencil_factor = np.sqrt(1 / (1 + np.sqrt(2)))
dh = ps.create_data_handling(domain_size=L, periodicity=True, default_target=ps.Target.CPU)
species = 2
n_fields = []
j_fields = []
r_flux_fields = []
for i in range(species):
n_fields.append(dh.add_array(f'n_{i}', values_per_cell=1))
j_fields.append(dh.add_array(f'j_{i}', values_per_cell=3 ** dh.dim // 2,
field_type=ps.FieldType.STAGGERED_FLUX))
r_flux_fields.append(dh.add_array(f'r_{i}', values_per_cell=1))
velocity_field = dh.add_array('v', values_per_cell=dh.dim)
D = 0.00666
time = 1000
r_order = [2.0, 0.0]
r_rate_const = 0.00001
r_coefs = [-2, 1]
def grad(f):
return sp.Matrix([ps.fd.diff(f, i) for i in range(dh.dim)])
flux_eq = - D * grad(n_fields[0])
fvm_eq = ps.fd.FVM1stOrder(n_fields[0], flux=flux_eq)
vof_adv = ps.fd.VOF(j_fields[0], velocity_field, n_fields[0])
continuity_assignments = fvm_eq.discrete_continuity(j_fields[0])
# merge calculation of advection and diffusion terms
flux = []
for adv, div in zip(vof_adv, fvm_eq.discrete_flux(j_fields[0])):
assert adv.lhs == div.lhs
flux.append(ps.Assignment(adv.lhs, adv.rhs + div.rhs))
flux = ps.AssignmentCollection(flux)
if(fluctuations):
rng_symbol_gen = random_symbol(flux.subexpressions, dim=dh.dim)
for i in range(len(flux.main_assignments)):
n = j_fields[0].staggered_stencil[i]
assert flux.main_assignments[i].lhs == j_fields[0].staggered_access(n)
# calculate mean density
dens = (n_fields[0].neighbor_vector(n) + n_fields[0].center_vector)[0] / 2
# multyply by smoothed haviside function so that fluctuation will not get bigger that the density
dens *= sp.Max(0,
sp.Min(1.0, n_fields[0].neighbor_vector(n)[0]) * sp.Min(1.0, n_fields[0].center_vector[0]))
# lenght of the vector
length = sp.sqrt(len(j_fields[0].staggered_stencil[i]))
# amplitude of the random fluctuations
fluct = sp.sqrt(2 * dens * D) * sp.sqrt(1 / length) * stencil_factor
# add fluctuations
fluct *= 2 * (next(rng_symbol_gen) - 0.5) * sp.sqrt(3)
flux.main_assignments[i] = ps.Assignment(flux.main_assignments[i].lhs, flux.main_assignments[i].rhs + fluct)
# Add the folding to the flux, so that the random numbers persist through the ghostlayers.
fold = {ps.astnodes.LoopOverCoordinate.get_loop_counter_symbol(i):
ps.astnodes.LoopOverCoordinate.get_loop_counter_symbol(i) % L[i] for i in range(len(L))}
flux.subs(fold)
r_flux = NodeCollection([SympyAssignment(j_fields[i].center, 0) for i in range(species)])
reaction = r_rate_const
for i in range(species):
reaction *= sp.Pow(n_fields[i].center, r_order[i])
new_assignments = []
if fluctuations:
rng_symbol_gen = random_symbol(new_assignments, dim=dh.dim)
reaction_fluctuations = sp.sqrt(sp.Abs(reaction)) * 2 * (next(rng_symbol_gen) - 0.5) * sp.sqrt(3)
reaction_fluctuations *= sp.Min(1, sp.Abs(reaction**2))
else:
reaction_fluctuations = 0.0
for i in range(species):
r_flux.all_assignments[i] = SympyAssignment(
r_flux_fields[i].center, (reaction + reaction_fluctuations) * r_coefs[i])
[r_flux.all_assignments.insert(0, new) for new in new_assignments]
continuity_assignments = [SympyAssignment(*assignment.args) for assignment in continuity_assignments]
continuity_assignments.append(SympyAssignment(n_fields[0].center, n_fields[0].center + r_flux_fields[0].center))
flux_kernel = ps.create_staggered_kernel(flux).compile()
reaction_kernel = ps.create_kernel(r_flux).compile()
config = ps.CreateKernelConfig(allow_double_writes=True)
pde_kernel = ps.create_kernel(continuity_assignments, config=config).compile()
sync_conc = dh.synchronization_function([n_fields[0].name, n_fields[1].name])
def f(t, r, n0, fac, fluctuations):
"""Calculates the amount of product created after a certain time of a reaction with form xA -> B
Args:
t: Time of the reation
r: Reaction rate constant
n0: Initial density of the
fac: Reaction order of A (this in most cases equals the stochometric coefficient x)
fluctuations: Boolian whether fluctuations were included during the reaction.
"""
if fluctuations:
return 1 / fac * (n0 + n0 / (n0 - (n0 + 1) * np.exp(fac * r * t)))
return 1 / fac * (n0 - (1 / (fac * r * t + (1 / n0))))
def run(density_init: float, velocity: np.ndarray, time: int):
for i in range(species):
dh.fill(n_fields[i].name, np.nan, ghost_layers=True, inner_ghost_layers=True)
dh.fill(j_fields[i].name, 0.0, ghost_layers=True, inner_ghost_layers=True)
dh.fill(r_flux_fields[i].name, 0.0, ghost_layers=True, inner_ghost_layers=True)
# set initial values for velocity and density
for i in range(dh.dim):
dh.fill(velocity_field.name, velocity[i], i, ghost_layers=True, inner_ghost_layers=True)
dh.fill(n_fields[0].name, density_init)
dh.fill(n_fields[1].name, 0.0)
measurement_intervall = 10
data = []
sync_conc()
for i in range(time):
if(i % measurement_intervall == 0):
data.append([i, dh.gather_array(n_fields[1].name).mean(), dh.gather_array(n_fields[0].name).mean()])
dh.run_kernel(reaction_kernel, seed=41, time_step=i)
for s_idx in range(species):
flux_kernel(n_0=dh.cpu_arrays[n_fields[s_idx].name],
j_0=dh.cpu_arrays[j_fields[s_idx].name],
v=dh.cpu_arrays[velocity_field.name], seed=42 + s_idx, time_step=i)
pde_kernel(n_0=dh.cpu_arrays[n_fields[s_idx].name],
j_0=dh.cpu_arrays[j_fields[s_idx].name],
r_0=dh.cpu_arrays[r_flux_fields[s_idx].name])
sync_conc()
data = np.array(data).transpose()
x = data[0]
analytical_value = f(x, r_rate_const, density_init, abs(r_coefs[0]), fluctuations)
# test mass conservation
np.testing.assert_almost_equal(
dh.gather_array(n_fields[0].name).mean() + 2 * dh.gather_array(n_fields[1].name).mean(), density_init)
r_tol = 2e-3
if fluctuations:
r_tol = 3e-2
np.testing.assert_allclose(data[1], analytical_value, rtol=r_tol)
return lambda density_init, v: run(density_init, np.array(v), time)
advection_diffusion_fluctuations.runners = {}
@pytest.mark.parametrize("velocity", list(product([0, 0.0041], [0, -0.0031])))
@pytest.mark.parametrize("density", [27.0, 56.5])
@pytest.mark.parametrize("fluctuations", [False, True])
@pytest.mark.longrun
def test_diffusion_reaction(fluctuations, density, velocity):
diffusion_reaction.runner = diffusion_reaction(fluctuations)
diffusion_reaction.runner(density, velocity)
def VOF2(j: ps.field.Field, v: ps.field.Field, ρ: ps.field.Field, simplify=True):
"""Volume-of-fluid discretization of advection
......@@ -203,7 +501,7 @@ def VOF2(j: ps.field.Field, v: ps.field.Field, ρ: ps.field.Field, simplify=True
@pytest.mark.parametrize("dim", [2, 3])
def test_advection(dim):
L = (8,) * dim
dh = ps.create_data_handling(L, periodicity=True, default_target='cpu')
dh = ps.create_data_handling(L, periodicity=True, default_target=ps.Target.CPU)
c = dh.add_array('c', values_per_cell=1)
j = dh.add_array('j', values_per_cell=3 ** dh.dim // 2, field_type=ps.FieldType.STAGGERED_FLUX)
u = dh.add_array('u', values_per_cell=dh.dim)
......@@ -235,7 +533,7 @@ def test_ek(stencil):
# data structures
dh = ps.create_data_handling(L, periodicity=True, default_target='cpu')
dh = ps.create_data_handling(L, periodicity=True, default_target=ps.Target.CPU)
c = dh.add_array('c', values_per_cell=1)
j = dh.add_array('j', values_per_cell=int(stencil[-1]) // 2, field_type=ps.FieldType.STAGGERED_FLUX)
Phi = dh.add_array('Φ', values_per_cell=1)
......@@ -288,7 +586,7 @@ def test_ek(stencil):
@pytest.mark.parametrize("derivative", [0, 1])
def test_flux_stencil(stencil, derivative):
L = (40, ) * int(stencil[1])
dh = ps.create_data_handling(L, periodicity=True, default_target='cpu')
dh = ps.create_data_handling(L, periodicity=True, default_target=ps.Target.CPU)
c = dh.add_array('c', values_per_cell=1)
j = dh.add_array('j', values_per_cell=int(stencil[3:]) // 2, field_type=ps.FieldType.STAGGERED_FLUX)
......@@ -311,7 +609,7 @@ def test_flux_stencil(stencil, derivative):
@pytest.mark.parametrize("stencil", ["D2Q5", "D2Q9", "D3Q7", "D3Q19", "D3Q27"])
def test_source_stencil(stencil):
L = (40, ) * int(stencil[1])
dh = ps.create_data_handling(L, periodicity=True, default_target='cpu')
dh = ps.create_data_handling(L, periodicity=True, default_target=ps.Target.CPU)
c = dh.add_array('c', values_per_cell=1)
j = dh.add_array('j', values_per_cell=int(stencil[3:]) // 2, field_type=ps.FieldType.STAGGERED_FLUX)
......@@ -324,3 +622,20 @@ def test_source_stencil(stencil):
assert len(diff.atoms(ps.field.Field.Access)) == 1
else:
assert len(diff.atoms(ps.field.Field.Access)) == 2
def test_fvm_staggered_simplification():
D = sp.Symbol("D")
data_type = "float64"
c = ps.fields(f"c: {data_type}[2D]", layout='fzyx')
j = ps.fields(f"j(2): {data_type}[2D]", layout='fzyx', field_type=ps.FieldType.STAGGERED_FLUX)
grad_c = sp.Matrix([ps.fd.diff(c, i) for i in range(c.spatial_dimensions)])
ek = ps.fd.FVM1stOrder(c, flux=-D * grad_c)
ast = ps.create_staggered_kernel(ek.discrete_flux(j))
code = ps.get_code_str(ast)
assert '_size_c_0 - 1 < _size_c_0 - 1' not in code
pystencils_tests/test_global_definitions.py → tests/test_global_definitions.py
View file @ 5600b6b6
......@@ -2,7 +2,7 @@ import sympy
import pystencils.astnodes
from pystencils.backends.cbackend import CBackend
from pystencils.data_types import TypedSymbol
from pystencils.typing import TypedSymbol
class BogusDeclaration(pystencils.astnodes.Node):
......@@ -95,7 +95,7 @@ def test_global_definitions_with_global_symbol():
z, x, y = pystencils.fields("z, y, x: [2d]")
normal_assignments = pystencils.AssignmentCollection([pystencils.Assignment(
z[0, 0], x[0, 0] * sympy.log(x[0, 0] * y[0, 0]))], [])
z[0, 0], x[0, 0] * x[0, 0] * y[0, 0])], [])
ast = pystencils.create_kernel(normal_assignments)
print(pystencils.show_code(ast))
......@@ -115,7 +115,7 @@ def test_global_definitions_without_global_symbol():
z, x, y = pystencils.fields("z, y, x: [2d]")
normal_assignments = pystencils.AssignmentCollection([pystencils.Assignment(
z[0, 0], x[0, 0] * sympy.log(x[0, 0] * y[0, 0]))], [])
z[0, 0], x[0, 0] * x[0, 0] * y[0, 0])], [])
ast = pystencils.create_kernel(normal_assignments)
print(pystencils.show_code(ast))
......
pystencils_tests/test_cudagpu.py → tests/test_gpu.py
View file @ 5600b6b6
import pytest
import numpy as np
import pycuda.autoinit
import pycuda.gpuarray as gpuarray
import sympy as sp
import math
from scipy.ndimage import convolve
from pystencils import Assignment, Field, fields
from pystencils.gpucuda import BlockIndexing, create_cuda_kernel, make_python_function
from pystencils.gpucuda.indexing import LineIndexing
from pystencils import Assignment, Field, fields, CreateKernelConfig, create_kernel, Target, get_code_str
from pystencils.gpu import BlockIndexing
from pystencils.simp import sympy_cse_on_assignment_list
from pystencils.slicing import add_ghost_layers, make_slice, remove_ghost_layers
from pystencils.slicing import add_ghost_layers, make_slice, remove_ghost_layers, normalize_slice
try:
import cupy as cp
device_numbers = range(cp.cuda.runtime.getDeviceCount())
except ImportError:
device_numbers = []
cp = None
def test_averaging_kernel():
pytest.importorskip('cupy')
size = (40, 55)
src_arr = np.random.rand(*size)
src_arr = add_ghost_layers(src_arr)
......@@ -22,13 +30,14 @@ def test_averaging_kernel():
update_rule = Assignment(dst_field[0, 0],
(src_field[0, 1] + src_field[0, -1] + src_field[1, 0] + src_field[-1, 0]) / 4)
ast = create_cuda_kernel(sympy_cse_on_assignment_list([update_rule]))
kernel = make_python_function(ast)
config = CreateKernelConfig(target=Target.GPU)
ast = create_kernel(sympy_cse_on_assignment_list([update_rule]), config=config)
kernel = ast.compile()
gpu_src_arr = gpuarray.to_gpu(src_arr)
gpu_dst_arr = gpuarray.to_gpu(dst_arr)
gpu_src_arr = cp.asarray(src_arr)
gpu_dst_arr = cp.asarray(dst_arr)
kernel(src=gpu_src_arr, dst=gpu_dst_arr)
gpu_dst_arr.get(dst_arr)
dst_arr = gpu_dst_arr.get()
stencil = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]]) / 4.0
reference = convolve(remove_ghost_layers(src_arr), stencil, mode='constant', cval=0.0)
......@@ -37,24 +46,26 @@ def test_averaging_kernel():
def test_variable_sized_fields():
pytest.importorskip('cupy')
src_field = Field.create_generic('src', spatial_dimensions=2)
dst_field = Field.create_generic('dst', spatial_dimensions=2)
update_rule = Assignment(dst_field[0, 0],
(src_field[0, 1] + src_field[0, -1] + src_field[1, 0] + src_field[-1, 0]) / 4)
ast = create_cuda_kernel(sympy_cse_on_assignment_list([update_rule]))
kernel = make_python_function(ast)
config = CreateKernelConfig(target=Target.GPU)
ast = create_kernel(sympy_cse_on_assignment_list([update_rule]), config=config)
kernel = ast.compile()
size = (3, 3)
src_arr = np.random.rand(*size)
src_arr = add_ghost_layers(src_arr)
dst_arr = np.zeros_like(src_arr)
gpu_src_arr = gpuarray.to_gpu(src_arr)
gpu_dst_arr = gpuarray.to_gpu(dst_arr)
gpu_src_arr = cp.asarray(src_arr)
gpu_dst_arr = cp.asarray(dst_arr)
kernel(src=gpu_src_arr, dst=gpu_dst_arr)
gpu_dst_arr.get(dst_arr)
dst_arr = gpu_dst_arr.get()
stencil = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]]) / 4.0
reference = convolve(remove_ghost_layers(src_arr), stencil, mode='constant', cval=0.0)
......@@ -63,6 +74,7 @@ def test_variable_sized_fields():
def test_multiple_index_dimensions():
pytest.importorskip('cupy')
"""Sums along the last axis of a numpy array"""
src_size = (7, 6, 4)
dst_size = src_size[:2]
......@@ -76,13 +88,14 @@ def test_multiple_index_dimensions():
update_rule = Assignment(dst_field[0, 0],
sum([src_field[offset[0], offset[1]](i) for i in range(src_size[-1])]))
ast = create_cuda_kernel([update_rule])
kernel = make_python_function(ast)
config = CreateKernelConfig(target=Target.GPU)
ast = create_kernel([update_rule], config=config)
kernel = ast.compile()
gpu_src_arr = gpuarray.to_gpu(src_arr)
gpu_dst_arr = gpuarray.to_gpu(dst_arr)
gpu_src_arr = cp.asarray(src_arr)
gpu_dst_arr = cp.asarray(dst_arr)
kernel(src=gpu_src_arr, dst=gpu_dst_arr)
gpu_dst_arr.get(dst_arr)
dst_arr = gpu_dst_arr.get()
reference = np.zeros_like(dst_arr)
gl = np.max(np.abs(np.array(offset, dtype=int)))
......@@ -94,6 +107,7 @@ def test_multiple_index_dimensions():
def test_ghost_layer():
pytest.importorskip('cupy')
size = (6, 5)
src_arr = np.ones(size)
dst_arr = np.zeros_like(src_arr)
......@@ -102,13 +116,15 @@ def test_ghost_layer():
update_rule = Assignment(dst_field[0, 0], src_field[0, 0])
ghost_layers = [(1, 2), (2, 1)]
ast = create_cuda_kernel([update_rule], ghost_layers=ghost_layers, indexing_creator=LineIndexing)
kernel = make_python_function(ast)
gpu_src_arr = gpuarray.to_gpu(src_arr)
gpu_dst_arr = gpuarray.to_gpu(dst_arr)
config = CreateKernelConfig(target=Target.GPU, ghost_layers=ghost_layers, gpu_indexing="line")
ast = create_kernel(sympy_cse_on_assignment_list([update_rule]), config=config)
kernel = ast.compile()
gpu_src_arr = cp.asarray(src_arr)
gpu_dst_arr = cp.asarray(dst_arr)
kernel(src=gpu_src_arr, dst=gpu_dst_arr)
gpu_dst_arr.get(dst_arr)
dst_arr = gpu_dst_arr.get()
reference = np.zeros_like(src_arr)
reference[ghost_layers[0][0]:-ghost_layers[0][1], ghost_layers[1][0]:-ghost_layers[1][1]] = 1
......@@ -116,25 +132,29 @@ def test_ghost_layer():
def test_setting_value():
pytest.importorskip('cupy')
arr_cpu = np.arange(25, dtype=np.float64).reshape(5, 5)
arr_gpu = gpuarray.to_gpu(arr_cpu)
arr_gpu = cp.asarray(arr_cpu)
iteration_slice = make_slice[:, :]
f = Field.create_generic("f", 2)
update_rule = [Assignment(f(0), sp.Symbol("value"))]
ast = create_cuda_kernel(update_rule, iteration_slice=iteration_slice, indexing_creator=LineIndexing)
kernel = make_python_function(ast)
config = CreateKernelConfig(target=Target.GPU, gpu_indexing="line", iteration_slice=iteration_slice)
ast = create_kernel(sympy_cse_on_assignment_list(update_rule), config=config)
kernel = ast.compile()
kernel(f=arr_gpu, value=np.float64(42.0))
np.testing.assert_equal(arr_gpu.get(), np.ones((5, 5)) * 42.0)
def test_periodicity():
from pystencils.gpucuda.periodicity import get_periodic_boundary_functor as periodic_gpu
pytest.importorskip('cupy')
from pystencils.gpu.periodicity import get_periodic_boundary_functor as periodic_gpu
from pystencils.slicing import get_periodic_boundary_functor as periodic_cpu
arr_cpu = np.arange(50, dtype=np.float64).reshape(5, 5, 2)
arr_gpu = gpuarray.to_gpu(arr_cpu)
arr_gpu = cp.asarray(arr_cpu)
periodicity_stencil = [(1, 0), (-1, 0), (1, 1)]
periodic_gpu_kernel = periodic_gpu(periodicity_stencil, (5, 5), 1, 2)
......@@ -143,23 +163,95 @@ def test_periodicity():
cpu_result = np.copy(arr_cpu)
periodic_cpu_kernel(cpu_result)
gpu_result = np.copy(arr_cpu)
periodic_gpu_kernel(pdfs=arr_gpu)
arr_gpu.get(gpu_result)
gpu_result = arr_gpu.get()
np.testing.assert_equal(cpu_result, gpu_result)
def test_block_indexing():
@pytest.mark.parametrize("device_number", device_numbers)
def test_block_indexing(device_number):
pytest.importorskip('cupy')
f = fields("f: [3D]")
bi = BlockIndexing(f, make_slice[:, :, :], block_size=(16, 8, 2), permute_block_size_dependent_on_layout=False)
s = normalize_slice(make_slice[:, :, :], f.spatial_shape)
bi = BlockIndexing(s, f.layout, block_size=(16, 8, 2),
permute_block_size_dependent_on_layout=False)
assert bi.call_parameters((3, 2, 32))['block'] == (3, 2, 32)
assert bi.call_parameters((32, 2, 32))['block'] == (16, 2, 8)
bi = BlockIndexing(f, make_slice[:, :, :], block_size=(32, 1, 1), permute_block_size_dependent_on_layout=False)
bi = BlockIndexing(s, f.layout, block_size=(32, 1, 1),
permute_block_size_dependent_on_layout=False)
assert bi.call_parameters((1, 16, 16))['block'] == (1, 16, 2)
bi = BlockIndexing(f, make_slice[:, :, :], block_size=(16, 8, 2), maximum_block_size="auto")
bi = BlockIndexing(s, f.layout, block_size=(16, 8, 2),
maximum_block_size="auto", device_number=device_number)
# This function should be used if number of needed registers is known. Can be determined with func.num_regs
blocks = bi.limit_block_size_by_register_restriction([1024, 1024, 1], 1000)
registers_per_thread = 1000
blocks = bi.limit_block_size_by_register_restriction([1024, 1024, 1], registers_per_thread)
if cp.cuda.runtime.is_hip:
max_registers_per_block = cp.cuda.runtime.deviceGetAttribute(71, device_number)
else:
device = cp.cuda.Device(device_number)
da = device.attributes
max_registers_per_block = da.get("MaxRegistersPerBlock")
assert np.prod(blocks) * registers_per_thread < max_registers_per_block
assert sum(blocks) < sum([1024, 1024, 1])
@pytest.mark.parametrize('gpu_indexing', ("block", "line"))
@pytest.mark.parametrize('layout', ("C", "F"))
@pytest.mark.parametrize('shape', ((5, 5, 5, 5), (3, 17, 387, 4), (23, 44, 21, 11)))
def test_four_dimensional_kernel(gpu_indexing, layout, shape):
pytest.importorskip('cupy')
n_elements = np.prod(shape)
arr_cpu = np.arange(n_elements, dtype=np.float64).reshape(shape, order=layout)
arr_gpu = cp.asarray(arr_cpu)
iteration_slice = make_slice[:, :, :, :]
f = Field.create_from_numpy_array("f", arr_cpu)
update_rule = [Assignment(f.center, sp.Symbol("value"))]
config = CreateKernelConfig(target=Target.GPU, gpu_indexing=gpu_indexing, iteration_slice=iteration_slice)
ast = create_kernel(update_rule, config=config)
kernel = ast.compile()
kernel(f=arr_gpu, value=np.float64(42.0))
np.testing.assert_equal(arr_gpu.get(), np.ones(shape) * 42.0)
@pytest.mark.parametrize('start', (1, 5))
@pytest.mark.parametrize('end', (-1, -2, -3, -4))
@pytest.mark.parametrize('step', (1, 2, 3, 4))
@pytest.mark.parametrize('shape', ([55, 60], [77, 101, 80], [44, 64, 66]))
def test_guards_with_iteration_slices(start, end, step, shape):
iter_slice = tuple([slice(start, end, step)] * len(shape))
kernel_config_gpu = CreateKernelConfig(target=Target.GPU, iteration_slice=iter_slice)
field_1 = fields(f"f(1) : double{list(shape)}")
assignment = Assignment(field_1.center, 1)
ast = create_kernel(assignment, config=kernel_config_gpu)
code_str = get_code_str(ast)
test_strings = list()
iteration_ranges = list()
for i, s in enumerate(iter_slice):
e = ((shape[i] + end) - s.start) / s.step
e = math.ceil(e) + s.start
test_strings.append(f"{s.start} < {e}")
a = s.start
counter = 0
while a < e:
a += 1
counter += 1
iteration_ranges.append(counter)
# check if the expected if statement is in the GPU code
for s in test_strings:
assert s in code_str
# check if these bounds lead to same lengths as the range function would produce
for i in range(len(iter_slice)):
assert iteration_ranges[i] == len(range(iter_slice[i].start, shape[i] + end, iter_slice[i].step))
tests/test_half_precision.py 0 → 100644
View file @ 5600b6b6
import pytest
import platform
import numpy as np
import pystencils as ps
@pytest.mark.parametrize('target', (ps.Target.CPU, ps.Target.GPU))
def test_half_precison(target):
if target == ps.Target.CPU:
if not platform.machine() in ['arm64', 'aarch64']:
pytest.xfail("skipping half precision test on non arm platform")
if 'clang' not in ps.cpu.cpujit.get_compiler_config()['command']:
pytest.xfail("skipping half precision because clang compiler is not used")
if target == ps.Target.GPU:
pytest.importorskip("cupy")
dh = ps.create_data_handling(domain_size=(10, 10), default_target=target)
f1 = dh.add_array("f1", values_per_cell=1, dtype=np.float16)
dh.fill("f1", 1.0, ghost_layers=True)
f2 = dh.add_array("f2", values_per_cell=1, dtype=np.float16)
dh.fill("f2", 2.0, ghost_layers=True)
f3 = dh.add_array("f3", values_per_cell=1, dtype=np.float16)
dh.fill("f3", 0.0, ghost_layers=True)
up = ps.Assignment(f3.center, f1.center + 2.1 * f2.center)
config = ps.CreateKernelConfig(target=dh.default_target, default_number_float=np.float32)
ast = ps.create_kernel(up, config=config)
kernel = ast.compile()
dh.run_kernel(kernel)
dh.all_to_cpu()
assert np.all(dh.cpu_arrays[f3.name] == 5.2)
assert dh.cpu_arrays[f3.name].dtype == np.float16