diff --git a/CMakeLists.txt b/CMakeLists.txt
index 15cdc03ad250bf7d9b197a8801a867f204ae4113..8dd89fe357e31548cbb51ee5f66bb84ccff03a99 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -10,7 +10,7 @@ include( PrepareSFG )
add_library( sfg_walberla INTERFACE )
target_sources( sfg_walberla
- INTERFACE include/GenericHbbBoundary.hpp
+ INTERFACE include/sfg/IndexVectors.hpp include/sfg/GenericHbbBoundary.hpp
)
target_include_directories( sfg_walberla INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/include )
diff --git a/cmake/PrepareSFG.cmake b/cmake/PrepareSFG.cmake
index c7551063b2ed0bb5938a6c2f69dd6b22da8216ae..ca0aa53570796e8f933751a04e55070ecd81cfd9 100644
--- a/cmake/PrepareSFG.cmake
+++ b/cmake/PrepareSFG.cmake
@@ -1,7 +1,7 @@
# Maybe set up private virtual environment
set(
- WALBERLA_CODEGEN_PRIVATE_VENV OFF
+ WALBERLA_CODEGEN_PRIVATE_VENV ON
CACHE BOOL
"Create a private virtual Python environment inside the build tree for code generation"
)
diff --git a/include/GenericHbbBoundary.hpp b/include/sfg/GenericHbbBoundary.hpp
similarity index 85%
rename from include/GenericHbbBoundary.hpp
rename to include/sfg/GenericHbbBoundary.hpp
index f68d5e12cbd304db15203ff9acc3ad878d5eeca3..f527606d3fd5f5bdedf276b56f40c33441df9b86 100644
--- a/include/GenericHbbBoundary.hpp
+++ b/include/sfg/GenericHbbBoundary.hpp
@@ -14,31 +14,13 @@
#include <memory>
#include <tuple>
+#include "IndexVectors.hpp"
+
namespace walberla
{
namespace sfg
{
-/**
- * @brief Index vectors for sparse-iteration boundary handling
- */
-template< typename LinkStruct >
-class BoundaryIndexVectors
-{
-public:
- using CpuIndexVector = std::vector< LinkStruct >;
-
- BoundaryIndexVectors() = default;
- bool operator==(const BoundaryIndexVectors& other) const { return other.links_ == links_; }
-
- CpuIndexVector& vector() { return links_; }
-
- LinkStruct* pointerCpu() { return links_.data(); }
-
-private:
- CpuIndexVector links_;
-};
-
struct HbbLink
{
int32_t x;
@@ -67,7 +49,7 @@ class GenericHbbBoundary
{
public:
using Stencil = Stencil_T;
- using IndexVectors = BoundaryIndexVectors< HbbLink >;
+ using IndexVectors = IndexListBuffer< HbbLink >;
GenericHbbBoundary(StructuredBlockForest& blocks)
{
diff --git a/include/sfg/IndexVectors.hpp b/include/sfg/IndexVectors.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9862e7662b63cceb337217392aa92cb55cf06739
--- /dev/null
+++ b/include/sfg/IndexVectors.hpp
@@ -0,0 +1,80 @@
+#pragma once
+
+#include <vector>
+#include "core/DataTypes.h"
+#include "core/cell/Cell.h"
+
+#include "domain_decomposition/BlockDataID.h"
+#include "blockforest/StructuredBlockForest.h"
+
+namespace walberla::sfg
+{
+ namespace internal
+ {
+ /**
+ * @brief Index vectors for sparse sweeps
+ */
+ template <typename IndexStruct>
+ class IndexListBuffer
+ {
+ public:
+ using CpuVector = std::vector<IndexStruct>;
+
+ IndexListBuffer() = default;
+ bool operator==(const IndexListBuffer &other) const { return other.vec_ == vec_; }
+
+ CpuVector &vector() { return vec_; }
+
+ IndexStruct *pointerCpu() { return vec_.data(); }
+
+ private:
+ CpuVector vec_;
+ };
+ }
+
+ struct CellIdx
+ {
+ int64_t x;
+ int64_t y;
+ int64_t z;
+
+ CellIdx(const Cell &c)
+ : x{int64_c(c.x())}, y{int64_c(c.y())}, z{int64_c(c.z())}
+ {
+ }
+
+ bool operator==(const CellIdx &other) const
+ {
+ return std::tie(x, y, z) == std::tie(other.x, other.y, other.z);
+ }
+ };
+
+ template <typename IndexStruct = CellIdx>
+ class SparseIndexList
+ {
+ public:
+ using Buffer_T = internal::IndexListBuffer<IndexStruct>;
+
+ explicit SparseIndexList(StructuredBlockForest &sbfs) : bufferId_{createBuffers(sbfs)} {}
+
+ std::vector<IndexStruct> &get(IBlock &block)
+ {
+ Buffer_T *buf = block.template getData<Buffer_T>(bufferId_);
+ return buf->vector();
+ }
+
+ const BlockDataID & bufferId() const {
+ return bufferId_;
+ }
+
+ private:
+ static BlockDataID createBuffers(StructuredBlockForest &sbfs)
+ {
+ auto createIdxVector = [](IBlock *const, StructuredBlockStorage *const)
+ { return new Buffer_T(); };
+ return sbfs.addStructuredBlockData<Buffer_T>(createIdxVector);
+ }
+
+ BlockDataID bufferId_;
+ };
+}
diff --git a/src/sfg_walberla/api.py b/src/sfg_walberla/api.py
index 2818b43f0ebe7bec7ef17859bc6857e7c40d8430..b6a1aeb7dcb177cd62adf1572580b74093511033 100644
--- a/src/sfg_walberla/api.py
+++ b/src/sfg_walberla/api.py
@@ -1,6 +1,6 @@
from __future__ import annotations
-from typing import Callable
+from typing import cast
from pystencils import Field
from pystencils.types import (
@@ -9,6 +9,7 @@ from pystencils.types import (
PsCustomType,
PsPointerType,
PsType,
+ PsStructType,
)
from pystencilssfg.lang import (
IFieldExtraction,
@@ -18,8 +19,9 @@ from pystencilssfg.lang import (
Ref,
ExprLike,
cpptype,
+ CppClass
)
-from pystencilssfg.lang.types import CppType
+from pystencilssfg.lang.types import CppTypeFactory, CppType
real_t = PsCustomType("walberla::real_t")
@@ -27,8 +29,8 @@ cell_idx_t = PsCustomType("walberla::cell_idx_t")
uint_t = PsCustomType("walberla::uint_t")
-class _CppClass(AugExpr):
- _type: Callable[..., CppType]
+class _PlainCppClass(AugExpr):
+ _type: CppTypeFactory
def __init__(self, const: bool = False, ref: bool = False):
dtype = self._type(const=const, ref=ref)
@@ -72,7 +74,7 @@ class Vector3(SrcVector):
return AugExpr(self._element_type).bind("{}[{}]", self, idx)
-class AABB(_CppClass):
+class AABB(_PlainCppClass):
_type = cpptype("walberla::AABB", "core/math/AABB.h")
def min(self) -> Vector3:
@@ -82,15 +84,15 @@ class AABB(_CppClass):
return Vector3(real_t).bind("{}.max()", self)
-class CellInterval(_CppClass):
+class CellInterval(_PlainCppClass):
_type = cpptype("walberla::CellInterval", "core/cell/CellInterval.h")
-class BlockDataID(_CppClass):
+class BlockDataID(_PlainCppClass):
_type = cpptype("walberla::BlockDataID", "domain_decomposition/BlockDataID.h")
-class IBlockPtr(_CppClass):
+class IBlockPtr(_PlainCppClass):
_type = cpptype("walberla::IBlock *", "domain_decomposition/IBlock.h")
def getData(self, dtype: str | PsType, id: BlockDataID) -> AugExpr:
@@ -103,34 +105,26 @@ class IBlockPtr(_CppClass):
return AugExpr.format("*{}", self)
+class SharedPtr(CppClass):
+ template = cpptype("std::shared_ptr< {} >", "<memory>")
+
+
class StructuredBlockForest(AugExpr):
+ typename = cpptype("walberla::StructuredBlockForest", "blockforest/StructuredBlockForest.h")
+
def __init__(self, ref: bool = True, const: bool = False):
- dtype = self.ref_type(const) if ref else self.cpp_type(const)
+ dtype = self.typename(const=const, ref=ref)
super().__init__(dtype)
@staticmethod
def shared_ptr(name: str = "blocks", const: bool = False):
- const_str = "const " if const else ""
- dtype = PsCustomType(
- f"std::shared_ptr< {const_str}walberla::StructuredBlockForest >", const=True
- )
- return AugExpr(dtype).var(name)
+ dtype = StructuredBlockForest.typename(const=const)
+ return SharedPtr(dtype).var(name)
@staticmethod
def shared_ptr_ref(name: str = "blocks", const: bool = False):
- const_str = "const " if const else ""
- dtype = PsCustomType(
- f"std::shared_ptr< {const_str}walberla::StructuredBlockForest >", const=True
- )
- return AugExpr(Ref(dtype)).var(name)
-
- @staticmethod
- def cpp_type(const: bool = False) -> PsType:
- return PsCustomType("walberla::StructuredBlockForest", const=const)
-
- @staticmethod
- def ref_type(const: bool = False) -> Ref:
- return Ref(StructuredBlockForest.cpp_type(const))
+ dtype = StructuredBlockForest.typename(const=const)
+ return SharedPtr(dtype, const=True, ref=True).var(name)
@staticmethod
def ref(name: str = "blocks", const: bool = False):
@@ -307,3 +301,56 @@ def glfield(field: Field, ci: str | AugExpr | None = None):
ci = AugExpr(PsCustomType("walberla::CellInterval", const=True)).var(ci)
return GhostLayerFieldExtraction(field_ptr, ci)
+
+
+class IndexListBufferPtr(SrcField):
+ _template = cpptype(
+ "walberla::sfg::internal::IndexListBuffer< {IndexStruct} >", "sfg/IndexVectors.hpp"
+ )
+
+ def __init__(self, idx_struct: PsStructType):
+ self._idx_struct = idx_struct
+ dtype = self._template(IndexStruct=idx_struct)
+ self._base_type = dtype
+ super().__init__(PsPointerType(dtype))
+
+ @property
+ def field_type(self) -> CppType:
+ return cast(CppType, self._base_type)
+
+ def pointerCpu(self) -> AugExpr:
+ return AugExpr().format("{}->pointerCpu()", self)
+
+ def vector(self) -> AugExpr:
+ return AugExpr().format("{}->vector()", self)
+
+ def get_extraction(self) -> IFieldExtraction:
+ class IdxVectorExtraction(IFieldExtraction):
+ def ptr(_) -> AugExpr:
+ return self.pointerCpu()
+
+ def size(_, coordinate: int) -> AugExpr | None:
+ if coordinate == 0:
+ return AugExpr.format("{}.size()", self.vector())
+ elif coordinate == 1:
+ return AugExpr.format("1")
+ else:
+ raise ValueError()
+
+ def stride(_, coordinate: int) -> AugExpr | None:
+ if coordinate <= 1:
+ return AugExpr.format("1")
+ else:
+ raise ValueError()
+
+ return IdxVectorExtraction()
+
+
+CellIdx = PsStructType(
+ (
+ ("x", create_type("int64_t")),
+ ("y", create_type("int64_t")),
+ ("z", create_type("int64_t"))
+ ),
+ "walberla::sfg::CellIdx"
+)
diff --git a/src/sfg_walberla/boundaries/hbb.py b/src/sfg_walberla/boundaries/hbb.py
index 2da0f1765def3493ae0c806c5a65d47f2c6f5ec2..e93ac79c3244e8683eda324596114a6d5ff2bd5d 100644
--- a/src/sfg_walberla/boundaries/hbb.py
+++ b/src/sfg_walberla/boundaries/hbb.py
@@ -2,26 +2,24 @@
from pystencils import Field, FieldType, TypedSymbol, Target, AssignmentCollection
-from pystencils.types import PsStructType, create_type, PsCustomType, PsPointerType
+from pystencils.types import PsStructType, create_type
from pystencilssfg import SfgComposer
from pystencilssfg.composer.class_composer import SfgClassComposer
from pystencilssfg.composer.custom import CustomGenerator
-from pystencilssfg.lang import SrcField, AugExpr
+from pystencilssfg.lang import AugExpr
from lbmpy.methods import AbstractLbMethod
from lbmpy.boundaries.boundaryconditions import LbBoundary
from lbmpy.advanced_streaming.indexing import Timestep
-from pystencilssfg.lang.expressions import IFieldExtraction
-
from ..sweep import (
FieldInfo,
SweepClassProperties,
BlockforestParameters,
combine_vectors,
)
-from ..api import GhostLayerFieldPtr, BlockDataID, IBlockPtr, StructuredBlockForest
+from ..api import GhostLayerFieldPtr, BlockDataID, IBlockPtr, StructuredBlockForest, IndexListBufferPtr
BoundaryIndexType = create_type("int32")
@@ -36,41 +34,6 @@ HbbLinkType = PsStructType(
)
-class BoundaryIndexVectorsPtr(SrcField):
- def __init__(self, link_struct: PsStructType):
- self._link_struct = link_struct
- typename = f"walberla::sfg::BoundaryIndexVectors< {link_struct.c_string()} >"
- self.vector_type = PsCustomType(typename)
- super().__init__(PsPointerType(self.vector_type))
-
- def pointerCpu(self) -> AugExpr:
- return AugExpr().format("{}->pointerCpu()", self)
-
- def vector(self) -> AugExpr:
- return AugExpr().format("{}->vector()", self)
-
- def get_extraction(self) -> IFieldExtraction:
- class IdxVectorExtraction(IFieldExtraction):
- def ptr(_) -> AugExpr:
- return self.pointerCpu()
-
- def size(_, coordinate: int) -> AugExpr | None:
- if coordinate == 0:
- return AugExpr.format("{}.size()", self.vector())
- elif coordinate == 1:
- return AugExpr.format("1")
- else:
- raise ValueError()
-
- def stride(_, coordinate: int) -> AugExpr | None:
- if coordinate <= 1:
- return AugExpr.format("1")
- else:
- raise ValueError()
-
- return IdxVectorExtraction()
-
-
class HbbBoundaryProperties(SweepClassProperties):
def _constructor(self, sfg: SfgComposer) -> SfgClassComposer.ConstructorBuilder:
blocks = StructuredBlockForest.shared_ptr_ref()
@@ -132,7 +95,7 @@ class SimpleHbbBoundary(CustomGenerator):
for s in sorted(parameters, key=lambda p: p.name):
props.add_property(s, setter=True, getter=True)
- idx_vector = BoundaryIndexVectorsPtr(HbbLinkType).var("indexVector")
+ idx_vector = IndexListBufferPtr(HbbLinkType).var("indexVector")
idx_vector_id = BlockDataID().bind("this->indexVectorsId_")
base_class = f"walberla::sfg::GenericHbbBoundary< walberla::stencil::{self._stencil.name} >"
@@ -148,16 +111,16 @@ class SimpleHbbBoundary(CustomGenerator):
),
# Get fields from block
*(
- sfg.init(fi.glfield_ptr)(
- block.getData(fi.glfield_ptr.field_type, fi.data_id)
+ sfg.init(fi.wlb_field_ptr)(
+ block.getData(fi.wlb_field_ptr.field_type, fi.data_id)
)
for fi in domain_fields
),
# Map GhostLayerFields to pystencils fields
- *(sfg.map_field(fi.field, fi.glfield_ptr) for fi in domain_fields),
+ *(sfg.map_field(fi.field, fi.wlb_field_ptr) for fi in domain_fields),
# Get index vector
sfg.init(idx_vector)(
- block.getData(idx_vector.vector_type, idx_vector_id)
+ block.getData(idx_vector.field_type, idx_vector_id)
),
# Map index vector to pystencils field
sfg.map_field(self._index_field, idx_vector),
diff --git a/src/sfg_walberla/sweep.py b/src/sfg_walberla/sweep.py
index f5960b8cbf9c7da58bc83728ba67c5fb2724b091..19fcc8d6b5f6f0374963b2c0d83d3b06276ed405 100644
--- a/src/sfg_walberla/sweep.py
+++ b/src/sfg_walberla/sweep.py
@@ -15,9 +15,11 @@ from pystencils import (
AssignmentCollection,
CreateKernelConfig,
Field,
+ FieldType,
Target,
+ create_type,
)
-from pystencils.types import PsType, constify, deconstify, PsCustomType
+from pystencils.types import PsType, constify, deconstify, PsCustomType, PsStructType
from pystencilssfg import SfgComposer
from pystencilssfg.lang import (
@@ -34,11 +36,13 @@ from .api import (
GenericWalberlaField,
GhostLayerFieldPtr,
GpuFieldPtr,
+ IndexListBufferPtr,
IBlockPtr,
BlockDataID,
Vector2,
Vector3,
CellInterval,
+ CellIdx,
)
@@ -134,7 +138,7 @@ class SweepClassProperties:
if var not in ctor.parameters:
if (
deconstify(strip_ptr_ref(var.dtype))
- == StructuredBlockForest.cpp_type()
+ == StructuredBlockForest.typename()
):
ctor.add_param(var, 0)
else:
@@ -257,9 +261,27 @@ class BlockforestParameters:
@dataclass
class FieldInfo:
field: Field
- glfield_ptr: GenericWalberlaField
+ wlb_field_ptr: GenericWalberlaField | IndexListBufferPtr
data_id: BlockDataID
+ def as_glfield(self) -> GenericWalberlaField:
+ if not isinstance(self.wlb_field_ptr, GenericWalberlaField):
+ raise AttributeError(
+ "This FieldInfo does not encapsulate a ghost layer field"
+ )
+ return self.wlb_field_ptr
+
+ def is_glfield(self) -> bool:
+ return isinstance(self.wlb_field_ptr, GenericWalberlaField)
+
+ def as_index_list(self) -> IndexListBufferPtr:
+ if not isinstance(self.wlb_field_ptr, IndexListBufferPtr):
+ raise AttributeError("This FieldInfo does not encapsulate an index list")
+ return self.wlb_field_ptr
+
+ def is_index_list(self) -> bool:
+ return isinstance(self.wlb_field_ptr, IndexListBufferPtr)
+
class ShadowFieldCache:
"""
@@ -277,12 +299,12 @@ class ShadowFieldCache:
def get_shadow(self, original_name: str) -> AugExpr:
original = self._originals[original_name]
return AugExpr.format(
- "this->{}({})", self._getter(original.field.name), original.glfield_ptr
+ "this->{}({})", self._getter(original.field.name), original.wlb_field_ptr
)
def perform_swap(self, original_name: str, shadow_field: FieldInfo) -> AugExpr:
original = self._originals[original_name]
- return original.glfield_ptr.swapDataPointers(shadow_field.glfield_ptr)
+ return original.as_glfield().swapDataPointers(shadow_field.wlb_field_ptr)
def render(self, sfg: SfgComposer):
if not self._originals:
@@ -295,22 +317,22 @@ class ShadowFieldCache:
for orig_name, orig in self._originals.items():
unique_ptr_type = PsCustomType(
- f"std::unique_ptr< {orig.glfield_ptr.field_type} >"
+ f"std::unique_ptr< {orig.as_glfield().field_type} >"
)
- cache_ptr_name = self._cache_ptr(orig.glfield_ptr)
+ cache_ptr_name = self._cache_ptr(orig.as_glfield())
cache_ptrs.append(sfg.var(cache_ptr_name, unique_ptr_type))
getters.append(
sfg.method(
self._getter(orig_name),
- returns=orig.glfield_ptr.get_dtype(),
+ returns=orig.wlb_field_ptr.get_dtype(),
inline=True,
)(
sfg.branch(f"{cache_ptr_name} == nullptr")(
AugExpr.format(
"{}.reset({});",
cache_ptr_name,
- orig.glfield_ptr.cloneUninitialized(),
+ orig.as_glfield().cloneUninitialized(),
)
),
f"return {cache_ptr_name}.get();",
@@ -410,15 +432,16 @@ class Sweep(CustomGenerator):
config = CreateKernelConfig(ghost_layers=0)
self._name = name
+
if isinstance(assignments, AssignmentCollection):
self._assignments = assignments
else:
self._assignments = AssignmentCollection(assignments)
self._gen_config = config
- if self._gen_config.target == Target.CUDA:
+ if self._gen_config.get_target() == Target.CUDA:
self._glfield_type = GpuFieldPtr
- elif self._gen_config.target.is_cpu():
+ elif self._gen_config.get_target().is_cpu():
self._glfield_type = GhostLayerFieldPtr
else:
raise ValueError(
@@ -428,6 +451,28 @@ class Sweep(CustomGenerator):
# Map from shadow field to shadowed field
self._shadow_fields: dict[Field, Field] = dict()
+ @property
+ def sparse(self) -> bool:
+ return self._gen_config.index_field is not None
+
+ @sparse.setter
+ def sparse(self, sparse_iteration: bool):
+ if sparse_iteration:
+ if self._gen_config.index_dtype != create_type("int64"):
+ raise ValueError(
+ "Sparse sweeps require `int64_t` as their index data type. Check your code generator config."
+ )
+
+ self._gen_config.index_field = Field.create_generic(
+ "indexList", 1, CellIdx, field_type=FieldType.INDEXED
+ )
+ self._gen_config.ghost_layers = None
+ else:
+ self._gen_config.index_field = None
+ self._gen_config.ghost_layers = 0
+
+ # CONFIGURATION
+
def swap_fields(self, field: Field, shadow_field: Field):
if field in self._shadow_fields:
raise ValueError(f"Field swap for {field} was already registered.")
@@ -440,6 +485,20 @@ class Sweep(CustomGenerator):
self._shadow_fields[shadow_field] = field
+ # CODE GENERATION
+
+ def _walberla_field(self, f: Field) -> GenericWalberlaField | IndexListBufferPtr:
+ match f.field_type:
+ case FieldType.GENERIC | FieldType.CUSTOM:
+ return self._glfield_type.create(f)
+ case FieldType.INDEXED:
+ assert isinstance(f.dtype, PsStructType)
+ return IndexListBufferPtr(f.dtype).var(f.name)
+ case _:
+ raise ValueError(
+ f"Unable to map field {f} of type {f.field_type} to a waLBerla field."
+ )
+
def generate(self, sfg: SfgComposer) -> None:
knamespace = sfg.kernel_namespace(f"{self._name}_kernels")
@@ -448,7 +507,7 @@ class Sweep(CustomGenerator):
all_fields: dict[str, FieldInfo] = {
f.name: FieldInfo(
- f, self._glfield_type.create(f), BlockDataID().var(f"{f.name}Id")
+ f, self._walberla_field(f), BlockDataID().var(f"{f.name}Id")
)
for f in khandle.fields
}
@@ -487,22 +546,29 @@ class Sweep(CustomGenerator):
*(sfg.init(fi.data_id)(props.get(fi.data_id)) for fi in block_fields),
# Extract regular fields from block
*(
- sfg.init(fi.glfield_ptr)(
- block.getData(fi.glfield_ptr.field_type, fi.data_id)
+ sfg.init(fi.wlb_field_ptr)(
+ block.getData(fi.wlb_field_ptr.field_type, fi.data_id)
)
for fi in block_fields
),
# Get shadow fields from cache
*(
- sfg.init(shadow_info.glfield_ptr)(
+ sfg.init(shadow_info.wlb_field_ptr)(
shadows_cache.get_shadow(orig_name)
)
for orig_name, shadow_info in swaps.items()
),
# Map GhostLayerFields to pystencils fields
*(
- sfg.map_field(fi.field, fi.glfield_ptr.with_cell_interval(ci))
+ sfg.map_field(fi.field, fi.as_glfield().with_cell_interval(ci))
+ for fi in all_fields.values()
+ if fi.is_glfield()
+ ),
+ # Map index lists to pystencils fields
+ *(
+ sfg.map_field(fi.field, fi.as_index_list())
for fi in all_fields.values()
+ if fi.is_index_list()
),
# Get parameters from class
*(sfg.init(param)(props.get(param)) for param in parameters),
@@ -521,13 +587,17 @@ class Sweep(CustomGenerator):
),
]
- sfg.klass(self._name)(
- *props.render(sfg),
- sfg.public(
- sfg.method("operator()")(*render_runmethod()),
+ methods = [sfg.method("operator()")(*render_runmethod())]
+
+ if not self.sparse:
+ methods.append(
sfg.method("runOnCellInterval")(
*render_runmethod(CellInterval(const=True, ref=True).var("ci"))
- ),
- ),
+ )
+ )
+
+ sfg.klass(self._name)(
+ *props.render(sfg),
+ sfg.public(*methods),
*shadows_cache.render(sfg),
)
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b75a729e6afd2fc0c610b87bdc10d746850297df
--- /dev/null
+++ b/tests/CMakeLists.txt
@@ -0,0 +1,19 @@
+cmake_minimum_required( VERSION 3.24 )
+project( sfg-walberla-testsuite )
+
+include(FetchContent)
+
+FetchContent_Declare(
+ walberla
+ GIT_REPOSITORY https://i10git.cs.fau.de/walberla/walberla.git
+)
+
+message( STATUS "Fetching waLBerla sources (this might take a while)..." )
+FetchContent_MakeAvailable(walberla)
+
+add_subdirectory(${CMAKE_SOURCE_DIR}/.. ${CMAKE_BINARY_DIR}/sfg-walberla)
+
+# Test Directories
+include(CTest)
+
+add_subdirectory( SparseIteration )
diff --git a/tests/CMakePresets.json b/tests/CMakePresets.json
new file mode 100644
index 0000000000000000000000000000000000000000..11b6c5797a4a848db2ef82c3e601e403848fb15a
--- /dev/null
+++ b/tests/CMakePresets.json
@@ -0,0 +1,29 @@
+{
+ "version": 6,
+ "cmakeMinimumRequired": {
+ "major": 3,
+ "minor": 23,
+ "patch": 0
+ },
+ "configurePresets": [
+ {
+ "name": "testsuite-dbg",
+ "binaryDir": "${sourceDir}/build/${presetName}",
+ "generator": "Ninja",
+ "cacheVariables": {
+ "WALBERLA_BUILD_TESTS": false
+ }
+ }
+ ],
+ "testPresets": [
+ {
+ "name": "SFG Testsuite",
+ "configurePreset": "testsuite-dbg",
+ "filter": {
+ "include": {
+ "name": "TestSparseSweeps"
+ }
+ }
+ }
+ ]
+}
diff --git a/tests/SparseIteration/CMakeLists.txt b/tests/SparseIteration/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..23666fabf4d41684ad3489000ac15993ab841fac
--- /dev/null
+++ b/tests/SparseIteration/CMakeLists.txt
@@ -0,0 +1,6 @@
+
+add_executable( TestSparseSweeps TestSparseSweeps.cpp )
+walberla_generate_sources( TestSparseSweeps SCRIPTS SparseSweep.py )
+target_link_libraries( TestSparseSweeps PRIVATE core blockforest field sfg_walberla )
+
+add_test( NAME TestSparseSweeps COMMAND TestSparseSweeps )
diff --git a/tests/SparseIteration/SparseSweep.py b/tests/SparseIteration/SparseSweep.py
new file mode 100644
index 0000000000000000000000000000000000000000..4489d3beee63edae16f324fb82c50a568134b22b
--- /dev/null
+++ b/tests/SparseIteration/SparseSweep.py
@@ -0,0 +1,20 @@
+from pystencilssfg import SourceFileGenerator
+from sfg_walberla import Sweep
+from sfg_walberla.symbolic import cell
+from pystencils import fields, Assignment
+
+with SourceFileGenerator() as sfg:
+ sfg.namespace("TestSparseSweeps::gen")
+ f = fields("f(3): [3D]", layout="fzyx")
+
+ asms = [
+ Assignment(f.center(0), cell.x()),
+ Assignment(f.center(1), cell.y()),
+ Assignment(f.center(2), cell.z()),
+ ]
+
+ sweep = Sweep("SparseSetCoordinates", asms)
+ sweep.sparse = True
+
+ sfg.generate(sweep)
+
diff --git a/tests/SparseIteration/TestSparseSweeps.cpp b/tests/SparseIteration/TestSparseSweeps.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..272930b78f835fe387a86806662b1275e706ac15
--- /dev/null
+++ b/tests/SparseIteration/TestSparseSweeps.cpp
@@ -0,0 +1,94 @@
+#include "core/all.h"
+#include "blockforest/all.h"
+#include "field/GhostLayerField.h"
+#include "field/AddToStorage.h"
+#include "sfg/IndexVectors.hpp"
+
+#include <array>
+
+#include "gen/TestSparseSweeps/SparseSweep.hpp"
+
+namespace TestSparseSweeps
+{
+ using namespace walberla;
+
+ using std::array;
+ using std::vector;
+
+ using PointField = field::GhostLayerField<real_t, 3>;
+
+ Vector3<real_t> pointOnSpiral(real_t z)
+ {
+ real_t two_pi_z{2.0 * math::pi * z};
+ return Vector3<real_t>{cos(two_pi_z), sin(two_pi_z), z};
+ }
+
+ void main(int argc, char **argv)
+ {
+
+ AABB domainAabb{0., 0., 0., 40., 40., 120.};
+ array<uint_t, 3> gridSize{4, 4, 12};
+ array<uint_t, 3> cpb{10, 10, 10};
+
+ auto blocks = blockforest::createUniformBlockGrid(
+ domainAabb,
+ gridSize[0], gridSize[1], gridSize[2],
+ cpb[0], cpb[1], cpb[2]);
+
+ BlockDataID fID = field::addToStorage<PointField>(blocks, "f", real_c(0.0), field::fzyx);
+ sfg::SparseIndexList indexList{*blocks};
+
+ auto isOnSpiral = [&](Cell globalCell) -> bool
+ {
+ Vector3<real_t> globalPoint{blocks->getCellCenter(globalCell)};
+ Vector3<real_t> spiralPoint{pointOnSpiral(globalPoint[2])};
+ return (globalPoint - spiralPoint).sqrLength() < 2.0;
+ };
+
+ CellInterval allCells{{0, 0, 0}, {cell_idx_c(cpb[0]), cell_idx_c(cpb[1]), cell_idx_c(cpb[2])}};
+
+ for (auto &b : *blocks)
+ {
+ vector<sfg::CellIdx> idxVec{indexList.get(b)};
+ for (Cell c : allCells)
+ {
+ Cell globalCell{c};
+ blocks->transformBlockLocalToGlobalCell(c, b);
+ if (isOnSpiral(globalCell))
+ {
+ idxVec.emplace_back(c);
+ }
+ }
+ }
+
+ gen::SparseSetCoordinates sweep{blocks, fID, indexList.bufferId()};
+
+ for (auto &b : *blocks)
+ {
+ sweep(&b);
+ }
+
+ for (auto &b : *blocks)
+ {
+ PointField & f = *b.getData<PointField>(fID);
+ for (Cell c : allCells)
+ {
+ Cell globalCell{c};
+ blocks->transformBlockLocalToGlobalCell(c, b);
+ if (isOnSpiral(c))
+ {
+ Vector3< real_t > cellCenter { blocks->getCellCenter(globalCell) };
+ WALBERLA_CHECK_FLOAT_EQUAL(f.get(c, 0), cellCenter[0]);
+ WALBERLA_CHECK_FLOAT_EQUAL(f.get(c, 1), cellCenter[1]);
+ WALBERLA_CHECK_FLOAT_EQUAL(f.get(c, 2), cellCenter[2]);
+ }
+ }
+ }
+ }
+}
+
+int main(int argc, char **argv)
+{
+ TestSparseSweeps::main(argc, argv);
+ return 0;
+}
\ No newline at end of file