Just rename from lib to modules

87 files changed, 6559 insertions, 0 deletions

diff --git a/modules/core/.nix/derivation.nix b/modules/core/.nix/derivation.nix
new file mode 100644
index 0000000..55a414f
--- /dev/null
+++ b/modules/core/.nix/derivation.nix
@@ -0,0 +1,36 @@
+{ lib
+, stdenv
+, scons
+, clang-tools
+, forstio
+, pname
+, version
+}:
+
+stdenv.mkDerivation {
+	inherit pname version;
+	src = ./..;
+
+  nativeBuildInputs = [
+    scons
+    clang-tools
+  ];
+
+  buildInputs = [
+    forstio.core
+    forstio.async
+    forstio.codec
+    forstio.codec-unit
+    forstio.codec-json
+  ];
+
+  doCheck = true;
+  checkPhase = ''
+    scons test
+    ./bin/tests
+	'';
+
+	preferLocalBuild = true;
+
+  outputs = [ "out" "dev" ];
+}
diff --git a/modules/core/SConstruct b/modules/core/SConstruct
new file mode 100644
index 0000000..8b3ab01
--- /dev/null
+++ b/modules/core/SConstruct
@@ -0,0 +1,75 @@
+#!/usr/bin/env python3
+
+import sys
+import os
+import os.path
+import glob
+import re
+
+
+if sys.version_info < (3,):
+    def isbasestring(s):
+        return isinstance(s,basestring)
+else:
+    def isbasestring(s):
+        return isinstance(s, (str,bytes))
+
+def add_kel_source_files(self, sources, filetype, lib_env=None, shared=False, target_post=""):
+
+    if isbasestring(filetype):
+        dir_path = self.Dir('.').abspath
+        filetype = sorted(glob.glob(dir_path+"/"+filetype))
+
+    for path in filetype:
+        target_name = re.sub( r'(.*?)(\.cpp|\.c\+\+)', r'\1' + target_post, path )
+        if shared:
+            target_name+='.os'
+            sources.append( self.SharedObject( target=target_name, source=path ) )
+        else:
+            target_name+='.o'
+            sources.append( self.StaticObject( target=target_name, source=path ) )
+    pass
+
+def isAbsolutePath(key, dirname, env):
+	assert os.path.isabs(dirname), "%r must have absolute path syntax" % (key,)
+
+env_vars = Variables(
+	args=ARGUMENTS
+)
+
+env_vars.Add('prefix',
+	help='Installation target location of build results and headers',
+	default='/usr/local/',
+	validator=isAbsolutePath
+)
+
+env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[],
+    CPPDEFINES=['SAW_UNIX'],
+    CXXFLAGS=[
+        '-std=c++20',
+        '-g',
+        '-Wall',
+        '-Wextra'
+    ],
+    LIBS=[
+        'forstio-core'
+    ]
+);
+env.__class__.add_source_files = add_kel_source_files
+env.Tool('compilation_db');
+env.cdb = env.CompilationDatabase('compile_commands.json');
+
+env.objects = [];
+env.sources = [];
+env.headers = [];
+env.targets = [];
+
+Export('env')
+SConscript('c++/SConscript')
+SConscript('tests/SConscript')
+
+env.Alias('cdb', env.cdb);
+env.Alias('all', [env.targets]);
+env.Default('all');
+
+env.Alias('install', '$prefix')
diff --git a/modules/core/c++/SConscript b/modules/core/c++/SConscript
new file mode 100644
index 0000000..91f5b3e
--- /dev/null
+++ b/modules/core/c++/SConscript
@@ -0,0 +1,35 @@
+#!/bin/false
+
+import os
+import os.path
+import glob
+
+
+Import('env')
+
+dir_path = Dir('.').abspath
+
+# Environment for base library
+core_env = env.Clone();
+
+core_env.sources = sorted(glob.glob(dir_path + "/*.cpp"));
+core_env.headers = sorted(glob.glob(dir_path + "/*.hpp"));
+
+core_env.particle_headers = sorted(glob.glob(dir_path + "/particle/*.hpp"));
+core_env.particle_geometry_headers = sorted(glob.glob(dir_path + "/particle/geometry/*.hpp"));
+
+core_env.math_headers = sorted(glob.glob(dir_path + "/math/*.hpp"));
+
+env.sources += core_env.sources;
+env.headers += core_env.headers;
+
+## Static lib
+objects = []
+core_env.add_source_files(objects, core_env.sources, shared=False);
+env.library_static = core_env.StaticLibrary('#build/kel-lbm', [objects]);
+
+env.Install('$prefix/lib/', env.library_static);
+env.Install('$prefix/include/kel/lbm/', core_env.headers);
+env.Install('$prefix/include/kel/lbm/particle/', core_env.particle_headers);
+env.Install('$prefix/include/kel/lbm/particle/geometry/', core_env.particle_geometry_headers);
+env.Install('$prefix/include/kel/lbm/math/', core_env.math_headers);
diff --git a/modules/core/c++/abstract/data.hpp b/modules/core/c++/abstract/data.hpp
new file mode 100644
index 0000000..ed23268
--- /dev/null
+++ b/modules/core/c++/abstract/data.hpp
@@ -0,0 +1,52 @@
+#pragma once
+
+#include "templates.hpp"
+
+namespace kel {
+namespace sch {
+struct Void {};
+
+struct UnsignedInteger {};
+struct SignedInteger {};
+struct FloatingPoint {};
+
+template<typename StorageT, typename InterfaceT>
+struct MixedPrecision {
+	using Meta = Void;
+	using StorageType = StorageT;
+	using InterfaceType = InterfaceT;
+};
+
+template<typename PrimType, uint64_t N>
+struct Primitive {
+	using Meta = Void;
+	using Type = PrimType;
+	static constexpr uint64_t Bytes = N;
+};
+
+template<typename T, uint64_t... Dims>
+struct FixedArray {
+	using Meta = Void;
+	using Inner = T;
+	static constexpr std::array<uint64_t,sizeof...(Dims)> Dimensions{Dims...};
+};
+
+template<typename T, uint64_t Dims>
+struct Array {
+	using Meta = FixedArray<UInt64,Dims>;
+	using Inner = T;
+	static constexpr std::array<uint64_t,sizeof...(Dims)> Dimensions{Dims};
+};
+
+template<typename... T>
+struct Tuple {
+};
+
+}
+
+template<typename Sch>
+struct schema {
+	using Type = Sch;
+};
+
+}
diff --git a/modules/core/c++/abstract/templates.hpp b/modules/core/c++/abstract/templates.hpp
new file mode 100644
index 0000000..d9d9faa
--- /dev/null
+++ b/modules/core/c++/abstract/templates.hpp
@@ -0,0 +1,8 @@
+#pragma once
+
+namespace kel {
+namespace lbm {
+template<typename... T>
+constexpr bool always_false = false;
+}
+}
diff --git a/modules/core/c++/args.hpp b/modules/core/c++/args.hpp
new file mode 100644
index 0000000..a0aa941
--- /dev/null
+++ b/modules/core/c++/args.hpp
@@ -0,0 +1,27 @@
+#pragma once
+
+#include <forstio/codec/data.hpp>
+#include <forstio/codec/args.hpp>
+
+namespace kel {
+namespace lbm {
+namespace sch {
+using namespace saw::schema;
+}
+
+template<typename ArgSchema>
+saw::error_or<saw::data<ArgSchema>> setup_lbm_env(int argc, char** argv){
+	saw::data<ArgSchema,saw::encode::Args> args_data{argc,argv};
+
+	saw::codec<ArgSchema,saw::encode::Args> args_codec;
+
+	saw::data<ArgSchema> args_decoded;
+	auto eov = args_codec.decode(args_data,args_decoded);
+	if(eov.is_error()){
+		return eov;
+	}
+
+	return args_decoded;
+}
+}
+}
diff --git a/modules/core/c++/boundary.hpp b/modules/core/c++/boundary.hpp
new file mode 100644
index 0000000..01ae7b5
--- /dev/null
+++ b/modules/core/c++/boundary.hpp
@@ -0,0 +1,249 @@
+#pragma once
+
+#include "component.hpp"
+#include "equilibrium.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct BounceBack {};
+
+template<uint64_t i>
+struct AntiBounceBack {};
+
+template<bool East>
+struct ZouHeHorizontal{};
+
+struct Equilibrium {};
+
+template<bool North>
+struct ZouHePressureY{};
+
+template<bool East>
+struct ZouHeVelocityX {};
+}
+
+/**
+ * Full-Way BounceBack.
+ */
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::BounceBack, Encode> {
+private:
+public:
+	component() = default;
+
+	using Component = cmpt::BounceBack;
+
+	/**
+	 * Thoughts regarding automating order setup
+	 */
+	using access = cmpt::access_tuple<
+		cmpt::access<"dfs", 1, true, true, true>
+	>;
+
+	/**
+	 * Thoughts regarding automating order setup
+	 */
+	static constexpr saw::string_literal name = "full_way_bounce_back";
+	static constexpr saw::string_literal after = "";
+	static constexpr saw::string_literal before = "streaming";
+
+	/**
+	 * Raw setup
+	 */
+	template<typename CellFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step)const{
+		bool is_even = ((time_step.get() % 2u) == 0u);
+
+		// This is a ref
+
+		// auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">();
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+
+		using dfi = df_info<T,Descriptor>;
+
+		// Technically use .copy()
+		auto& dfs_old = dfs_old_f.at(index);
+		auto df_cpy = dfs_old;
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			dfs_old.at({i}) = df_cpy.at({dfi::opposite_index.at(i)});
+		}
+	}
+};
+
+/**
+ * Anti Full-Way BounceBack.
+ */
+template<typename T, typename Descriptor, uint64_t D, typename Encode>
+class component<T, Descriptor, cmpt::AntiBounceBack<D>, Encode> {
+private:
+public:
+	component() = default;
+
+	using Component = cmpt::AntiBounceBack<D>;
+
+	/**
+	 * Thoughts regarding automating order setup
+	 */
+	using access = cmpt::access_tuple<
+		cmpt::access<"dfs", 1, true, true, true>
+	>;
+
+	/**
+	 * Thoughts regarding automating order setup
+	 */
+	static constexpr saw::string_literal name = "full_way_anti_bounce_back";
+	static constexpr saw::string_literal after = "";
+	static constexpr saw::string_literal before = "streaming";
+
+	/**
+	 * Raw setup
+	 */
+	template<typename CellFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step)const{
+		bool is_even = ((time_step.get() % 2u) == 0u);
+
+		// This is a ref
+
+		// auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">();
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+
+		using dfi = df_info<T,Descriptor>;
+
+		// Technically use .copy()
+		auto& dfs_old = dfs_old_f.at(index);
+		auto df_cpy = dfs_old;
+		
+		static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly");
+
+		dfs_old.at({0u}) = df_cpy.at({0u});
+		dfs_old.at({1u}) = df_cpy.at({1u});
+		dfs_old.at({2u}) = df_cpy.at({2u});
+		dfs_old.at({3u}) = df_cpy.at({4u});
+		dfs_old.at({4u}) = df_cpy.at({3u});
+		dfs_old.at({5u}) = df_cpy.at({7u});
+		dfs_old.at({6u}) = df_cpy.at({8u});
+		dfs_old.at({7u}) = df_cpy.at({5u});
+		dfs_old.at({8u}) = df_cpy.at({6u});
+	}
+};
+
+template<typename FP, typename Descriptor, typename Encode>
+class component<FP, Descriptor, cmpt::Equilibrium, Encode> final {
+private:
+	saw::data<sch::Scalar<FP>> density_;
+	saw::data<sch::Vector<FP,Descriptor::D>> velocity_;
+public:
+	component(
+		saw::data<sch::Scalar<FP>> density__,
+		saw::data<sch::Vector<FP,Descriptor::D>> velocity__
+	):
+		density_{density__},
+		velocity_{velocity__}
+	{}
+
+	template<typename CellFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step) const {
+		
+		bool is_even = ((time_step.get() % 2u) == 0u);
+		using dfi = df_info<FP,Descriptor>;
+		
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+
+		auto eq = equilibrium<FP,Descriptor>(density_,velocity_);
+
+		dfs_old_f.at(index) = eq;
+	}
+};
+
+/**
+ * This is massively hacky and expects a lot of conditions
+ * Either this or mirrored along the horizontal line works
+ *
+ * 0 - 2 - 2
+ * 0 - 3 - 1
+ * 0 - 3 - 1
+ * .........
+ * 0 - 3 - 1
+ * 0 - 2 - 2
+ *
+ */
+template<typename FP, typename Descriptor, bool Dir, typename Encode>
+class component<FP, Descriptor, cmpt::ZouHeHorizontal<Dir>, Encode> final {
+private:
+	saw::data<FP> rho_setting_;
+public:
+	component(const saw::data<FP>& rho_setting__):
+		rho_setting_{rho_setting__}
+	{}
+
+	template<typename CellFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+		using dfi = df_info<FP,Descriptor>;
+
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		auto& info_f = field.template get<"info">();
+		
+		// auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">();
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+
+		auto& dfs_old = dfs_old_f.at(index);
+
+		auto rho_vel_x = [&]() -> saw::data<FP> {
+			if constexpr (Dir){
+				auto S = dfs_old.at({0u})
+					+ dfs_old.at({3u}) + dfs_old.at({4u})
+					+ (dfs_old.at({1u}) + dfs_old.at({5u}) + dfs_old.at({7u})) * 2u;
+				return rho_setting_ - S;
+			}
+			else if constexpr (not Dir) {
+				auto S = dfs_old.at({0u})
+					+ dfs_old.at({3u}) + dfs_old.at({4u})
+					+ (dfs_old.at({2u}) + dfs_old.at({6u}) + dfs_old.at({8u})) * 2u;
+				return S - rho_setting_;
+			}
+			return {};
+		}();
+
+		static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly");
+
+		if constexpr (Dir) {
+			dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data<FP>{2.0 / 3.0} * rho_vel_x;
+			dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u}));
+			dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u}));
+		}else if constexpr (not Dir){
+			dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data<FP>{2.0 / 3.0} * rho_vel_x;
+			dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u}));
+			dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u}));
+		}
+	}
+};
+
+/*
+template<typename FP, typename Descriptor, bool East, typename Encode>
+class component<FP, Descriptor, cmpt::ZouHeVelocityX<East>, Encode> final {
+private:
+	saw::data<sch::Vector<FP,Descriptor::D>> velocity_;
+public:
+	component(
+		saw::data<sch::Vector<FP,Descriptor::D>> velocity__
+	):
+		velocity_{velocity__}
+	{}
+
+	template<typename CellFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step) const {
+		
+		bool is_even = ((time_step.get() % 2u) == 0u);
+		using dfi = df_info<FP,Descriptor>;
+		
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+
+	}
+};
+*/
+
+}
+}
diff --git a/modules/core/c++/chunk.hpp b/modules/core/c++/chunk.hpp
new file mode 100644
index 0000000..635af91
--- /dev/null
+++ b/modules/core/c++/chunk.hpp
@@ -0,0 +1,103 @@
+#pragma once
+
+#include "common.hpp"
+#include "flatten.hpp"
+
+namespace kel {
+namespace lbm {
+namespace sch {
+namespace impl {
+template<typename Sch, uint64_t Ghost, typename LeftG, typename RightG = saw::tmpl_value_group<uint64_t>>
+struct chunk_schema_type_helper;
+
+template<typename Sch, uint64_t Ghost, uint64_t Side0, uint64_t... Sides, uint64_t... AddedSides>
+struct chunk_schema_type_helper<Sch, Ghost, saw::tmpl_value_group<uint64_t,Side0,Sides...>, saw::tmpl_value_group<uint64_t,AddedSides...>> final {
+	using Schema = typename chunk_schema_type_helper<Sch,Ghost,saw::tmpl_value_group<uint64_t,Sides...>, saw::tmpl_value_group<uint64_t,AddedSides...,(Side0+2u*(Ghost))>>::Schema;
+};
+
+template<typename Sch, uint64_t Ghost, uint64_t... AddedSides>
+struct chunk_schema_type_helper<Sch, Ghost, saw::tmpl_value_group<uint64_t>, saw::tmpl_value_group<uint64_t,AddedSides...>> final {
+	using Schema = FixedArray<Sch,AddedSides...>;
+};
+}
+
+
+template<typename Sch, uint64_t Ghost, uint64_t... Sides>
+struct Chunk {
+	using InnerSchema = typename impl::chunk_schema_type_helper<Sch, Ghost, saw::tmpl_value_group<uint64_t,Sides...>>::Schema;
+	using StoredValueSchema = Sch;
+};
+
+// Not needed for now
+template<typename ChunkSchema, uint64_t Dim>
+using SuperChunk = Array<ChunkSchema,Dim>;
+}
+}
+}
+
+namespace saw {
+
+template<typename Sch, uint64_t Ghost, uint64_t... Sides, typename Encode>
+class data<kel::lbm::sch::Chunk<Sch,Ghost,Sides...>,Encode> final {
+public:
+	using Schema = kel::lbm::sch::Chunk<Sch,Ghost,Sides...>;
+private:
+	using InnerSchema = typename Schema::InnerSchema;
+	using ValueSchema = typename InnerSchema::ValueType;
+
+	data<InnerSchema, Encode> values_;
+public:
+	data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
+		return values_.at(index);
+	}
+
+	const data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
+		return values_.at(index);
+	}
+
+	static constexpr auto get_ghost_dims() {
+		return data<InnerSchema,Encode>::get_dims();
+	}
+
+	static constexpr auto to_ghost_index(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
+		std::decay_t<decltype(index)> ind;
+		for(uint64_t i = 0u; i < sizeof...(Sides); ++i){
+			ind.at({i}) = index.at({i}) + Ghost;
+		}
+		return ind;
+	}
+
+	data<ValueSchema, Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
+		std::decay_t<decltype(index)> ind = to_ghost_index(index);
+		return values_.at(ind);
+	}
+
+	const data<ValueSchema, Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
+		std::decay_t<decltype(index)> ind = to_ghost_index(index);
+		return values_.at(ind);
+	}
+
+	/*
+	const data<ValueSchema, Encode>& neighbour_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
+	}
+	*/
+
+	static constexpr auto get_dims(){
+		return data<schema::FixedArray<schema::UInt64, sizeof...(Sides)>,Encode>{{Sides...}};
+	}
+
+	auto flat_data(){
+		return values_.flat_data();
+	}
+
+	static constexpr auto flat_size() {
+		return data<InnerSchema,Encode>::flat_dims();
+	}
+};
+
+template<typename Sch, uint64_t Ghost, uint64_t... Sides>
+struct meta_schema<kel::lbm::sch::Chunk<Sch,Ghost,Sides...>> {
+	using MetaSchema = typename meta_schema<Sch>::MetaSchema;
+};
+
+}
diff --git a/modules/core/c++/collision.hpp b/modules/core/c++/collision.hpp
new file mode 100644
index 0000000..023f61f
--- /dev/null
+++ b/modules/core/c++/collision.hpp
@@ -0,0 +1,201 @@
+#pragma once
+
+#include "macroscopic.hpp"
+#include "component.hpp"
+#include "equilibrium.hpp"
+#include "chunk.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct BGK {};
+struct BGKGuo {};
+}
+
+/**
+ * Standard BGK collision operator for LBM
+ */
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::BGK, Encode> {
+private:
+	saw::data<T> relaxation_;
+	saw::data<T> frequency_;
+public:
+	component(typename saw::native_data_type<T>::type relaxation__):
+		relaxation_{{relaxation__}},
+		frequency_{saw::data<T>{1} / relaxation_}
+	{}
+	
+	component(const saw::data<T>& relaxation__):
+		relaxation_{relaxation__},
+		frequency_{saw::data<T>{1} / relaxation_}
+	{}
+
+	using Component = cmpt::BGK;
+
+	/**
+	 * Thoughts regarding automating order setup
+	 */
+	using access = cmpt::access_tuple<
+		cmpt::access<"dfs", 1, true, true, true>
+	>;
+
+	/**
+	 * Thoughts regarding automating order setup
+	 */
+	static constexpr saw::string_literal name = "collision";
+	static constexpr saw::string_literal after = "";
+	static constexpr saw::string_literal before = "streaming";
+
+	/**
+	 * Raw setup
+	 */
+	template<typename CellFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		// auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">();
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+
+		saw::data<sch::Scalar<T>> rho;
+		saw::data<sch::Vector<T,Descriptor::D>> vel;
+		compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel);
+		auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i}));
+		}
+	}
+
+	template<typename CellFieldSchema, typename MacroFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		// auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">();
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+
+		auto& rho_f = macros.template get<"density">();
+		auto& vel_f = macros.template get<"velocity">();
+
+		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
+		saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);
+
+		compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel);
+
+		auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i}));
+		}
+	}
+};
+
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::BGKGuo, Encode> {
+private:
+	typename saw::native_data_type<T>::type relaxation_;
+	saw::data<T> frequency_;
+	saw::data<sch::Vector<T,Descriptor::D>> global_force_;
+public:
+	component(typename saw::native_data_type<T>::type relaxation__):
+		relaxation_{relaxation__},
+		frequency_{typename saw::native_data_type<T>::type(1) / relaxation_},
+		global_force_{}
+	{}
+
+	component(
+		typename saw::native_data_type<T>::type relaxation__,
+		saw::data<sch::Vector<T,Descriptor::D>>& global_force__
+	):
+		relaxation_{relaxation__},
+		frequency_{typename saw::native_data_type<T>::type(1) / relaxation_},
+		global_force_{global_force__}
+	{}
+
+	using Component = cmpt::BGKGuo;
+
+	using access = cmpt::access_tuple<
+		cmpt::access<"dfs", 1, true, true, true>,
+		cmpt::access<"force", 0, true, false>
+	>;
+
+	/**
+	 * Thoughts regarding automating order setup
+	 */
+	static constexpr saw::string_literal name = "collision";
+	static constexpr saw::string_literal after = "";
+	static constexpr saw::string_literal before = "streaming";
+
+	template<typename CellFieldSchema, typename MacroFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+	// void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, saw::data<sch::UInt64> time_step){
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+		auto& dfs = dfs_old_f.at(index);
+
+		auto& rho_f = macros.template get<"density">();
+		auto& vel_f = macros.template get<"velocity">();
+	
+		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
+		
+		auto& force_f = macros.template get<"force">();
+		auto& force = force_f.at(index);
+
+		auto total_force = force + global_force_;
+
+		saw::data<sch::Scalar<T>> half;
+		half.at({}).set(0.5);
+		auto& vel = vel_f.at(index);
+		vel = vel + total_force * ( half / rho );
+
+		compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel);
+		auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+		using dfi = df_info<T,Descriptor>;
+
+		
+		saw::data<sch::Scalar<T>> dfi_inv_cs2;
+		dfi_inv_cs2.at({}).set(dfi::inv_cs2);
+			
+
+		// auto vel = vel_f.at(index);
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			// saw::data<T> ci_min_u{0};
+			saw::data<sch::Vector<T,Descriptor::D>> ci;
+			for(uint64_t d = 0u; d < Descriptor::D; ++d){
+				ci.at({{d}}).set(static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d]));
+			}
+			auto ci_dot_u = saw::math::dot(ci,vel);
+
+			// saw::data<sch::Vector<T,Descriptor::D>> F_i;
+			// F_i = f * ((c_i - u) * ics2  + <c_i,u> * c_i * ics2 * ics2) * w_i;
+			saw::data<sch::Scalar<T>> w;
+			w.at({}).set(dfi::weights[i]);
+
+			/*
+			saw::data<sch::Scalar<T>> F_i_sum;
+			for(uint64_t d = 0u; d < Descriptor::D; ++d){
+				saw::data<sch::Scalar<T>> F_i_d;
+				F_i_d.at({}) = F_i.at({{d}});
+				F_i_sum = F_i_sum + F_i_d;
+			}
+			*/
+			auto term1 = (ci-vel) * dfi_inv_cs2;
+			auto term2 = ci * (ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2);
+
+			auto force_projection = saw::math::dot(term1 + term2, total_force);
+
+			auto F_i = w * force_projection;
+
+			dfs.at({i}) = dfs.at({i})
+										+ frequency_
+											* (eq.at(i) - dfs.at({i}) )
+										+ F_i.at({})
+											* (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_);
+		}
+	}
+};
+}
+}
diff --git a/modules/core/c++/common.hpp b/modules/core/c++/common.hpp
new file mode 100644
index 0000000..280d953
--- /dev/null
+++ b/modules/core/c++/common.hpp
@@ -0,0 +1,12 @@
+#pragma once
+
+#include <forstio/codec/data.hpp>
+#include <forstio/codec/math.hpp>
+
+namespace kel {
+namespace lbm {
+namespace sch {
+using namespace saw::schema;
+}
+}
+}
diff --git a/modules/core/c++/component.hpp b/modules/core/c++/component.hpp
new file mode 100644
index 0000000..1e5dbbf
--- /dev/null
+++ b/modules/core/c++/component.hpp
@@ -0,0 +1,38 @@
+#pragma once
+
+#include "descriptor.hpp"
+
+namespace kel {
+namespace lbm {
+
+namespace cmpt {
+
+/**
+ * Maybe for the future
+ */
+template<saw::string_literal Name, uint64_t Dist, bool Read, bool Write, bool SkipSync = false>
+class access {
+public:
+	static constexpr saw::string_literal name = Name;
+	static constexpr uint64_t distance = Dist;
+	static constexpr bool read = Read;
+	static constexpr bool write = Write;
+	static constexpr bool skip_sync = SkipSync;
+};
+
+/**
+ * Maybe for the future
+ */
+template<typename... Acc>
+class access_tuple {
+public:
+};
+}
+
+/**
+ * Component class which forms the basis of each processing step
+ */
+template<typename T, typename Descriptor, typename Cmpt, typename Encode = saw::encode::Native>
+class component {};
+}
+}
diff --git a/modules/core/c++/config.hpp b/modules/core/c++/config.hpp
new file mode 100644
index 0000000..64f7a0f
--- /dev/null
+++ b/modules/core/c++/config.hpp
@@ -0,0 +1,53 @@
+#pragma once
+
+#include <forstio/codec/data.hpp>
+#include <forstio/codec/json/json.hpp>
+
+#include <fstream>
+#include <sstream>
+#include <string_view>
+#include <string>
+
+namespace kel {
+namespace lbm {
+namespace sch {
+using namespace saw::schema;
+template<typename T, typename Desc>
+using LbmConfig = Struct<
+    Member<T, "delta_x">,
+    Member<T, "kinematic_viscosity">,
+    Member<T, "delta_t">
+>;
+}
+
+template<typename T, typename Desc>
+saw::error_or<saw::data<sch::LbmConfig<T,Desc>>> load_lbm_config(std::string_view file_name){
+    std::ifstream file{std::string{file_name}};
+
+    if(!file.is_open()){
+        return saw::make_error<saw::err::not_found>("Couldn't open file");
+    }
+
+    saw::data<sch::LbmConfig<T,Desc>, saw::encode::Json> lbm_json_conf{saw::heap<saw::array_buffer>(1u)};
+
+    uint8_t ele{};
+    while(file.readsome(reinterpret_cast<char*>(&ele), 1u) > 0u){
+        auto err = lbm_json_conf.get_buffer().push(ele,1u);
+        if(err.failed()){
+            return err;
+        }
+    }
+
+    saw::data<sch::LbmConfig<T,Desc>> lbm_conf;
+    saw::codec<sch::LbmConfig<T,Desc>, saw::encode::Json> json_codec;
+    {
+        auto eov = json_codec.decode(lbm_json_conf, lbm_conf);
+        if(eov.is_error()){
+            return std::move(eov.get_error());
+        }
+    }
+
+    return lbm_conf;
+}
+}
+}
diff --git a/modules/core/c++/converter.hpp b/modules/core/c++/converter.hpp
new file mode 100644
index 0000000..4370a2c
--- /dev/null
+++ b/modules/core/c++/converter.hpp
@@ -0,0 +1,87 @@
+#pragma once
+
+#include "lbm_unit.hpp"
+#include "descriptor.hpp"
+
+namespace kel {
+namespace lbm {
+
+namespace sch {
+using namespace saw::schema;
+}
+
+/**
+ * Helps converting from SI types to LBM types
+ */
+template<typename T>
+class converter {
+private:
+	saw::data<typename saw::unit_division<sch::SiMeter<T>, sch::LbmMeter<T> >::Schema > meter_conv_; 
+	saw::data<typename saw::unit_division<sch::SiSecond<T>, sch::LbmSecond<T> >::Schema > second_conv_; 
+public:
+	converter() = delete;
+	converter(
+		saw::data<typename saw::unit_division<sch::SiMeter<T>, sch::LbmMeter<T> >::Schema > meter_conv__,
+		saw::data<typename saw::unit_division<sch::SiSecond<T>, sch::LbmSecond<T> >::Schema > second_conv__
+	):
+		meter_conv_{meter_conv__},
+		second_conv_{second_conv__}
+	{}
+
+	/**
+	 * Get the conversion parameter with the conversion type
+	 */
+	auto conversion_x() const {
+		return meter_conv_;
+	}
+
+	/**
+	 * Get the conversion parameter with the conversion type
+	 */
+	auto conversion_t() const {
+		return second_conv_;
+	}
+
+	auto delta_x() const {
+		return meter_conv_*saw::data<sch::LbmMeter<T>>{1.0};
+	}
+
+	auto delta_t() const {
+		return second_conv_*saw::data<sch::LbmSecond<T>>{1.0};
+	}
+
+	auto delta_v() const {
+		return (meter_conv_ / second_conv_) * saw::data<sch::LbmVelocity<T>>{1.0};
+	}
+	
+	auto delta_a() const {
+		return (meter_conv_ / (second_conv_*second_conv_)) * saw::data<sch::LbmAcceleration<T>>{1.0};
+	}
+
+	saw::data<sch::LbmMeter<T>> meter_si_to_lbm(const saw::data<sch::SiMeter<T>>& m_si) const {
+		return m_si / meter_conv_;
+	}
+	
+	saw::data<sch::LbmSecond<T>> second_si_to_lbm(const saw::data<sch::SiSecond<T>>& s_si) const {
+		return s_si / second_conv_;
+	}
+	
+	saw::data<sch::LbmVelocity<T>> velocity_si_to_lbm(const saw::data<sch::SiVelocity<T>>& vel_si) const {
+		return vel_si * second_conv_ / meter_conv_;
+	}
+	
+	saw::data<sch::LbmAcceleration<T>> acceleration_si_to_lbm(const saw::data<sch::SiAcceleration<T>>& acc_si) const {
+		return acc_si * (second_conv_ * second_conv_) / meter_conv_;
+	}
+
+	saw::data<sch::LbmKinematicViscosity<T>> kinematic_viscosity_si_to_lbm (const saw::data<sch::SiKinematicViscosity<T>>& kin_si) const {
+		return (kin_si / (meter_conv_ * meter_conv_)) * second_conv_;
+	}
+
+	template<typename Desc>
+	saw::data<sch::Pure<T>> kinematic_viscosity_si_to_tau(const saw::data<sch::SiKinematicViscosity<T>>& kin_si) const {
+		return saw::data<sch::Pure<T>>{saw::data<typename saw::unit_division<sch::Pure<T>, sch::LbmKinematicViscosity<T>>::Schema >{df_info<T,Desc>::inv_cs2} * kinematic_viscosity_si_to_lbm(kin_si) + saw::data<sch::Pure<T>>{0.5}};
+	}
+};
+}
+}
diff --git a/modules/core/c++/descriptor.hpp b/modules/core/c++/descriptor.hpp
new file mode 100644
index 0000000..9f7399a
--- /dev/null
+++ b/modules/core/c++/descriptor.hpp
@@ -0,0 +1,438 @@
+#pragma once
+
+#include <forstio/codec/data.hpp>
+#include <forstio/codec/data_math.hpp>
+#include <forstio/codec/schema_factory.hpp>
+
+namespace kel {
+namespace lbm {
+namespace sch {
+using namespace saw::schema;
+
+template<uint64_t DV, uint64_t QV>
+struct Descriptor {
+	static constexpr uint64_t D = DV;
+	static constexpr uint64_t Q = QV;
+};
+
+using D2Q9 = Descriptor<2u,9u>;
+//using D2Q5 = Descriptor<2u,5u>;
+using D3Q27 = Descriptor<3u,27u>;
+
+template<typename Sch, typename Desc, uint64_t SC_V, uint64_t DC_V, uint64_t QC_V>
+struct Cell {
+	using Descriptor = Desc;
+	static constexpr uint64_t SC = SC_V;
+	static constexpr uint64_t DC = DC_V;
+	static constexpr uint64_t QC = QC_V;
+	static constexpr uint64_t Size = SC + Desc::D * DC + Desc::Q * QC;
+};
+
+template<typename Desc, typename Cell>
+struct CellField{
+	using Descriptor = Desc;
+};
+
+template<typename Desc, typename... CellFieldTypes, saw::string_literal... CellFieldNames>
+struct CellField<
+	Desc,
+	Struct<
+		Member<CellFieldTypes, CellFieldNames>...
+	>
+> {
+	using Descriptor = Desc;
+};
+
+template<typename Desc, typename... CellFieldMembers>
+struct CellFieldStruct {
+	using Descriptor = Desc;
+	// using MetaSchema = FixedArray<UInt64, Desc::D>;
+};
+
+}
+
+template<typename T, typename Desc>
+class df_info{};
+
+/*
+namespace impl {
+template<typename Desc>
+struct df_ct_helper {
+	template<uint64_t i>
+	static constexpr uint64_t apply_i(const std::array<std::array<int32_t,Desc::D>,Desc::Q>& dirs, const std::array<T,Desc::D+1u>& inp){
+		if constexpr ( i < Desc::Q ){
+			for(uint64_t j = 0u; j < Desc::D; ++j){
+			}
+		}
+		return 0u;
+	}
+};
+}
+*/
+
+template<typename T> 
+class df_info<T,sch::Descriptor<1,3>> {
+public:
+	using Descriptor = sch::Descriptor<1,3>;
+
+	static constexpr uint64_t D = 1u;
+	static constexpr uint64_t Q = 3u;
+
+	static constexpr std::array<std::array<int32_t,D>,Q> directions = {{
+		{ 0},
+		{-1},
+		{ 1}
+	}};
+
+	static constexpr std::array<typename saw::native_data_type<T>::type, Q> weights = {
+		2./3.,
+		1./6.,
+		1./6.
+	};
+
+	static constexpr std::array<uint64_t,Q> opposite_index = {
+		0,2,1
+	};
+	
+	static constexpr typename saw::native_data_type<T>::type inv_cs2 = 3.0;
+	static constexpr typename saw::native_data_type<T>::type cs2 = 1./3.;
+};
+
+/**
+ * D2Q5 Descriptor
+ */
+template<typename T>
+class df_info<T,sch::Descriptor<2, 5>> {
+public:
+	using Descriptor = sch::Descriptor<2,5>;
+
+	static constexpr uint64_t D = 2u;
+	static constexpr uint64_t Q = 5u;
+
+	static constexpr std::array<std::array<int32_t, D>, Q> directions = {{
+		{ 0, 0},
+		{-1, 0},
+		{ 1, 0},
+		{ 0,-1},
+		{ 0, 1},
+	}};
+
+	static constexpr std::array<typename saw::native_data_type<T>::type,Q> weights = {
+		1./3.,
+		1./6.,
+		1./6.,
+		1./6.,
+		1./6.
+	};
+
+	static constexpr std::array<uint64_t,Q> opposite_index = {
+		0,
+		2,
+		1,
+		4,
+		3
+	};
+
+	static constexpr typename saw::native_data_type<T>::type inv_cs2 = 3.0;
+	static constexpr typename saw::native_data_type<T>::type cs2 = 1./3.;
+};
+
+template<typename T>
+class df_info<T,sch::Descriptor<2, 9>> {
+public:
+	using Descriptor = sch::Descriptor<2,9>;
+	
+	static constexpr uint64_t D = 2u;
+	static constexpr uint64_t Q = 9u;
+
+	static constexpr std::array<std::array<int32_t, D>, Q> directions = {{
+		{ 0, 0}, // 0
+		{-1, 0}, // 1
+		{ 1, 0}, // 2
+		{ 0,-1}, // 3
+		{ 0, 1}, // 4
+		{-1,-1}, // 5
+		{ 1, 1}, // 6
+		{-1, 1}, // 7
+		{ 1,-1}  // 8
+	}};
+
+	static constexpr std::array<typename saw::native_data_type<T>::type,Q> weights = {
+		4./9.,
+		1./9.,
+		1./9.,
+		1./9.,
+		1./9.,
+		1./36.,
+		1./36.,
+		1./36.,
+		1./36.
+	};
+
+	static constexpr std::array<uint64_t,Q> opposite_index = {
+		0,
+		2,
+		1,
+		4,
+		3,
+		6,
+		5,
+		8,
+		7
+	};
+	
+	static constexpr typename saw::native_data_type<T>::type inv_cs2 = 3.0;
+	static constexpr typename saw::native_data_type<T>::type cs2 = 1./3.;
+};
+
+
+template<typename T>
+class df_info<T,sch::Descriptor<3, 27>> {
+public:
+	using Descriptor = sch::Descriptor<3,27>;
+
+	static constexpr uint64_t D = 3u;
+	static constexpr uint64_t Q = 27u;
+
+	static constexpr std::array<std::array<int32_t, D>, Q> directions = {{
+		{ 0, 0, 0},  // 0
+		// Into 1D
+		{-1, 0, 0},  // 1
+		{ 1, 0, 0},  // 2
+		// Expand into 2D
+		{ 0, -1, 0}, // 3
+		{-1, -1, 0}, // 4
+		{ 1, -1, 0}, // 5
+		{ 0, 1, 0},  // 6
+		{-1, 1, 0},  // 7
+		{ 1, 1, 0},  // 8
+		// Expand into 3D
+		{ 0, 0, -1}, // 9
+		{-1, 0, -1}, // 10
+		{ 1, 0, -1}, // 11
+		{ 0, -1, -1},// 12
+		{-1, -1, -1},// 13
+		{ 1, -1, -1},// 14
+		{ 0, 1, -1}, // 15
+		{-1, 1, -1}, // 16
+		{ 1, 1, -1}, // 17
+		{ 0, 0, 1},  // 18
+		{-1, 0, 1},  // 19
+		{ 1, 0, 1},  // 20
+		{ 0, -1, 1}, // 21
+		{-1, -1, 1}, // 22
+		{ 1, -1, 1}, // 23
+		{ 0, 1, 1},  // 24
+		{-1, 1, 1},  // 25
+		{ 1, 1, 1}   // 26
+	}};
+
+	static constexpr std::array<typename saw::native_data_type<T>::type,Q> weights = {
+		8./27.,
+		// 1D
+		1./36.,
+		1./36.,
+		// 2D
+		1./36.,
+		1./54.,
+		1./54.,
+		1./36.,
+		1./54.,
+		1./54.,
+		// 3D
+		1./36.,
+		1./54.,
+		1./54.,
+		1./54.,
+		1./216.,
+		1./216.,
+		1./54.,
+		1./216.,
+		1./216.,
+		1./36.,
+		1./54.,
+		1./54.,
+		1./54.,
+		1./216.,
+		1./216.,
+		1./54.,
+		1./216.,
+		1./216.
+	};
+
+	static constexpr std::array<uint64_t,Q> opposite_index = {
+		0,2,1,
+		6,8,7,3,5,4,
+		18,20,19,24,26,25,21,23,22,9,11,10,15,17,16,12,14,13
+	};
+
+	static constexpr typename saw::native_data_type<T>::type inv_cs2 = 3.0;
+	static constexpr typename saw::native_data_type<T>::type cs2 = 1./3.;
+};
+
+
+template<typename Schema>
+class cell_schema_builder {
+private:
+	saw::schema_factory<Schema> factory_struct_;
+public:
+	cell_schema_builder() = default;
+
+	cell_schema_builder(saw::schema_factory<Schema> inp):
+		factory_struct_{inp}
+	{}
+
+	/*
+	template<typename TA, saw::string_literal KA>
+	constexpr auto require() const noexcept {
+		return {factory_struct_.add_maybe()};
+	}
+	*/
+};
+ 
+}
+}
+
+namespace saw {
+template<typename T, typename Desc, uint64_t S, uint64_t D, uint64_t Q>
+struct meta_schema<kel::lbm::sch::Cell<T,Desc,S,D,Q>> {
+	using MetaSchema = schema::Void;
+	using Schema = kel::lbm::sch::Cell<T,Desc,S,D,Q>;
+};
+
+template<typename Desc, typename CellT>
+struct meta_schema<kel::lbm::sch::CellField<Desc, CellT>> {
+	using MetaSchema = schema::FixedArray<schema::UInt64,Desc::D>;
+	using Schema = kel::lbm::sch::CellField<Desc, CellT>;
+};
+
+template<typename Desc, typename... CellFieldsT>
+struct meta_schema<kel::lbm::sch::CellFieldStruct<Desc,schema::Struct<CellFieldsT...>>> {
+	using MetaSchema = schema::FixedArray<schema::UInt64,Desc::D>;
+	using Schema = kel::lbm::sch::CellFieldStruct<Desc,schema::Struct<CellFieldsT...>>;
+};
+
+template<typename Sch, typename Desc, uint64_t S, uint64_t D, uint64_t Q, typename Encode>
+class data<kel::lbm::sch::Cell<Sch, Desc, S, D, Q>, Encode> final {
+public:
+	using Schema = kel::lbm::sch::Cell<Sch,Desc,S,D,Q>;
+	using MetaSchema = typename meta_schema<Schema>::MetaSchema;
+private:
+	data<schema::FixedArray<Sch, Schema::Size>, Encode> inner_;
+public:
+	data() = default;
+
+	data<Sch, Encode>& operator()(const data<schema::UInt64>& index){
+		return inner_.at(index);
+	}
+
+	const data<Sch, Encode>& operator()(const data<schema::UInt64>& index)const{
+		return inner_.at(index);
+	}
+
+	const data<kel::lbm::sch::Cell<Sch, Desc, S, D, Q>, Encode> copy() const {
+		return *this;
+	}
+};
+
+/*
+ * Create a cellfield
+ */
+template<typename Desc, typename CellT, typename Encode>
+class data<kel::lbm::sch::CellField<Desc, CellT>, Encode> final {
+public:
+	using Schema = kel::lbm::sch::CellField<Desc,CellT>;
+	using MetaSchema = typename meta_schema<Schema>::MetaSchema;
+private:
+	data<schema::Array<CellT,Desc::D>, Encode> inner_;
+public:
+	data() = default;
+	data(const data<MetaSchema,Encode>& inner_meta__):
+		inner_{inner_meta__}
+	{}
+
+	const data<MetaSchema, Encode> meta() const {
+		return inner_.dims();
+	}
+
+	template<uint64_t i>
+	data<schema::UInt64,Encode> get_dim_size() const {
+		static_assert(i < Desc::D, "Not enough dimensions");
+		return inner_.template get_dim_size<i>();
+	}
+
+	const data<CellT>& operator()(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index)const{
+		return inner_.at(index);
+	}
+
+	data<CellT>& operator()(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index){
+		return inner_.at(index);
+	}
+
+	const data<CellT>& at(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index)const{
+		return inner_.at(index);
+	}
+
+	data<CellT>& at(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index){
+		return inner_.at(index);
+	}
+
+	data<schema::UInt64,Encode> internal_size() const {
+		return inner_.internal_size();
+	}
+
+	data<CellT,Encode>* internal_data() {
+		return inner_.internal_data();
+	}
+};
+
+/**
+ * Is basically a struct, but additionally has members for meta() calls. It technically is a Field of a struct, but organized through a struct of fields.
+ */
+template<typename Desc, typename... CellFieldsT, typename Encode>
+class data<kel::lbm::sch::CellFieldStruct<Desc, schema::Struct<CellFieldsT...>>, Encode> final {
+public:
+	using Schema = kel::lbm::sch::CellFieldStruct<Desc,schema::Struct<CellFieldsT...>>;
+	/// @TODO Add MetaSchema to Schema
+	using MetaSchema = typename meta_schema<Schema>::MetaSchema;
+private:
+	static_assert(sizeof...(CellFieldsT) > 0u, "");
+	data<schema::Struct<CellFieldsT...>, Encode> inner_;
+
+	template<uint64_t i>
+	saw::error_or<void> helper_constructor(const data<schema::FixedArray<schema::UInt64,Desc::D>>& grid_size){
+		using MemT = saw::parameter_pack_type<i,CellFieldsT...>::type;
+		{
+			inner_.template get<MemT::KeyLiteral>() = {grid_size};
+		}
+		if constexpr (sizeof...(CellFieldsT) > (i+1u)){
+			return helper_constructor<i+1u>(grid_size);
+		}
+		return saw::make_void();
+	}
+public:
+	data() = delete;
+
+	data(const data<schema::FixedArray<schema::UInt64,Desc::D>>& grid_size__){
+		auto eov = helper_constructor<0u>(grid_size__);
+		(void)eov;
+	}
+
+	const data<MetaSchema, Encode> meta() const {
+		using MemT = saw::parameter_pack_type<0u,CellFieldsT...>::type;
+		return inner_.template get<MemT::KeyLiteral>().meta();
+	}
+
+	template<saw::string_literal Key>
+	auto& get() {
+		return inner_.template get<Key>();
+	}
+
+	template<saw::string_literal Key>
+	const auto& get() const {
+		return inner_.template get<Key>();
+	}
+};
+
+
+}
diff --git a/modules/core/c++/dfs.hpp b/modules/core/c++/dfs.hpp
new file mode 100644
index 0000000..175e9f2
--- /dev/null
+++ b/modules/core/c++/dfs.hpp
@@ -0,0 +1,9 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+// Make a view class which handles the AA push pull pattern, or alternatively the esoteric twist
+}
+}
diff --git a/modules/core/c++/environment.hpp b/modules/core/c++/environment.hpp
new file mode 100644
index 0000000..d8aa9ae
--- /dev/null
+++ b/modules/core/c++/environment.hpp
@@ -0,0 +1,80 @@
+#pragma once
+
+#include <forstio/error.hpp>
+#include <filesystem>
+#include <cstdlib>
+
+namespace kel {
+namespace lbm {
+
+struct environment {
+	std::string name;
+	
+	std::filesystem::path lbm_dir;
+	std::filesystem::path data_dir;
+};
+
+saw::error_or<environment> setup_lbm_env(const std::string_view name){
+	namespace fs = std::filesystem;
+
+	const char* home_dir = std::getenv("HOME");
+	if(not home_dir){
+		return saw::make_error<saw::err::not_found>("Couldn't find home dir");
+	}
+
+	environment env;
+
+	env.lbm_dir = std::filesystem::path{home_dir} / ".lbm";
+
+	{
+		// Set default data location to "${lbm_dir}/data"
+		env.data_dir = env.lbm_dir / "data";
+		// LBM Data Location
+		auto lbm_dir_config = env.lbm_dir / "data_dir_location.txt";
+		if( fs::exists(lbm_dir_config) && fs::is_regular_file(lbm_dir_config) ){
+			std::ifstream file{lbm_dir_config};
+
+			if(file.is_open()){
+				std::stringstream buffer;
+				buffer <<file.rdbuf();
+				env.data_dir = fs::path{buffer.str()};
+			}
+		}
+	}
+
+	return env;
+}
+
+/**
+ * Returns the default output directory.
+ * Located outside the project dir because dispatching build jobs with output data in the git directory
+ * also copies simulated data which takes a long time.
+ */
+saw::error_or<std::filesystem::path> output_directory(){
+	namespace fs = std::filesystem;
+
+	const char* home_dir = std::getenv("HOME");
+	if(not home_dir){
+		return saw::make_error<saw::err::not_found>("Couldn't find home dir");
+	}
+
+	auto lbm_dir = std::filesystem::path{home_dir} / ".lbm";
+
+	{
+		auto lbm_dir_config = lbm_dir / "data_dir_location.txt";
+		if( fs::exists(lbm_dir_config) && fs::is_regular_file(lbm_dir_config) ){
+			std::ifstream file{lbm_dir_config};
+
+			if(file.is_open()){
+				std::stringstream buffer;
+				buffer <<file.rdbuf();
+				return fs::path{buffer.str()};
+			}
+		}
+	}
+
+	return lbm_dir / "data";
+}
+
+}
+}
diff --git a/modules/core/c++/equilibrium.hpp b/modules/core/c++/equilibrium.hpp
new file mode 100644
index 0000000..eb2b043
--- /dev/null
+++ b/modules/core/c++/equilibrium.hpp
@@ -0,0 +1,50 @@
+#pragma once
+
+#include "descriptor.hpp"
+
+namespace kel {
+namespace lbm {
+template<typename T, typename Descriptor>
+saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<sch::Scalar<T>> rho, saw::data<sch::Vector<T,Descriptor::D>> vel){
+	using dfi = df_info<T, Descriptor>;
+
+	saw::data<sch::FixedArray<T,Descriptor::Q>> eq;
+	// Brain broken, here's an owl
+	//   ^
+	//  0.0
+	//  / \
+	//  | |
+	//
+	//  Velocity * Velocity meaning || vel ||_2^2 or <vel,vel>_2
+	saw::data<T> vel_vel{0.0};
+	for(uint64_t j = 0u; j < Descriptor::D; ++j){
+		vel_vel = vel_vel + vel.at({{j}}) * vel.at({{j}});
+	}
+
+	/**
+	 * Calculate equilibrium
+	 */
+	for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+		saw::data<T> vel_c{};
+		for(uint64_t j = 0u; j < Descriptor::D; ++j){
+			// <vel,c_i>_2
+			vel_c = vel_c + (vel.at({{j}}) * saw::data<T>{static_cast<saw::native_data_type<T>::type>(dfi::directions[i][j])});
+		}
+
+		auto vel_c_cs2 = vel_c * saw::data<T>{dfi::inv_cs2};
+
+		eq.at(i).set(
+			dfi::weights[i] * rho.at({}).get() * 
+			(
+			 1.0
+			 + vel_c_cs2.get()
+			 - dfi::inv_cs2 * 0.5 * vel_vel.get()
+			 + vel_c_cs2.get() * vel_c_cs2.get() * 0.5
+			)
+		);
+	}
+
+	return eq;
+}
+}
+}
diff --git a/modules/core/c++/flatten.hpp b/modules/core/c++/flatten.hpp
new file mode 100644
index 0000000..1609589
--- /dev/null
+++ b/modules/core/c++/flatten.hpp
@@ -0,0 +1,42 @@
+#pragma once
+
+#include <forstio/error.hpp>
+#include <forstio/codec/data.hpp>
+
+namespace kel {
+namespace lbm {
+namespace sch {
+using namespace saw::schema;
+}
+
+template<typename T, uint64_t D>
+struct flatten_index {
+public:
+	template<uint64_t i>
+	static constexpr saw::data<sch::UInt64> stride(const saw::data<sch::FixedArray<sch::UInt64,D>>& meta) {
+		if constexpr (i > 0u){
+			return stride<i-1u>(meta) * meta.at({i-1u});
+		}
+
+		return 1u;
+	}
+private:
+	/// 2,3,4 => 2,6,24
+	/// i + j * 2 + k * 3*2
+	/// 1 + 2 * 2 + 3 * 3*2 = 1+4+18 = 23
+	template<uint64_t i>
+	static void apply_i(saw::data<sch::UInt64>& flat_ind, const saw::data<sch::FixedArray<T,D>>& index, const saw::data<sch::FixedArray<T,D>>& meta){
+		if constexpr ( D > i ) {
+			flat_ind = flat_ind + index.at({i}) * stride<i>(meta);
+			apply_i<i+1u>(flat_ind,index,meta);
+		}
+	}
+public:
+	static saw::data<T> apply(const saw::data<sch::FixedArray<T,D>>& index, const saw::data<sch::FixedArray<T,D>>& meta){
+		saw::data<T> flat_ind{0u};
+		apply_i<0u>(flat_ind, index, meta);
+		return flat_ind;
+	}
+};
+}
+}
diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp
new file mode 100644
index 0000000..c974472
--- /dev/null
+++ b/modules/core/c++/fplbm.hpp
@@ -0,0 +1,191 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct FpLbm {};
+struct FpLbmOneParticleNoVelocity{};
+}
+
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::FpLbm, Encode> final {
+private:
+	saw::data<T> relaxation_;
+	saw::data<T> frequency_;
+public:
+	component(typename saw::native_data_type<T>::type relaxation__):
+		relaxation_{{relaxation__}},
+		frequency_{saw::data<T>{1} / relaxation_}
+	{}
+
+	component(const saw::data<T>& relaxation__):
+		relaxation_{relaxation__},
+		frequency_{saw::data<T>{1} / relaxation_}
+	{}
+
+	using Component = cmpt::FpLbm;
+
+
+	template<typename CellFieldSchema, typename MacroFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+	// void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, saw::data<sch::UInt64> time_step){
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+		auto& dfs = dfs_old_f.at(index);
+
+		auto& rho_f = macros.template get<"density">();
+		auto& vel_f = macros.template get<"velocity">();
+
+		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
+		saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);
+
+		compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel);
+		auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+		using dfi = df_info<T,Descriptor>;
+
+		auto& force_f = macros.template get<"force">();
+		auto& force = force_f.at(index);
+
+		auto& porosity_f = macros.template get<"porosity">();
+		auto& porosity = porosity_f.at(index);
+
+		saw::data<sch::Scalar<T>> dfi_inv_cs2;
+		dfi_inv_cs2.at({}).set(dfi::inv_cs2);
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			// saw::data<T> ci_min_u{0};
+			saw::data<sch::Vector<T,Descriptor::D>> ci;
+			for(uint64_t d = 0u; d < Descriptor::D; ++d){
+				ci.at({{d}}).set(static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d]));
+			}
+			auto ci_dot_u = saw::math::dot(ci,vel);
+
+			// saw::data<sch::Vector<T,Descriptor::D>> F_i;
+			// F_i = f * (c_i - u * ics2  + <c_i,u> * c_i * ics2 * ics2) * w_i;
+			saw::data<sch::Scalar<T>> w;
+			w.at({}).set(dfi::weights[i]);
+
+			auto F_i_d = saw::math::dot(force * w, (ci - vel * dfi_inv_cs2 + ci * ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2 ));
+			/*
+			saw::data<sch::Scalar<T>> F_i_sum;
+			for(uint64_t d = 0u; d < Descriptor::D; ++d){
+				saw::data<sch::Scalar<T>> F_i_d;
+				F_i_d.at({}) = F_i.at({{d}});
+				F_i_sum = F_i_sum + F_i_d;
+			}
+			*/
+
+			dfs.at({i}) = dfs.at({i}) + frequency_ * (eq.at(i) - dfs.at({i}) ) + F_i_d.at({}) * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_);
+		}
+	}
+};
+
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::FpLbmOneParticleNoVelocity, Encode> final {
+private:
+	saw::data<T> relaxation_;
+	saw::data<T> frequency_;
+public:
+	component(typename saw::native_data_type<T>::type relaxation__):
+		relaxation_{{relaxation__}},
+		frequency_{saw::data<T>{1} / relaxation_}
+	{}
+
+	component(const saw::data<T>& relaxation__):
+		relaxation_{relaxation__},
+		frequency_{saw::data<T>{1} / relaxation_}
+	{}
+
+	using Component = cmpt::FpLbmOneParticleNoVelocity;
+
+	template<typename CellFieldSchema, typename MacroFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+	// void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, saw::data<sch::UInt64> time_step){
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+		auto& dfs = dfs_old_f.at(index);
+
+		auto& rho_f = macros.template get<"density">();
+		auto& vel_f = macros.template get<"velocity">();
+		auto& por_f = macros.template get<"porosity">();
+
+		// Temporary find a better way
+		auto& force_f = macros.template get<"force">();
+
+
+		/**
+		 * The other methods all work with a flipped porosity.
+		 * For compat reasons this is also flipped
+		 */
+		auto porosity = por_f.at(index);
+		saw::data<sch::Scalar<T>> one;
+		one.at({}) = 1.0;
+		auto flip_porosity = one - porosity;
+
+		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
+
+		saw::data<sch::Scalar<T>> half;
+		half.at({}).set(0.5);
+		saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);// + total_force * ( half / rho );
+
+		compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel);
+		auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+		using dfi = df_info<T,Descriptor>;
+
+		saw::data<sch::Scalar<T>> min_two;
+		min_two.at({}).set(-2);
+
+		// Maybe ? 
+		auto& force = force_f.at(index);
+		force = vel * rho * min_two * flip_porosity;
+
+		saw::data<sch::Scalar<T>> dfi_inv_cs2;
+		dfi_inv_cs2.at({}).set(dfi::inv_cs2);
+
+		// auto vel = vel_f.at(index);
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			// saw::data<T> ci_min_u{0};
+			saw::data<sch::Vector<T,Descriptor::D>> ci;
+			for(uint64_t d = 0u; d < Descriptor::D; ++d){
+				ci.at({{d}}).set(static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d]));
+			}
+			auto ci_dot_u = saw::math::dot(ci,vel);
+
+			// saw::data<sch::Vector<T,Descriptor::D>> F_i;
+			// F_i = f * ((c_i - u) * ics2  + <c_i,u> * c_i * ics2 * ics2) * w_i;
+			saw::data<sch::Scalar<T>> w;
+			w.at({}).set(dfi::weights[i]);
+
+			/*
+			saw::data<sch::Scalar<T>> F_i_sum;
+			for(uint64_t d = 0u; d < Descriptor::D; ++d){
+				saw::data<sch::Scalar<T>> F_i_d;
+				F_i_d.at({}) = F_i.at({{d}});
+				F_i_sum = F_i_sum + F_i_d;
+			}
+			*/
+			auto term1 = (ci-vel) * dfi_inv_cs2;
+			auto term2 = ci * (ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2);
+
+			auto force_projection = saw::math::dot(term1 + term2, force);
+
+			auto F_i = w * force_projection;
+
+			dfs.at({i}) = dfs.at({i})
+									+ frequency_
+										* (eq.at(i) - dfs.at({i}) )
+									+ F_i.at({})
+										* (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_);
+		}
+	}
+
+};
+}
+}
diff --git a/modules/core/c++/geometry.hpp b/modules/core/c++/geometry.hpp
new file mode 100644
index 0000000..1c5d0a3
--- /dev/null
+++ b/modules/core/c++/geometry.hpp
@@ -0,0 +1,21 @@
+#pragma once
+
+#include "component.hpp"
+
+namespace kel {
+namespace lbm {
+/*
+template<typename Schema>
+struct geometry {
+	void apply(const saw::data<Schema>& field, const saw::data<sch::FixedArray<sch::UInt64,2u>>& start, const saw::data<sch::FixedArray<sch::UInt64,2u>>& end, const saw::data<sch::UInt8>& type){
+
+	}
+};
+*/
+
+// Ghost - 0
+// Wall  - 1
+// Fluid - 2
+// Other Stuff - 3-x
+}
+}
diff --git a/modules/core/c++/geometry/poiseulle_channel.hpp b/modules/core/c++/geometry/poiseulle_channel.hpp
new file mode 100644
index 0000000..f719ec4
--- /dev/null
+++ b/modules/core/c++/geometry/poiseulle_channel.hpp
@@ -0,0 +1,22 @@
+#pragma once
+
+#include "../geometry.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct PoiseulleChannel;
+}
+
+template<typename Schema, typename Desc, typename Encode>
+class component<Schema, Desc, cmpt::PoiseulleChannel, Encode> final {
+private:
+public:
+	template<typename CellFieldSchema>
+	void apply(saw::data<CellFieldSchema,Encode>& field, const saw::data<sch::FixedArraysch::UInt64,Desc::D>){
+		auto& info_f = field.template get<"info">();
+	}
+};
+
+}
+}
diff --git a/modules/core/c++/grid.hpp b/modules/core/c++/grid.hpp
new file mode 100644
index 0000000..c9a3b05
--- /dev/null
+++ b/modules/core/c++/grid.hpp
@@ -0,0 +1,34 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+
+/**
+ * I'm mixing up the geometry values regarding inflow etc. And a bit of logic
+ */
+template<typename InfoSchema>
+class domain_registry final {
+private:
+	saw::data<sch::Array<InfoSchema>> infos_;
+public:
+	template<typename T>
+	saw::data<InfoSchema> get_template_id() {
+		static saw::data<InfoSchema> id{std::numeric_limits<typename saw::native_data_type<InfoSchema>::type>::max()};
+
+		if( id == std::numeric_limits<typename saw::native_data_type<InfoSchema>::type>::max() ){
+			auto err_or_id = search_or_register_id(T::name);
+			if(err_or_id.is_error()){
+				// Unsure about recovery from this
+				exit(-1);
+			}
+		}
+	}
+};
+
+template<typename Schema>
+void clean_grid(saw::data<Schema>& info_field){
+}
+}
+}
diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp
new file mode 100644
index 0000000..18c32a5
--- /dev/null
+++ b/modules/core/c++/hlbm.hpp
@@ -0,0 +1,168 @@
+#pragma once
+
+#include "macroscopic.hpp"
+#include "component.hpp"
+#include "equilibrium.hpp"
+
+#include "particle/particle.hpp"
+
+#include <iostream>
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct HlbmInit {}; 
+struct Hlbm {};
+struct HlbmParticle {};
+
+struct HlbmOneParticleMomentumExchange {};
+}
+
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::HlbmInit, Encode> final {
+private:
+	typename saw::native_data_type<T>::type relaxation_;
+	saw::data<T> frequency_;
+public:
+	component(typename saw::native_data_type<T>::type relaxation__):
+		relaxation_{relaxation__},
+		frequency_{typename saw::native_data_type<T>::type(1) / relaxation_}
+	{}
+
+	template<typename CellFieldSchema, typename MacroFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+		auto& porosity_f = macros.template get<"porosity">();
+		auto& particle_N_f = field.template get<"particle_N">();
+		auto& particle_D_f = field.template get<"particle_D">();
+
+		auto& por = porosity_f.at(index);
+		por = {};
+
+		auto& pnf = particle_N_f.at(index);
+		pnf = {};
+
+		auto& pnd = particle_D_f.at(index);
+		pnd = {};
+	}
+};
+
+/**
+ * HLBM collision operator for LBM
+ */
+template<typename T, typename Desc, typename Encode>
+class component<T, Desc, cmpt::Hlbm, Encode> final {
+private:
+	typename saw::native_data_type<T>::type relaxation_;
+	saw::data<T> frequency_;
+public:
+	component(typename saw::native_data_type<T>::type relaxation__):
+		relaxation_{relaxation__},
+		frequency_{typename saw::native_data_type<T>::type(1) / relaxation_}
+	{}
+
+	template<typename CellFieldSchema, typename MacroFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Desc::D>> index, saw::data<sch::UInt64> time_step) const {
+		
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+		auto& particle_N_f = field.template get<"particle_N">();
+		auto& particle_D_f = field.template get<"particle_D">();
+
+		auto& porosity_f = macros.template get<"porosity">();
+		auto& rho_f = macros.template get<"density">();
+		auto& vel_f = macros.template get<"velocity">();
+
+		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
+		saw::data<sch::Vector<T,Desc::D>>& vel = vel_f.at(index);
+
+		compute_rho_u<T,Desc>(dfs_old_f.at(index), rho, vel);
+
+		auto& porosity = porosity_f.at(index);
+		saw::data<sch::Scalar<T>> one;
+		one.at({}) = 1.0;
+		auto flip_porosity = one - porosity;
+
+		auto& N = particle_N_f.at(index);
+		auto& D = particle_D_f.at(index);
+		// Convex combination of velocities
+		vel = vel * porosity + [&]() -> saw::data<sch::Vector<T,Desc::D>> {
+			return (D.at({}).get() > 0.0 ? N * flip_porosity / D : N);
+		}();
+		// Equilibrium
+		auto eq = equilibrium<T,Desc>(rho,vel);
+
+		for(uint64_t i = 0u; i < Desc::Q; ++i){
+			dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i}));
+		}
+
+		porosity.at({}) = 1.0;
+		D.at({}) = 0.0;
+		N = {};
+	}
+};
+
+namespace impl {
+
+}
+
+template<typename T, typename Desc, typename Encode>
+class component<T, Desc, cmpt::HlbmParticle, Encode> final {
+private:
+	template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema, uint64_t i>
+	void apply_i(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& part_groups, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step) const {
+		
+	}	
+public:
+	/*
+	template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema, uint64_t i>
+	void apply_i()
+	*/
+
+	template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& part_group, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step) const {
+		/// Figure out how to access the particle list
+		// auto& p = particles.at(i);
+
+		/// Iterate over the grid bounds
+		// auto& grid = p.template get<"grid">();
+
+		auto& part_spheroid_group = part_group;
+		auto& mvel = macros.template get<"velocity">();
+		{
+			auto& parts = part_spheroid_group.template get<"particles">();
+			auto parts_size = parts.meta().at({0u});
+
+			auto& pi = parts.at(index);
+			auto& pirb = pi.template get<"rigid_body">();
+			auto& pirb_pos = pirb.template get<"position">();
+
+			saw::data<sch::FixedArray<sch::UInt64,Desc::D>> start;
+			saw::data<sch::FixedArray<sch::UInt64,Desc::D>> stop;
+
+			auto aabb = particle_aabb<ParticleSchema>::calculate(part_spheroid_group,index,mvel.meta());
+			/// Ok, I iterate over the space which covers our particle? So lower bounds to upper bounds
+
+			iterator<Desc::D>::apply([&](const auto& index){
+				// ask for the d_k value here.
+				// For every value im iterating over I need sth
+					std::cout<<"Pos: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl;
+			},start,stop);
+
+			// Check
+		}
+		
+
+	}
+};
+
+template<typename T, typename Desc, typename Encode>
+class component<T, Desc, cmpt::HlbmOneParticleMomentumExchange, Encode> final {
+public:
+	template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& part_group, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step) const {
+		//
+	}
+};
+}
+}
diff --git a/modules/core/c++/index.hpp b/modules/core/c++/index.hpp
new file mode 100644
index 0000000..00e597e
--- /dev/null
+++ b/modules/core/c++/index.hpp
@@ -0,0 +1,30 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+
+template<typename T, uint64_t D>
+saw::data<sch::FixedArray<sch::UInt64,D>> lower_index_from_pos(const saw::data<sch::Vector<T,D>>& pos){
+	saw::data<sch::FixedArray<sch::UInt64,D>> ind;
+
+	for(saw::data<sch::UInt64> i{0u}; i < ind.size(); ++i){
+		ind.at(i).set(std::max(pos.at(i).template cast_to<sch::UInt64>().get(),0));
+	}
+
+	return ind;
+}
+
+template<typename T, uint64_t D>
+saw::data<sch::FixedArray<sch::UInt64,D>> upper_index_from_pos(const saw::data<sch::Vector<T,D>>& pos){
+	auto ind = lower_index_from_pos(pos);
+
+	for(saw::data<sch::UInt64> i{0u}; i < ind.size(); ++i){
+		++ind.at(i);
+	}
+	return ind;
+}
+
+}
+}
diff --git a/modules/core/c++/iterator.hpp b/modules/core/c++/iterator.hpp
new file mode 100644
index 0000000..7fd6f58
--- /dev/null
+++ b/modules/core/c++/iterator.hpp
@@ -0,0 +1,99 @@
+#pragma once
+
+#include "descriptor.hpp"
+
+namespace kel {
+namespace lbm {
+template<typename Func>
+void iterate_over(Func&& func, const saw::data<sch::FixedArray<sch::UInt64,2u>>& start, const saw::data<sch::FixedArray<sch::UInt64,2u>>& end, const saw::data<sch::FixedArray<sch::UInt64,2u>>& dist = {{{0u,0u}}}){
+	// static_assert(D == 2u, "Currently a lazy implementation for AND combinations of intervalls.");
+	for(saw::data<sch::UInt64> i{start.at({0u}) + dist.at({0u})}; (i+dist.at({0u})) < end.at({0u}); ++i){
+		for(saw::data<sch::UInt64> j{start.at({1u}) + dist.at({1u})}; (j+dist.at({1u})) < end.at({1u}); ++j){
+			func({{i,j}});
+		}
+	}
+	return;
+}
+
+/**
+ *
+ */
+template<uint64_t D>
+struct iterator {
+private:
+	template<uint64_t i, typename Func>
+	static void iterate_over_i(Func& func,
+		const saw::data<sch::FixedArray<sch::UInt64,D>>& start,
+		const saw::data<sch::FixedArray<sch::UInt64,D>>& end,
+		const saw::data<sch::FixedArray<sch::UInt64,D>>& dist,
+		saw::data<sch::FixedArray<sch::UInt64,D>>& iter
+	){
+		static_assert(i <= D, "Eh. Too tired to think of a good message");
+		if constexpr ( i == D ){
+			func(iter);
+		}else{
+			for(iter.at({i}) = start.at({i})+dist.at({i}); (iter.at({i})+dist.at({i}) < end.at({i})); ++iter.at({i})){
+				iterate_over_i<i+1u,Func>(func,start,end,dist,iter);
+			}
+		}
+	}
+public:
+	template<typename Func>
+	static void apply(
+		Func&& func,
+		const saw::data<sch::FixedArray<sch::UInt64,D>>& start,
+		const saw::data<sch::FixedArray<sch::UInt64,D>>& end,
+		const saw::data<sch::FixedArray<sch::UInt64,D>>& dist = {}
+	){
+		saw::data<sch::FixedArray<sch::UInt64,D>> iter;
+		iterate_over_i<0u,Func>(func, start, end, dist, iter);
+	}
+};
+
+template<typename Tupl>
+struct ct_tuple_iterator;
+
+template<typename... Tup>
+struct ct_tuple_iterator<sch::Tuple<Tup...>> final {
+private:
+	template<typename Func, uint64_t i>
+	static constexpr saw::error_or<void> apply_i(
+		Func& func,
+		saw::data<sch::Tuple<Tup...>>& tup
+	){
+		if constexpr ( i < sizeof...(Tup) ){
+			auto eov = func(tup.template get<i>());
+			if(eov.is_error()){
+				return eov;
+			}
+			return apply_i<Func,i+1u>(func,tup);
+		}
+
+		return saw::make_void();
+	}
+public:
+	template<typename Func>
+	static constexpr saw::error_or<void> apply(
+		Func&& func,
+		saw::data<sch::Tuple<Tup...>>& tup
+	){
+		return apply_i<Func,0u>(func,tup);
+	}
+};
+
+/* Ambiguous
+template<typename Func>
+void iterate_over(Func&& func, const saw::data<sch::FixedArray<sch::UInt64,3u>>& start, const saw::data<sch::FixedArray<sch::UInt64,3u>>& end, const saw::data<sch::FixedArray<sch::UInt64,3u>>& dist = {{{0u,0u,0u}}}){
+	// static_assert(D == 2u, "Currently a lazy implementation for AND combinations of intervalls.");
+	for(saw::data<sch::UInt64> i{start.at({0u}) + dist.at({0u})}; (i+dist.at({0u})) < end.at({0u}); ++i){
+		for(saw::data<sch::UInt64> j{start.at({1u}) + dist.at({1u})}; (j+dist.at({1u})) < end.at({1u}); ++j){
+			for(saw::data<sch::UInt64> k{start.at({2u}) + dist.at({2u})}; (j+dist.at({2u})) < end.at({2u}); ++j){
+				func({{k,j,i}});
+			}
+		}
+	}
+	return;
+}
+*/
+}
+}
diff --git a/modules/core/c++/lbm.hpp b/modules/core/c++/lbm.hpp
new file mode 100644
index 0000000..23d30b1
--- /dev/null
+++ b/modules/core/c++/lbm.hpp
@@ -0,0 +1,65 @@
+#pragma once
+
+#include "args.hpp"
+#include "schema.hpp"
+#include "descriptor.hpp"
+#include "boundary.hpp"
+#include "converter.hpp"
+#include "config.hpp"
+#include "collision.hpp"
+#include "component.hpp"
+#include "environment.hpp"
+#include "equilibrium.hpp"
+#include "flatten.hpp"
+#include "fplbm.hpp"
+#include "chunk.hpp"
+#include "iterator.hpp"
+#include "hlbm.hpp"
+#include "macroscopic.hpp"
+#include "momentum_gather.hpp"
+#include "memory.hpp"
+#include "psm.hpp"
+#include "stream.hpp"
+#include "write_csv.hpp"
+#include "write_vtk.hpp"
+#include "util.hpp"
+
+#include "math/math.hpp"
+
+#include <forstio/codec/unit/unit_print.hpp>
+#include <iostream>
+
+namespace kel {
+namespace lbm {
+template<typename T, typename Desc>
+void print_lbm_meta(
+		const converter<T>& conv,
+		const saw::data<sch::SiKinematicViscosity<T>>& kin_vis_si,
+		const saw::data<sch::SiVelocity<T>>& char_vel,
+		const saw::data<sch::SiMeter<T>>& char_len
+){
+	std::cout
+		<<"[Meta]\n"
+		<<"======\n"
+		<<"Re:     "<<(char_vel * char_len / kin_vis_si)<<"\n"
+		<<"Ma:     "<<(char_vel * saw::data<typename saw::unit_division<sch::Pure<T>, sch::SiVelocity<T>>::Schema>{std::sqrt(df_info<T,Desc>::inv_cs2)})<<"\n"
+		<<"\n"
+		<<"[SI]\n"
+		<<"====\n"
+		<<"Δx:     "<<conv.delta_x()<<"\n"
+		<<"Δt:     "<<conv.delta_t()<<"\n"
+		<<"Δv:     "<<conv.delta_v()<<"\n"
+		<<"Δa:     "<<conv.delta_a()<<"\n"
+		<<"KinVis: "<<kin_vis_si<<"\n"
+		<<"CharV:  "<<char_vel<<"\n"
+		<<"CharL:  "<<char_len<<"\n"
+		<<"\n"
+		<<"[LBM]\n"
+		<<"=====\n"
+		<<"KinVis: "<<conv.kinematic_viscosity_si_to_lbm(kin_vis_si)<<"\n"
+		<<"τ:      "<<(saw::data<typename saw::unit_division<sch::Pure<T>, sch::LbmKinematicViscosity<T>>::Schema >{df_info<T,Desc>::inv_cs2} * conv.kinematic_viscosity_si_to_lbm(kin_vis_si) + saw::data<sch::Pure<T>>{0.5})<<"\n"
+		<<std::endl
+	;
+}
+}
+}
diff --git a/modules/core/c++/lbm_unit.hpp b/modules/core/c++/lbm_unit.hpp
new file mode 100644
index 0000000..2d90652
--- /dev/null
+++ b/modules/core/c++/lbm_unit.hpp
@@ -0,0 +1,70 @@
+#pragma once
+	/**
+	 * Get the conversion parameter with the conversion type
+	 */
+
+#include <forstio/codec/unit/unit.hpp>
+
+#include <string_view>
+
+namespace kel {
+namespace lbm {
+namespace lbm_type {
+struct meter {
+	static constexpr std::string_view name = "meter_lbm";
+	static constexpr std::string_view short_name = "m_lbm";
+};
+
+struct second {
+	static constexpr std::string_view name = "second_lbm";
+	static constexpr std::string_view short_name = "s_lbm";
+};
+}
+
+namespace si_type {
+struct meter {
+	static constexpr std::string_view name = "meter_si";
+	static constexpr std::string_view short_name = "m_si";
+};
+
+struct second {
+	static constexpr std::string_view name = "second_si";
+	static constexpr std::string_view short_name = "s_si";
+};
+}
+
+namespace sch {
+using namespace saw::schema;
+template<typename S>
+using SiMeter = Unit<S, UnitElement<si_type::meter, 1>>;
+
+template<typename S>
+using LbmMeter = Unit<S, UnitElement<lbm_type::meter, 1>>;
+
+template<typename S>
+using SiSecond = Unit<S, UnitElement<si_type::second, 1>>;
+
+template<typename S>
+using LbmSecond = Unit<S, UnitElement<lbm_type::second, 1>>;
+
+template<typename S>
+using SiVelocity = Unit<S, UnitElement<si_type::meter, 1>, UnitElement<si_type::second, -1>>;
+
+template<typename S>
+using LbmVelocity = Unit<S, UnitElement<lbm_type::meter, 1>, UnitElement<lbm_type::second, -1>>;
+
+template<typename S>
+using SiAcceleration = Unit<S, UnitElement<si_type::meter, 1>, UnitElement<si_type::second, -2>>;
+
+template<typename S>
+using LbmAcceleration = Unit<S, UnitElement<lbm_type::meter, 1>, UnitElement<lbm_type::second, -2>>;
+
+template<typename S>
+using SiKinematicViscosity = Unit<S, UnitElement<si_type::meter, 2>, UnitElement<si_type::second, -1>>;
+
+template<typename S>
+using LbmKinematicViscosity = Unit<S, UnitElement<lbm_type::meter, 2>, UnitElement<lbm_type::second, -1>>;
+
+}
+}
+}
diff --git a/modules/core/c++/macroscopic.hpp b/modules/core/c++/macroscopic.hpp
new file mode 100644
index 0000000..19ab3e9
--- /dev/null
+++ b/modules/core/c++/macroscopic.hpp
@@ -0,0 +1,36 @@
+#pragma once
+
+#include "descriptor.hpp"
+
+namespace kel {
+namespace lbm {
+/**
+ * Calculate the macroscopic variables rho and u in Lattice Units.
+ */
+template<typename T, typename Desc>
+void compute_rho_u (
+		const saw::data<sch::FixedArray<T,Desc::Q>>& dfs,
+		saw::ref<saw::data<sch::Scalar<T>>> rho,
+		saw::ref<saw::data<sch::Vector<T,Desc::D>>> vel
+	)
+{
+	using dfi = df_info<T, Desc>;
+
+	rho().at({}).set(0);
+	for(uint64_t i = 0; i < Desc::D; ++i){
+		vel().at({{i}}).set(0);
+	}
+
+	for(size_t j = 0; j < Desc::Q; ++j){
+		rho().at({}) = rho().at({}) + dfs.at({j});
+		for(size_t i = 0; i < Desc::D; ++i){
+			vel().at({{i}}) = vel().at({{i}}) + saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[j][i])} * dfs.at({j});
+		}
+	}
+
+	for(size_t i = 0; i < Desc::D; ++i){
+		vel().at({{i}}) = vel().at({{i}}) / rho().at({});
+	}
+}
+}
+}
diff --git a/modules/core/c++/math/math.hpp b/modules/core/c++/math/math.hpp
new file mode 100644
index 0000000..3920bec
--- /dev/null
+++ b/modules/core/c++/math/math.hpp
@@ -0,0 +1,4 @@
+#pragma once
+
+#include "n_linear.hpp"
+#include "n_closest.hpp"
diff --git a/modules/core/c++/math/n_closest.hpp b/modules/core/c++/math/n_closest.hpp
new file mode 100644
index 0000000..ac0fe2f
--- /dev/null
+++ b/modules/core/c++/math/n_closest.hpp
@@ -0,0 +1,54 @@
+#pragma once
+
+#include "../common.hpp"
+#include "../iterator.hpp"
+
+namespace kel {
+namespace lbm {
+
+template<typename FieldSchema, typename Encode, typename T, uint64_t D>
+saw::data<typename FieldSchema::StoredValueSchema> n_closest_read(const saw::data<sch::Ptr<FieldSchema>,Encode>& f, const saw::data<sch::Vector<T,D>>& frac_ind){
+
+	auto shift_frac_ind = frac_ind;
+	for(uint64_t i{0u}; i < D; ++i){
+		
+		shift_frac_ind.at({{i}}) = shift_frac_ind.at({{i}})  + saw::data<T>{0.5};
+		if(shift_frac_ind.at({{i}}).get() < 0){
+			shift_frac_ind.at({{i}}) = {};
+		}
+	}
+
+	saw::data<sch::FixedArray<sch::UInt64,D>> shift_ind;
+	for(uint64_t i{0u}; i < D; ++i){
+		shift_ind.at({i}) = frac_ind.at({{i}}).template cast_to<sch::UInt64>();
+	}
+
+	return f.at(shift_ind);
+}
+
+template<typename FieldSchema, typename Encode, typename T, uint64_t D>
+void n_closest_add(const saw::data<sch::Ptr<FieldSchema>,Encode>& f, const saw::data<sch::Vector<T,D>>& frac_ind, const saw::data<typename FieldSchema::StoredValueSchema>& val){
+	auto shift_frac_ind = frac_ind;
+	for(uint64_t i{0u}; i < D; ++i){
+		
+		shift_frac_ind.at({{i}}) = shift_frac_ind.at({{i}})  + saw::data<T>{0.5};
+		if(shift_frac_ind.at({{i}}).get() < 0){
+			shift_frac_ind.at({{i}}) = {};
+		}
+	}
+
+	auto f_meta = f.meta();
+	saw::data<sch::FixedArray<sch::UInt64,D>> shift_ind;
+	for(uint64_t i{0u}; i < D; ++i){
+		shift_ind.at({i}) = frac_ind.at({{i}}).template cast_to<sch::UInt64>();
+		if(shift_ind.at({i}) < f_meta.at({i})){
+			shift_ind.at({i}) = f_meta.at({i}) - 1u;
+		}
+	}
+	auto& f_i = f.at(shift_ind);
+
+	f_i = f_i + val;
+}
+
+}
+}
diff --git a/modules/core/c++/math/n_linear.hpp b/modules/core/c++/math/n_linear.hpp
new file mode 100644
index 0000000..b378440
--- /dev/null
+++ b/modules/core/c++/math/n_linear.hpp
@@ -0,0 +1,139 @@
+#pragma once
+
+#include "../common.hpp"
+#include "../iterator.hpp"
+
+namespace kel {
+namespace lbm {
+namespace impl {
+template<typename FieldSchema, typename T, uint64_t D>
+struct n_linear_interpolate_helper final {
+	template<uint64_t i = 0u>
+	auto apply(const saw::data<FieldSchema>& field, const saw::data<sch::Vector<T,D>>& pos){
+		return pos;
+	}
+};
+}
+
+template<typename T, uint64_t D>
+saw::data<sch::Tuple<sch::Vector<sch::UInt64,D>,sch::Vector<T,D>>> position_to_index_and_fraction(const saw::data<sch::Vector<T,D>>& pos){
+	saw::data<sch::Tuple<sch::Vector<sch::UInt64,D>,sch::Vector<T,D>>> sep;
+
+	auto& ind = sep.template get<0u>();
+	auto& frac = sep.template get<1u>();
+
+	auto pos_cpy = pos;
+	// Guarantee that the pos is at least 0
+	for(uint64_t i = 0u; i < D; ++i){
+		pos_cpy.at({{i}}).set(std::max(pos.at({{i}}).get(), static_cast<typename saw::native_data_type<T>::type>(0)));
+	}
+
+	// Now we can cast to uint64_t
+	for(uint64_t i = 0u; i < D; ++i){
+		ind.at({{i}}) = pos_cpy.at({{i}}).template cast_to<sch::UInt64>();
+	}
+
+	frac = pos_cpy - ind.template cast_to<T>();
+
+	return sep;
+}
+
+template<typename T, uint64_t D>
+auto floor_index_from_position(const saw::data<sch::Vector<T,D>>& pos){
+	return position_to_index_and_fraction(pos).template get<0u>();
+}
+
+template<typename T, uint64_t D>
+saw::data<sch::Tuple<sch::Vector<sch::UInt64,D>,sch::Vector<T,D>>> position_to_index_and_fraction_bounded(
+	const saw::data<sch::Vector<T,D>>& pos,
+	const saw::data<sch::Vector<sch::UInt64,D>>& bound)
+{
+	auto infr = position_to_index_and_fraction(pos);
+	auto& ind = infr.template get<0u>();
+	auto& fra = infr.template get<1u>();
+	for(uint64_t i = 0u; i < D; ++i){
+		// If index is higher than bound. Set to bound and reset fraction
+		if((ind.at({{i}}).get()+1u) >= bound.at({{i}}).get()){
+			ind.at({{i}}).set(bound.at({{i}}).get()-1u);
+			fra.at({{i}}) = {};
+		}
+	}
+	return infr;
+}
+
+
+template<typename FieldSchema, typename T, uint64_t D>
+auto n_linear_interpolate(
+		const saw::data<FieldSchema>& field, const saw::data<sch::Vector<T,D>>& pos){
+
+	// Pos
+	auto pos_bound = pos;
+	
+	// Dimensions
+	auto meta = field.dims();
+
+	// Lower Index
+	saw::data<sch::FixedArray<sch::UInt64,D>> ind;
+
+	for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){
+		// Native Positive i
+		auto npos_i = pos.at({i}).get();
+
+		{
+			// Ok I want to stay in bounds
+			npos_i = std::min(npos_i,meta.at(i).get()-1.0);
+			npos_i = std::max(npos_i,1.0);
+		}
+
+		// Native Index i
+		auto nind_i = static_cast<uint64_t>(std::floor(npos_i))-1ul;
+
+		// Set index to i
+		ind.at(i).set(nind_i);
+	}
+	saw::data<sch::Vector<T,D>> pos_frac;
+	for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){
+		pos_frac.at({i}) = pos_bound.at({i}) - ind.at(i).template cast_to<T>();
+	}
+
+	// Base value
+	saw::data<typename FieldSchema::ValueType> res;
+
+	// constexpr uint64_t d_corners = 1ul << D;
+
+	saw::data<sch::FixedArray<sch::UInt64,D>> ones_ind;
+	for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){
+		ones_ind.at({i}).set(1u);
+	}
+
+	iterator<D>::apply([&](auto ind){
+		// Iterates over (0,0,0) to (1,1,1)
+		saw::data<T> weight{1.0};
+
+		for(saw::data<sch::UInt64> d{0u}; d < saw::data<sch::UInt64>{D}; ++d){
+
+			saw::data<T> t = pos_frac.at({d});
+
+			if(ind.at(d).get() == 0u){
+				weight = weight * (saw::data<T>{1} - t);
+			}else{
+				weight = weight * t;
+			}
+		}
+
+	}, {}, ones_ind);
+
+	/// TODO I need to actually calc stuff
+	return field.at({});
+}
+
+template<typename FieldSchema, typename T>
+saw::data<sch::Vector<T,2u>> bilinear_interpolate(const saw::data<FieldSchema>& field, const saw::data<sch::Vector<T,2u>>& pos){
+	saw::data<sch::Vector<T,2u>> res;
+
+	{
+	}
+	return {};
+}
+}
+}
diff --git a/modules/core/c++/math/round.hpp b/modules/core/c++/math/round.hpp
new file mode 100644
index 0000000..d3a2586
--- /dev/null
+++ b/modules/core/c++/math/round.hpp
@@ -0,0 +1,26 @@
+#pragma once
+
+#include "../common.hpp"
+
+namespace kel {
+namespace lbm {
+
+template<typename T, uint64_t D>
+saw::data<sch::FixedArray<sch::UInt64,D>> round_to_unsigned(const saw::data<sch::Vector<T,D>>& inp){
+	saw::data<sch::FixedArray<sch::UInt64,D>> rv;
+				
+	auto zero = static_cast<saw::native_data_type<T>::type>(0);
+	auto half = static_cast<saw::native_data_type<T>::type>(0.5);
+
+	for(uint64_t i{0u}; i < D; ++i){
+		auto val = inp.at({{i}}).get()+half;
+		val = std::max(zero,val);
+
+		rv.at({i}).set(static_cast<uint64_t>(val));
+	}
+
+	return rv;
+}
+
+}
+}
diff --git a/modules/core/c++/memory.hpp b/modules/core/c++/memory.hpp
new file mode 100644
index 0000000..e97c7fc
--- /dev/null
+++ b/modules/core/c++/memory.hpp
@@ -0,0 +1,82 @@
+#pragma once
+
+#include "common.hpp"
+#include "chunk.hpp"
+
+namespace kel {
+namespace lbm {
+
+namespace impl {
+
+template<typename Sch>
+struct memory_size_helper;
+
+template<typename Sch, uint64_t N>
+struct memory_size_helper<sch::Primitive<Sch,N>> {
+	static constexpr uint64_t bytes = N;
+};
+
+template<typename Sch, uint64_t... N>
+struct memory_size_helper<sch::FixedArray<Sch,N...>> {
+	static constexpr uint64_t bytes = saw::ct_multiply<uint64_t,N...>::value * memory_size_helper<Sch>::bytes;
+};
+
+template<typename Sch, uint64_t... N>
+struct memory_size_helper<sch::Tensor<Sch,N...>> {
+	static constexpr uint64_t bytes = saw::ct_multiply<uint64_t,N...>::value * memory_size_helper<Sch>::bytes;
+};
+
+template<typename Sch, uint64_t Ghost, uint64_t... N>
+struct memory_size_helper<sch::Chunk<Sch,Ghost,N...>> {
+	static constexpr uint64_t bytes = memory_size_helper<typename sch::Chunk<Sch,Ghost,N...>::InnerSchema>::bytes;
+};
+
+template<typename... Members>
+struct memory_size_helper<sch::Struct<Members...>> {
+
+	template<uint64_t i>
+	static constexpr uint64_t apply_i() {
+		if constexpr ( i < sizeof...(Members) ){
+			using M_I = typename saw::parameter_pack_type<i,Members...>::type;
+			return apply_i<i+1u>() + memory_size_helper<typename M_I::ValueType>::bytes;
+		}
+		return 0u;
+	}
+
+	static constexpr uint64_t bytes = apply_i<0u>();
+};
+
+template<typename... T>
+class memory_estimate_helper final {
+private:
+
+	template<uint64_t i>
+	static void apply_i(saw::data<sch::UInt64>& bytes){
+
+		if constexpr ( i < sizeof...(T)){
+			using T_I = typename saw::parameter_pack_type<i,T...>::type;
+
+			bytes.set(bytes.get() + memory_size_helper<T_I>::bytes);
+
+			apply_i<i+1u>(bytes);
+		}
+	}
+
+public:
+
+	static void apply(saw::data<sch::UInt64>& bytes){
+		apply_i<0u>(bytes);
+	}
+};
+}
+
+template<typename... T>
+saw::data<sch::UInt64> memory_estimate(){
+	saw::data<sch::UInt64> bytes;
+
+	impl::memory_estimate_helper<T...>::apply(bytes);
+
+	return bytes;
+}
+}
+}
diff --git a/modules/core/c++/momentum_gather.hpp b/modules/core/c++/momentum_gather.hpp
new file mode 100644
index 0000000..7b1fc96
--- /dev/null
+++ b/modules/core/c++/momentum_gather.hpp
@@ -0,0 +1,54 @@
+#pragma once
+
+#include "common.hpp"
+#include "macroscopic.hpp"
+#include "component.hpp"
+#include "equilibrium.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+// Gather the force, Luke!
+struct ForceGather {};
+}
+
+template<typename T, typename Descriptor, typename Encode>
+class component<T,Descriptor,cmpt::ForceGather, Encode> {
+private:
+public:
+	template<typename CellFieldSchema, typename MacroFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros,
+		saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+
+		using dfi = df_info<T,Descriptor>;
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+
+		auto& porous_f = macros.template get<"porosity">();
+		auto& porous = porous_f.at(index);
+
+		auto& force_f = macros.template get<"force">();
+		auto& force = force_f.at(index);
+
+		for(uint64_t k{0u}; k < Descriptor::D; ++k){
+			force.at({{k}}).set(0);
+		}
+
+		auto& dfs = dfs_old_f.at(index);
+		for(uint64_t i{0u}; i < Descriptor::Q; ++i){
+			uint64_t i_opp = dfi::opposite_index[i];
+			auto dfs_diff = dfs.at({i}) - dfs.at({i_opp});
+			for(uint64_t k{0u}; k < Descriptor::D; ++k){
+				force.at({{k}}) = force.at({{k}}) + dfs_diff * saw::data<T>{dfi::directions[i][k]};
+			}
+		}
+
+		saw::data<sch::Scalar<T>> one;
+		one.at({}) = 1.0;
+		auto flip_porous = one - porous;
+		force = force * flip_porous;
+	}
+};
+}
+}
diff --git a/modules/core/c++/particle.hpp b/modules/core/c++/particle.hpp
new file mode 100644
index 0000000..691a74b
--- /dev/null
+++ b/modules/core/c++/particle.hpp
@@ -0,0 +1,37 @@
+#pragma once
+
+#include "component.hpp"
+#include "particle/particle.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct Particle {};
+}
+
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::Particle, Encode> {
+private:
+	saw::data<sch::Scalar<T>> dt_;
+public:
+	component(saw::data<sch::Scalar<T>> dt__):
+		dt_{dt__}
+	{}
+
+	template<typename ParticleSchema, typename MacroFieldSchema>
+	void apply(const saw::data<ParticleSchema, Encode>& particles, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::UInt64> index, saw::data<sch::UInt64> time_step) const {
+
+		auto& p = particles.at(index);
+
+		// Compute forces
+
+		// Update particle velocity
+		verlet_step_lambda<T,Descriptor::D>(p,{1.0});  
+
+		// Update porosity over lattice nodes
+	
+		
+	}
+};
+}
+}
diff --git a/modules/core/c++/particle/aabb.hpp b/modules/core/c++/particle/aabb.hpp
new file mode 100644
index 0000000..1773dea
--- /dev/null
+++ b/modules/core/c++/particle/aabb.hpp
@@ -0,0 +1,55 @@
+#pragma once
+
+#include "common.hpp"
+#include "schema.hpp"
+
+namespace kel {
+namespace lbm {
+template<typename PGroup>
+class particle_aabb final {
+	static_assert(saw::always_false<PGroup>, "Not supported");
+};
+
+template<typename T, uint64_t D>
+class particle_aabb<
+	sch::ParticleGroup<T,D,coll::Spheroid<T>>
+> final {
+public:
+	using Schema = sch::ParticleGroup<T,D,coll::Spheroid<T>>;
+
+	using AABB = sch::Struct<
+		sch::Member<sch::FixedArray<sch::UInt64,D>,"a">,
+		sch::Member<sch::FixedArray<sch::UInt64,D>,"b">
+	>;
+public:
+	template<typename Encode>
+	static constexpr saw::data<AABB> calculate(const saw::data<Schema,Encode>& p_grp, const saw::data<sch::FixedArray<sch::UInt64,1u>>& i, const saw::data<sch::FixedArray<sch::UInt64,D>>& meta){
+
+		saw::data<AABB> aabb;
+		auto& parts = p_grp.template get<"particles">();
+
+		auto& pi = parts.at(i);
+		auto& pirb = pi.template get<"rigid_body">();
+		auto& pirb_pos = pirb.template get<"position">();
+
+		auto& a = aabb.template get<"a">();
+		auto& b = aabb.template get<"b">();
+
+		const saw::data<sch::Scalar<T>>& rad_d = p_grp.template get<"collision">().template get<"radius">().at({0u});
+
+		saw::data<sch::Vector<T,D>> lower;
+		saw::data<sch::Vector<T,D>> upper;
+
+		for(uint64_t i{0u}; i < D; ++i){
+			lower.at({{i}}) = pirb_pos.at({{i}}) >= rad_d.at({}) ? (pirb_pos.at({{i}}) - rad_d.at({})) : saw::data<T>{0};
+			a.at({i}) = lower.at({{i}}).template cast_to<sch::UInt64>();
+			upper.at({{i}}) = pirb_pos.at({{i}}) + rad_d.at({});
+			b.at({i}) = (upper.at({{i}})+saw::data<T>{1}).template cast_to<sch::UInt64>();
+		}
+		
+		return aabb;
+
+	}
+};
+}
+}
diff --git a/modules/core/c++/particle/blur.hpp b/modules/core/c++/particle/blur.hpp
new file mode 100644
index 0000000..b7a1988
--- /dev/null
+++ b/modules/core/c++/particle/blur.hpp
@@ -0,0 +1,37 @@
+#pragma once
+
+namespace kel {
+namespace lbm {
+template<typename T, uint64_t D>
+void blur_mask(saw::data<sch::Array<T,D>>& p_mask){
+
+	saw::data<T> mid{2.0/3.0};
+	saw::data<T> edge{1.0/6.0};
+
+	// saw::data<sch::FixedArray<sch::UInt64>,2u> blur_weights{{2.0/3.0},{1.0/6.0}};
+
+	auto meta = p_mask.dims();
+	saw::data<sch::Array<T,D>> blurred_mask{meta};
+
+	/* 1D blur into N-D Blur*/
+	for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){
+		iterator<D>::apply([&](const auto& index){
+			blurred_mask.at(index) = p_mask.at(index) * mid;
+
+			if(index.at({i}).get() > 0u){
+				auto l_ind = index;
+				l_ind.at({i}) = ind.at({i}) - 1u;
+				blurred_mask.at(index) = blurred_mask.at(index) + p_mask.at(l_ind) * edge;
+			}
+			if((index.at({i}).get() + 1u) < meta.at({i})){
+				auto r_ind = index;
+				r_ind.at({i}) = ind.at({i}) + 1u;
+				blurred_mask.at(index) = blurred_mask.at(index) + p_mask.at(r_ind) * edge;
+			}
+		},{},meta);
+
+		p_mask = blurred_mask;
+	}
+}
+}
+}
diff --git a/modules/core/c++/particle/common.hpp b/modules/core/c++/particle/common.hpp
new file mode 100644
index 0000000..9e673c2
--- /dev/null
+++ b/modules/core/c++/particle/common.hpp
@@ -0,0 +1,3 @@
+#pragma once
+
+#include "../common.hpp"
diff --git a/modules/core/c++/particle/geometry/circle.hpp b/modules/core/c++/particle/geometry/circle.hpp
new file mode 100644
index 0000000..cf9b87e
--- /dev/null
+++ b/modules/core/c++/particle/geometry/circle.hpp
@@ -0,0 +1,85 @@
+#pragma once
+
+#include "../particle.hpp"
+
+namespace kel {
+namespace lbm {
+
+template<typename T, uint64_t D>
+class particle_circle_geometry {
+private:
+public:
+	particle_circle_geometry()
+	{}
+
+	template<typename MT = T>
+	saw::data<sch::ParticleMask<MT,D>> generate_mask(uint64_t resolution, uint64_t boundary_nodes = 0) const {
+		saw::data<sch::ParticleMask<MT,D>> mask;
+
+		auto& grid = mask.template get<"grid">();
+		auto& com = mask.template get<"center_of_mass">();
+		com = {};
+
+		//uint64_t rad_i = static_cast<uint64_t>(resolution * radius_.get())+1u;
+		uint64_t diameter_i = resolution;
+		uint64_t size = diameter_i + 2*boundary_nodes;
+		grid = {{{size,size}}};
+
+		saw::data<T> delta_x{static_cast<typename saw::native_data_type<T>::type>(2.0  / static_cast<double>(diameter_i))};
+		saw::data<T> center = (saw::data<T>{1.0} + saw::data<T>{2.0} * saw::data<T>{static_cast<saw::native_data_type<T>::type>(boundary_nodes)/diameter_i});
+
+		saw::data<sch::FixedArray<sch::UInt64,D>> boundary_arr;
+		for(uint64_t i = 0u; i < D; ++i){
+			boundary_arr.set(boundary_nodes);
+		}
+
+		iterator<D>::apply([&](const auto& index){
+			grid.at(index).set(0);
+		}, {}, grid.dims());
+		
+		iterator<D>::apply([&](const auto& index){
+			saw::data<sch::Vector<T,D>> f_ind;
+			for(uint64_t i = 0u; i < D; ++i){
+				f_ind.at({{i}}) = (index.at({{i}}) + 0.5) * delta_x - center;
+			}
+
+			auto norm_f_ind = saw::math::norm(f_ind);
+		}, {}, grid.dims(),boundary_arr);
+
+		for(uint64_t i = 0; i < size; ++i){
+			for(uint64_t j = 0; j < size; ++j){
+				grid.at({{i,j}}).set(0);
+				if(i >= boundary_nodes and j >= boundary_nodes and i < (diameter_i + boundary_nodes) and j < (diameter_i + boundary_nodes) ){
+					saw::data<T> fi = saw::data<T>{static_cast<saw::native_data_type<T>::type>(0.5+i)} * delta_x - center;
+					saw::data<T> fj = saw::data<T>{static_cast<saw::native_data_type<T>::type>(0.5+j)} * delta_x - center;
+
+					auto norm_f_ij = fi*fi + fj*fj;
+					if(norm_f_ij.get() <= 1){
+						grid.at({{i,j}}).set(1);
+					}
+				}
+			}
+		}
+
+		// I don't fully remember what this did
+		saw::data<T> total_mass{};
+		iterator::apply<D>::apply([&](const saw::data<sch::FixedArray<sch::UInt64,D>>& index){
+			auto ind_vec = saw::math::vectorize_data(index).template cast_to<T>();
+			for(uint64_t i{0u}; i < D; ++i){
+				ind_vec.at({{i}}) = ind_vec.at({{i}}) * grid.at(index);
+			}
+			com = com + ind_vec;
+
+			total_mass = total_mass + grid.at(index);
+		},{}, grid.dims());
+
+		for(uint64_t i{0u}; i < D; ++i){
+			com.at({{i}}) = com.at({{i}}) / total_mass;
+		}
+
+		return mask;
+	}
+};
+
+}
+}
diff --git a/modules/core/c++/particle/geometry/cube.hpp b/modules/core/c++/particle/geometry/cube.hpp
new file mode 100644
index 0000000..6392de8
--- /dev/null
+++ b/modules/core/c++/particle/geometry/cube.hpp
@@ -0,0 +1,22 @@
+#pragma once
+
+#include "../particle.hpp"
+
+namespace kel {
+namespace lbm {
+template<typename T, uint64_t D>
+class particle_cubic_geometry {
+private:
+public:
+	template<typename MT = T>
+	saw::data<sch::ParticleMask<MT,D>> generate_mask(uint64_t resolution, uint64_t bound_nodes = 0u){
+
+		saw::data<sch::ParticleMask<MT,D>> mask;
+
+		auto& grid = mask.template get<"grid">();
+		auto& com = mask.template get<"center_of_mass">();
+
+
+	}
+}
+}
diff --git a/modules/core/c++/particle/particle.hpp b/modules/core/c++/particle/particle.hpp
new file mode 100644
index 0000000..8e75e5a
--- /dev/null
+++ b/modules/core/c++/particle/particle.hpp
@@ -0,0 +1,246 @@
+#pragma once
+
+#include <forstio/codec/math.hpp>
+#include <forstio/codec/data_math.hpp>
+#include <forstio/codec/data.hpp>
+
+#include "../iterator.hpp"
+
+#include "schema.hpp"
+#include "aabb.hpp"
+#include "particle_opa.hpp"
+
+namespace kel {
+namespace lbm {
+
+template<typename T, uint64_t D>
+saw::data<sch::ParticleGroup<T,D, coll::Spheroid<T>>> create_spheroid_particle_group(
+		saw::data<sch::Scalar<T>> radius_p,
+		saw::data<sch::Scalar<T>> density_p,
+		const saw::data<sch::UInt64>& mask_resolution
+){
+	saw::data<sch::ParticleGroup<T,D,coll::Spheroid<T>>> part;
+
+	auto& rad_s = part.template get<"collision">().at({0u}).template get<"radius">();
+	rad_s = radius_p;
+
+	auto& mask = part.template get<"mask">();
+	auto& density = part.template get<"density">().at({{0u}});
+
+	auto& total_mass = part.template get<"total_mass">().at({{0u}});
+	// Paranoia
+	total_mass.at({}) = {};
+
+	static_assert(D >= 1u and D <= 3u, "Dimensions only supported for Dim 1,2 & 3.");
+	density = density_p;
+
+	saw::data<sch::FixedArray<sch::UInt64,D>> mask_dims;
+	for(uint64_t i = 0u; i < D; ++i){
+		mask_dims.at({i}) = mask_resolution;
+	}
+	saw::data<T> rad_d = radius_p.at({});
+	saw::data<T> dia_d = rad_d * 2;
+
+	mask = {mask_dims};
+
+	auto& mask_step = part.template get<"mask_step">().at({{0u}});
+	mask_step.at({}) = dia_d / mask_resolution.template cast_to<T>();
+
+	auto& com = part.template get<"center_of_mass">().at({{0u}});
+	// Paranoia
+	for(uint64_t i = 0u; i < D; ++i){
+		com.at({{i}}) = {};
+	}
+	saw::data<sch::UInt64> ele_ctr{0u};
+
+	// Radius ^ 2
+	saw::data<sch::Scalar<T>> rad_2_d;
+	rad_2_d.at({}) = rad_d * rad_d;
+
+	saw::data<sch::Vector<T,D>> center;
+	for(uint64_t i = 0u; i < D; ++i){
+		center.at({{i}}) = rad_d;
+	}
+
+	iterator<D>::apply([&](const auto& index){
+		++ele_ctr;
+
+		saw::data<sch::Vector<T,D>> offset_index = saw::math::vectorize_data(index).template cast_to<T>() - center;
+
+		auto& dpi = mask.at(index);
+
+		for(uint64_t i = 0u; i < D; ++i){
+			com.at({{i}}) = com.at({{i}}) + index.at({i}).template cast_to<T>() * dpi;
+		}
+
+		total_mass.at({}) = total_mass.at({}) + dpi;
+
+	},{},mask_dims);
+
+	for(uint64_t i = 0u; i < D; ++i){
+		com.at({{i}}) = com.at({{i}}) / total_mass.at({});
+	}
+	return part;
+}
+
+/*
+template<typename T, uint64_t D>
+saw::data<sch::Particle<T,D, sch::ParticleCollisionSpheroid<T>>> create_spheroid_particle(
+		saw::data<sch::Vector<T,D>> pos_p,
+		saw::data<sch::Vector<T,D>> vec_p,
+		saw::data<sch::Vector<T,D>> acc_p,
+		saw::data<sch::Vector<T,D>> rot_pos_p,
+		saw::data<sch::Vector<T,D>> rot_vel_p,
+		saw::data<sch::Vector<T,D>> rot_acc_p,
+		saw::data<sch::Scalar<T>> rad_p,
+		saw::data<sch::Scalar<T>> density_p,
+		saw::data<sch::Scalar<T>> dt
+	){
+
+	saw::data<sch::Particle<T,D>> part;
+	auto& body = part.template get<"rigid_body">();
+	auto& mass = part.template get<"mass">();
+	
+	auto& pos = body.template get<"position">();
+	auto& pos_old = body.template get<"position_old">();
+	auto& acc = body.template get<"acceleration">();
+
+	auto& rot = body.template get<"rotation">();
+	auto& rot_old = body.template get<"rotation_old">();
+
+	auto& coll = part.template get<"collision">();
+	auto& rad = coll.template get<"radius">();
+
+	pos = pos_p;
+	pos_old = pos - vec_p * dt;
+	acc = acc_p;
+	rad = rad_p;
+
+	if constexpr ( D == 1u){
+		saw::data<sch::Scalar<T>> c;
+		c.at({}).set(2.0);
+		mass = rad_p * c * density_p;
+	} else if constexpr ( D == 2u){
+		saw::data<sch::Scalar<T>> pi;
+		pi.at({}).set(3.14159);
+		mass = rad_p * rad_p * pi * density_p;
+	} else if constexpr ( D == 3u ){
+		saw::data<sch::Scalar<T>> c;
+		c.at({}).set(3.14159 * 4.0 / 3.0);
+		mass = rad_p * rad_p * rad_p * c * density_p;
+	} else {
+		static_assert(D == 0u or D > 3u, "Dimensions only supported for Dim 1,2 & 3.");
+	}
+
+	return part;
+}
+*/
+
+template<typename T,uint64_t D>
+constexpr auto verlet_step_lambda = [](saw::data<sch::Particle<T,D>>& particle, saw::data<sch::Scalar<T>> time_step_delta){
+	auto& body = particle.template get<"rigid_body">();
+
+	auto& pos = body.template get<"position">();
+	auto& pos_old = body.template get<"position_old">();
+
+	auto& pos_acc = body.template get<"acceleration">();
+	
+	auto& rot = body.template get<"rotation">();
+	auto& rot_old = body.template get<"rotation_old">();
+
+	auto& rot_acc = body.template get<"angular_acceleration">();
+
+	auto tsd_squared = time_step_delta * time_step_delta;
+
+	saw::data<sch::Vector<T,D>> pos_new;
+	// Actual step
+	saw::data<sch::Scalar<T>> two;
+	two.at({}).set(2.0);
+	pos_new = pos * two - pos_old + pos_acc * tsd_squared;
+
+	// Angular
+	saw::data<typename sch::impl::rotation_type_helper<T,D>::Schema> rot_new;
+	rot_new = rot * two - rot_old + rot_acc * tsd_squared;
+
+	// Swap - Could be std::swap?
+	pos_old = pos;
+	pos = pos_new;
+
+	rot_old = rot;
+	rot = rot_new;
+};
+
+template<typename T, uint64_t D>
+constexpr auto handle_collision = [](
+		saw::data<sch::Particle<T,D>>& left, const saw::data<sch::Scalar<T>>& mass_l,
+		saw::data<sch::Particle<T,D>>& right, const saw::data<sch::Scalar<T>>& mass_r,
+		saw::data<sch::Vector<T,D>> unit_pos_rel, saw::data<sch::Vector<T,D>> vel_rel,
+		saw::data<sch::Scalar<T>> d_t
+){
+	auto& rb_l = left.template get<"rigid_body">();
+	auto& pos_l = rb_l.template get<"position">();
+	auto& pos_old_l = rb_l.template get<"position_old">();
+	auto vel_l = (pos_l-pos_old_l)/d_t;
+
+	auto& rb_r = right.template get<"rigid_body">();
+	auto& pos_r = rb_r.template get<"position">();
+	auto& pos_old_r = rb_r.template get<"position_old">();
+	auto vel_r = (pos_r-pos_old_r)/d_t;
+
+	auto vel_pos_rel_dot = saw::math::dot(unit_pos_rel,vel_rel);
+
+  if( vel_pos_rel_dot.at({{0u}}).get() < 0.0 ){
+		pos_l = pos_l + vel_rel * unit_pos_rel * d_t;
+		pos_r = pos_r - vel_rel * unit_pos_rel * d_t;
+	}
+};
+
+
+template<typename T, uint64_t D, typename Collision>
+constexpr auto broadphase_collision_distance_squared = [](
+	saw::data<sch::Particle<T,D>>& left,
+	const saw::data<Collision>& coll_l,
+	saw::data<sch::Particle<T,D>>& right,
+	const saw::data<Collision>& coll_r
+) -> std::pair<bool,saw::data<sch::Scalar<T>>>{
+
+	auto rad_l = coll_l.template get<"radius">();
+	auto rad_r = coll_r.template get<"radius">();
+
+	auto& rb_l = left.template get<"rigid_body">();
+	auto& rb_r = right.template get<"rigid_body">();
+
+	auto& pos_l = rb_l.template get<"position">();
+	auto& pos_r = rb_r.template get<"position">();
+
+	auto pos_dist = pos_l - pos_r;
+
+	auto norm_2 = saw::math::dot(pos_dist,pos_dist);
+
+	saw::data<sch::Scalar<T>> two;
+	two.at({}) = 2.0;
+	auto rad_ab_2 = rad_l * rad_l + rad_r * rad_r + rad_r * rad_l * two;
+
+	return std::make_pair((norm_2.at({}).get() < rad_ab_2.at({}).get()), norm_2);
+};
+/**
+*
+*
+*/
+template<typename T, uint64_t D, typename Collision>
+constexpr auto broadphase_collision_check = [](
+	saw::data<sch::Particle<T,D>>& left,
+	const saw::data<Collision>& coll_l,
+	saw::data<sch::Particle<T,D>>& right,
+	const saw::data<Collision>& coll_r
+) -> bool{
+	return broadphase_collision_distance_squared<T,D,Collision>(left,coll_l,right,coll_r).first;
+};
+
+
+
+namespace impl {
+}
+
+}
+}
diff --git a/modules/core/c++/particle/particle_opa.hpp b/modules/core/c++/particle/particle_opa.hpp
new file mode 100644
index 0000000..a5e6c28
--- /dev/null
+++ b/modules/core/c++/particle/particle_opa.hpp
@@ -0,0 +1,49 @@
+#pragma once
+
+#include "../component.hpp"
+#include "common.hpp"
+#include "schema.hpp"
+#include "porosity.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct OneParticleAt {};
+}
+
+template<typename T, typename Descriptor, typename Encode>
+
+class component<T,Descriptor,cmpt::OneParticleAt, Encode> final {
+private:
+	saw::data<sch::Vector<T,Descriptor::D>> pos_;
+	saw::data<sch::Scalar<T>> rad_;
+	saw::data<sch::Scalar<T>> eps_;
+public:
+	component(
+		const saw::data<sch::Vector<T,Descriptor::D>> pos__,
+		const saw::data<sch::Scalar<T>> rad__,
+		const saw::data<sch::Scalar<T>> eps__
+	):
+		pos_{pos__},
+		rad_{rad__},
+		eps_{eps__}
+	{}
+
+	template<typename MacroFieldSchema>
+	void apply(const saw::data<MacroFieldSchema, Encode>& macros, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+		using dfi = df_info<T,Descriptor>;
+
+		auto& porous_f = macros.template get<"porosity">();
+
+		auto& porous = porous_f.at(index);
+
+		auto pos_ind = saw::math::vectorize_data(index);
+
+		// Write out value
+		auto diff = pos_ind.template cast_to<T>() - pos_;
+		auto diff_dot = saw::math::dot(diff,diff);
+		porous = particle_porosity<T,Descriptor::D,por::ParticleSpheroid<T>>::calculate(diff, rad_, eps_);
+	}
+};
+}
+}
diff --git a/modules/core/c++/particle/porosity.hpp b/modules/core/c++/particle/porosity.hpp
new file mode 100644
index 0000000..2f39d2e
--- /dev/null
+++ b/modules/core/c++/particle/porosity.hpp
@@ -0,0 +1,129 @@
+#pragma once
+
+#include "common.hpp"
+#include "schema.hpp"
+#include "../math/n_closest.hpp"
+
+namespace kel {
+namespace lbm {
+namespace por {
+template<typename T>
+struct ParticleSpheroid {};
+}
+
+template<typename T, uint64_t D, typename Coll>
+class particle_porosity {
+public:
+
+	static saw::data<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,Coll>>& part_group, uint64_t p_i, const saw::data<sch::Vector<T,D>>& lbm_pos){
+		auto& mask = part_group.template get<"mask">();
+
+		auto& particles = part_group.template get<"particles">();
+		auto& part_i = particles.at({p_i});
+
+		auto& part_i_rb = part_i.template get<"rigid_body">();
+		auto& pirb = part_i_rb.template get<"position">();
+
+		auto dist = lbm_pos - pirb;
+
+		// index 0 is at 
+
+		return {};
+	}
+};
+
+template<typename T, uint64_t D>
+class particle_porosity<T, D, por::ParticleSpheroid<T>> final {
+public:
+	static saw::data<sch::Scalar<T>> calculate(const saw::data<sch::Vector<T,D>>& lbm_rel_dist, saw::data<sch::Scalar<T>> rad, saw::data<sch::Scalar<T>> eps){
+		saw::data<sch::Scalar<T>> por;
+		
+		auto s_dist_2 = saw::math::dot(lbm_rel_dist,lbm_rel_dist);
+		auto s_dist = saw::math::sqrt(s_dist_2);
+
+		saw::data<sch::Scalar<T>> eps_h;
+		eps_h.at({}) = eps.at({}).get() / 2;
+
+		auto rad_low = rad - eps_h;
+		auto rad_high = rad + eps_h;
+
+		if(s_dist.at({}) <= rad_low.at({})){
+			por.at({}).set(0);
+			return por;
+		}
+		if(s_dist.at({}) >= rad_high.at({})){
+			por.at({}).set(1);
+			return por;
+		}
+
+		{
+			typename saw::native_data_type<T>::type inner = (std::numbers::pi / ( 2 * eps.at({}).get() )) * (s_dist - rad_low).at({}).get();
+			auto cos_inner = std::sin(inner);
+			auto cos_i_2 = cos_inner * cos_inner;
+			por.at({}).set(cos_i_2);
+		}
+
+		return por;
+	}
+};
+
+
+template<typename T, uint64_t D>
+class particle_porosity<T, D, coll::Spheroid<T>> final {
+public:
+	static saw::data<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,coll::Spheroid<T> > >& part_group, uint64_t i, const saw::data<sch::Vector<T,D>>& lbm_pos) {
+		saw::data<sch::Scalar<T>> por;
+
+		auto& parts = part_group.template get<"particles">();
+		auto& pi = parts.at({i});
+		auto& pi_rb = pi.template get<"rigid_body">();
+
+		auto& pi_rb_pos = pi_rb.template get<"position">();
+
+		// Basically the queried position minus the center of the particle
+		auto dist = lbm_pos - pi_rb_pos;
+
+		saw::data<sch::Scalar<T>> dist_len = saw::math::dot(dist,dist);
+		dist_len.at({}).set(std::sqrt(dist_len.at({}).get()));
+
+		auto& coll = part_group.template get<"collision">();
+		const saw::data<sch::Scalar<T>>& rad_d = coll.template get<"radius">();
+		const auto& eps = part_group.template get<"epsilon">();
+		// Move this somewhere
+		
+		saw::data<sch::Scalar<T>> eps_h;
+		eps_h.at({}) = eps.at({}).get() / 2;
+
+		saw::data<sch::Scalar<T>> rad_d_eps_p = rad_d + eps_h;
+		saw::data<sch::Scalar<T>> rad_d_eps_n = rad_d - eps_h;
+
+		// saw::data<sch::Scalar<T>> rad_d_eps_p_2 = rad_d_eps_p * rad_d_eps_p;
+		// saw::data<sch::Scalar<T>> rad_d_eps_n_2 = rad_d_eps_n * rad_d_eps_n;
+		// saw::data<sch::Scalar<T>> rad_2;
+		// rad_2 = rad_d * rad_d;
+
+		saw::data<sch::Scalar<T>> inside;
+		if(dist_len.at({}) <= rad_d_eps_n.at({})){
+			inside.at({}).set(0);
+			return inside;
+		}
+		if(dist_len.at({}) > rad_d_eps_p.at({})){
+			inside.at({}).set(1);
+			return inside;
+		}
+
+		/*
+		 * cos^2 ( ||x-X(t)||_2 - (R-eps/2) )
+		 */
+		{
+			typename saw::native_data_type<T>::type inner = (std::numbers::pi / (eps)) * (dist_len - rad_d_eps_n).at({}).get();
+			auto cos_inner = std::cos(inner);
+			auto cos_i_2 = cos_inner * cos_inner;
+			inside.at({}).set(cos_i_2);
+		}
+		return inside;
+	}
+};
+
+}
+}
diff --git a/modules/core/c++/particle/schema.hpp b/modules/core/c++/particle/schema.hpp
new file mode 100644
index 0000000..714a16f
--- /dev/null
+++ b/modules/core/c++/particle/schema.hpp
@@ -0,0 +1,68 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+
+namespace coll {
+template<typename T>
+struct Spheroid {
+	using ValueSchema = T;
+	using Schema = sch::Struct<
+		sch::Member<sch::Scalar<ValueSchema>,"radius">
+	>;
+};
+}
+
+namespace sch {
+using namespace saw::schema;
+
+namespace impl {
+template<typename T,uint64_t D>
+struct rotation_type_helper;
+
+template<typename T>
+struct rotation_type_helper<T,2u> {
+	using Schema = Scalar<T>;
+};
+
+template<typename T>
+struct rotation_type_helper<T,3u> {
+	using Schema = Vector<T,3u>;
+};
+}
+
+template<typename T, uint64_t D>
+using ParticleRigidBody = Struct<
+	Member<Vector<T,D>, "position">,
+	Member<Vector<T,D>, "position_old">,
+	Member<typename impl::rotation_type_helper<T,D>::Schema, "rotation">,
+	Member<typename impl::rotation_type_helper<T,D>::Schema, "rotation_old">,
+
+	Member<Vector<T,D>, "acceleration">,
+	Member<typename impl::rotation_type_helper<T,D>::Schema, "angular_acceleration">
+>;
+
+
+template<typename T, uint64_t D>
+using Particle = Struct<
+	Member<ParticleRigidBody<T,D>, "rigid_body">
+	// Problem is that dynamic data would two layered
+	// Member<Array<Float64,D>, "mask">,
+>;
+
+template<typename T, uint64_t D, typename CollisionType = coll::Spheroid<T>>
+using ParticleGroup = Struct<
+	Member<Array<T,D>, "mask">,
+	Member<FixedArray<typename CollisionType::Schema,1u>, "collision">,
+	Member<FixedArray<Scalar<T>,1u>, "epsilon">,
+	Member<FixedArray<Scalar<T>,1u>, "mask_step">,
+	Member<FixedArray<Scalar<T>,1u>, "density">,
+	Member<FixedArray<Vector<T,D>,1u>, "center_of_mass">,
+	Member<FixedArray<Scalar<T>,1u>, "total_mass">,
+	Member<Array<Particle<T,D>,1u>, "particles">
+>;
+}
+}
+}
diff --git a/modules/core/c++/psm.hpp b/modules/core/c++/psm.hpp
new file mode 100644
index 0000000..6dc146f
--- /dev/null
+++ b/modules/core/c++/psm.hpp
@@ -0,0 +1,83 @@
+#pragma once
+
+#include "macroscopic.hpp"
+#include "component.hpp"
+#include "equilibrium.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct PSM {};
+}
+
+/**
+ * PSM collision operator for LBM
+ */
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::PSM, Encode> {
+private:
+	saw::data<T> relaxation_;
+	saw::data<T> frequency_;
+public:
+	component(
+		typename saw::native_data_type<T>::type relaxation__
+	):
+		relaxation_{relaxation__}
+	{
+		saw::data<T> one;
+		one = 1.0;
+		frequency_ = one / relaxation_;
+	}
+
+	template<typename CellFieldSchema, typename MacroFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+	
+		using dfi = df_info<T,Descriptor>;
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+		auto& porous_f = macros.template get<"porosity">();
+
+		auto& rho_f = macros.template get<"density">();
+		auto& vel_f = macros.template get<"velocity">();
+		
+		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
+		saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);
+
+		compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel);
+		
+		auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+		saw::data<T> one{1.0};
+		auto& porous = porous_f.at(index);
+		auto flip_porous = one - porous.at({});
+
+		auto& dfs = dfs_old_f.at(index);
+
+		auto dfs_cpy = dfs;
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			uint64_t i_opp = dfi::opposite_index[i];
+			dfs.at({i}) = dfs_cpy.at({i}) + frequency_ * (eq.at(i) - dfs_cpy.at({i})) * porous.at({}) + (dfs_cpy.at({i_opp}) - dfs_cpy.at({i}) ) * flip_porous;
+		}
+
+		auto& force_f = macros.template get<"force">();
+		auto& force = force_f.at(index);
+		for(uint64_t k{0u}; k < Descriptor::D; ++k){
+			force.at({{k}}).set(0);
+		}
+
+		for(uint64_t i{0u}; i < Descriptor::Q; ++i){
+			uint64_t i_opp = dfi::opposite_index[i];
+			auto dfs_diff = dfs.at({i}) - dfs.at({i_opp});
+			for(uint64_t k{0u}; k < Descriptor::D; ++k){
+				force.at({{k}}) = force.at({{k}}) + dfs_diff;
+			}
+		}
+
+		force = force * porous;
+	}
+};
+
+}
+}
diff --git a/modules/core/c++/rar.hpp b/modules/core/c++/rar.hpp
new file mode 100644
index 0000000..dc2b61c
--- /dev/null
+++ b/modules/core/c++/rar.hpp
@@ -0,0 +1,25 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+saw::error_or<void> run_and_record(int argc, char** argv){
+	using RarGit = Struct<
+		Member<String, "commit">,
+		Member<String, "origin">
+	>;
+	using RarCommand = Struct<
+		Member<String, "command">,
+		Member<Array<String>, "arguments">
+	>;
+	using Rar = Struct<
+		Member<Array<RarCommand>, "commands">,
+		Member<UInt64, "runtime_ns">,
+		Member<RarGit, "git">
+	>;
+
+	return saw::make_void();
+}
+}
+}
diff --git a/modules/core/c++/schema.hpp b/modules/core/c++/schema.hpp
new file mode 100644
index 0000000..7712f99
--- /dev/null
+++ b/modules/core/c++/schema.hpp
@@ -0,0 +1,11 @@
+#pragma once
+
+#include <forstio/codec/schema.hpp>
+
+namespace kel {
+namespace lbm {
+namespace sch {
+using namespace saw::schema;
+}
+}
+}
diff --git a/modules/core/c++/simulation/common.hpp b/modules/core/c++/simulation/common.hpp
new file mode 100644
index 0000000..f675a99
--- /dev/null
+++ b/modules/core/c++/simulation/common.hpp
@@ -0,0 +1,7 @@
+#pragma once
+
+namespace kel {
+namespace lbm {
+
+}
+}
diff --git a/modules/core/c++/simulation/hlbm.hpp b/modules/core/c++/simulation/hlbm.hpp
new file mode 100644
index 0000000..93df9a7
--- /dev/null
+++ b/modules/core/c++/simulation/hlbm.hpp
@@ -0,0 +1,12 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+
+class simulation {};
+
+
+}
+}
diff --git a/modules/core/c++/statistics.hpp b/modules/core/c++/statistics.hpp
new file mode 100644
index 0000000..c07ccb7
--- /dev/null
+++ b/modules/core/c++/statistics.hpp
@@ -0,0 +1,11 @@
+#pragma once
+
+namespace kel {
+namespace lbm {
+template<typename T>
+class statistics {
+private:
+public:
+};
+}
+}
diff --git a/modules/core/c++/stream.hpp b/modules/core/c++/stream.hpp
new file mode 100644
index 0000000..dc7cfb3
--- /dev/null
+++ b/modules/core/c++/stream.hpp
@@ -0,0 +1,61 @@
+#pragma once
+
+#include "component.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct Stream {};
+}
+
+template<typename T, typename Descriptor, typename Encode>
+class component<T,Descriptor, cmpt::Stream, Encode> final {
+private:
+public:
+	static constexpr saw::string_literal name = "streaming";
+	static constexpr saw::string_literal after = "collide";
+	static constexpr saw::string_literal before = "";
+
+	template<typename CellFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+		using dfi = df_info<T,Descriptor>;
+		auto g_index = index;
+
+		for(uint64_t i = 0u; i < Descriptor::D; ++i){
+			++g_index.at({{i}});
+		}
+
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+		auto& dfs_new_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">();
+		auto info_f = field.template get<"info">();
+
+		auto info_meta = info_f.get_dims();
+
+		/*
+		bool border = false;
+		for(uint64_t i = 0u; i < Descriptor::D; ++i){
+			auto ind_i = index.at({i});
+			border |= (ind_i.get()) == 0u or (ind_i == info_meta.at({i}));
+		}
+
+		if (not border){
+		}
+		*/
+
+		auto& dfs_new = dfs_new_f.ghost_at(g_index);
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			auto dir = dfi::directions[dfi::opposite_index[i]];
+			auto g_index_nb = g_index;
+			for(uint64_t z = 0u; z < Descriptor::D; ++z){
+				g_index_nb.at({z}).set(g_index.at({z}).get() + dir[z]);
+			}
+
+			dfs_new.at({i}) = dfs_old_f.ghost_at(g_index_nb).at({i});
+		}
+	}
+};
+}
+}
diff --git a/modules/core/c++/util.hpp b/modules/core/c++/util.hpp
new file mode 100644
index 0000000..0bdebd1
--- /dev/null
+++ b/modules/core/c++/util.hpp
@@ -0,0 +1,93 @@
+#pragma once
+
+#include <forstio/string_literal.hpp>
+#include <forstio/codec/data.hpp>
+
+#include <iostream>
+#include <iomanip>
+#include <sys/ioctl.h>
+
+namespace kel {
+namespace lbm {
+/*
+template<typename T, typename Descriptor, typename CellFieldSchema>
+struct is_neighbour {
+	static bool operator()(saw::data<CellFieldSchema>& latt,  saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>& index) {
+		using dfi = df_info<T,Descriptor>;
+
+		if(index.at({0u}).get() == 0u or index.at({1u}).get() == 0u){
+			return false;
+		}
+
+		for(saw::data<sch::UInt64> k{0u}; k.get() < Descriptor::Q; ++k){
+			// TODO
+			saw::data<sch::FixedArray<sch::UInt64,2u>>
+		}
+
+		auto& cell = latt(index);
+	}
+};
+*/
+
+/* Might be stupidly complex
+class cell_info_registry final {
+public:
+	static uint8_t next_reg_id = 0u;
+
+	template<string_literal T>
+	static uint8_t get_registration_id() {
+		static uint8_t reg_id = std::numeric_limit<uint8_t>::max();
+
+		if(reg_id == std::numeric_limit<uint8_t>::max()){
+			reg_id = next_reg_id;
+			++next_reg_id;
+		}
+
+		return reg_id;
+	}
+};
+*/
+
+void print_progress_bar(const uint32_t progress, const uint32_t progress_target){
+	std::cout<<"\r";
+	//	<<"Progress: "
+	//	<<((100 * progress) / progress_target)
+	//	<<"% [";
+
+	const uint32_t progress_min = std::min(progress,progress_target);
+	constexpr uint64_t max_perc_progress = 100u;
+	uint64_t perc_prog = (max_perc_progress*progress_min) / progress_target;
+
+	std::string_view prog_str{"Progress: "};
+	// Progress string + (100 width and perc char) + ([] brackets) + space + pointer
+	uint64_t i{prog_str.size() + 4u + 2u + 1u + 1u};
+
+	std::cout<<prog_str;
+	std::cout<<std::setw(3u)<<perc_prog<<"%";
+	std::cout<<" ";
+	std::cout<<"[";
+
+	uint64_t max_display = []() -> uint64_t{
+		struct winsize w;
+		if(ioctl(STDOUT_FILENO, TIOCGWINSZ,&w) == -1){
+			// Standardized Terminal size
+			return 80u;
+		}
+		return w.ws_col;
+	}();
+	max_display = std::max(max_display,i) - i;
+	uint64_t progress_display = (max_display * progress_min) / progress_target;
+
+	for(i = 0u; i < progress_display; ++i){
+		std::cout<<"#";
+	}
+	for(; i < max_display; ++i){
+		std::cout<<"-";
+	}
+
+	std::cout<<"]";
+
+	std::cout<<std::flush;
+}
+}
+}
diff --git a/modules/core/c++/write_csv.hpp b/modules/core/c++/write_csv.hpp
new file mode 100644
index 0000000..a60e208
--- /dev/null
+++ b/modules/core/c++/write_csv.hpp
@@ -0,0 +1,184 @@
+#pragma once
+
+#include <forstio/error.hpp>
+
+#include <forstio/codec/data.hpp>
+#include <forstio/codec/data_math.hpp>
+
+#include "descriptor.hpp"
+#include "flatten.hpp"
+#include "chunk.hpp"
+
+#include <fstream>
+#include <filesystem>
+
+namespace kel {
+namespace lbm {
+namespace impl {
+
+template<typename CellFieldSchema>
+struct lbm_csv_writer {
+};
+
+template<typename T, uint64_t D>
+struct lbm_csv_writer<sch::Primitive<T,D>> {
+	static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::Primitive<T,D>>& field){
+		if constexpr (std::is_same_v<T,sch::UnsignedInteger> and D == 1u) {
+			csv_file<<field.template cast_to<sch::UInt16>().get();
+		}else{
+			csv_file<<field.get();
+		}
+		return saw::make_void();
+	}
+};
+
+template<typename T, uint64_t D>
+struct lbm_csv_writer<sch::FixedArray<T,D>> {
+	static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::FixedArray<T,D>>& field){
+		saw::data<sch::FixedArray<sch::UInt64,D>> index;
+		for(saw::data<sch::UInt64> it{0}; it.get() < D; ++it){
+			index.at({0u}).set(0u);
+		}
+
+		// csv_file<<"VECTORS "<<name<<" float\n";
+		for(uint64_t i = 0u; i < D; ++i){
+			if(i > 0){
+				csv_file<<",";
+			}
+			csv_file<<field.at({i}).get();
+		}
+		return saw::make_void();
+	}
+};
+
+template<typename T, uint64_t Ghost, uint64_t... D>
+struct lbm_csv_writer<sch::Chunk<T,Ghost,D...>> {
+
+	template<uint64_t d>
+	static saw::error_or<void> apply_d(std::ostream& csv_file, const saw::data<sch::Chunk<T,Ghost,D...>>& field, saw::data<sch::FixedArray<sch::UInt64,sizeof...(D)>>& index){
+		// VTK wants to iterate over z,y,x instead of x,y,z
+		// So we do the same with CSV to stay consistent for now
+		// We could reorder the dimensions, but eh
+		if constexpr ( d > 0u){
+			for(index.at({d-1u}) = 0u; index.at({d-1u}) < field.get_dims().at({d-1u}); ++index.at({d-1u})){
+				auto eov = apply_d<d-1u>(csv_file, field, index);
+			}
+		}else{
+			auto eov = lbm_csv_writer<T>::apply(csv_file, field.at(index));
+			csv_file<<"\n";
+			if(eov.is_error()) return eov;
+		}
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::Chunk<T,Ghost,D...>>& field, std::string_view name){
+		saw::data<sch::FixedArray<sch::UInt64,sizeof...(D)>> index;
+		for(saw::data<sch::UInt64> it{0}; it.get() < sizeof...(D); ++it){
+			index.at({0u}).set(0u);
+		}
+
+		{
+			auto eov = apply_d<sizeof...(D)>(csv_file, field, index);
+			if(eov.is_error()){
+				return eov;
+			}
+		}
+
+		return saw::make_void();
+	}
+};
+
+template<typename T, uint64_t D>
+struct lbm_csv_writer<sch::Vector<T,D>> {
+	static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::Vector<T,D>>& field){
+		static_assert(D > 0, "Non-dimensionality is bad for velocity.");
+
+		// csv_file<<"VECTORS "<<name<<" float\n";
+		for(uint64_t i = 0u; i < D; ++i){
+			if(i > 0){
+				csv_file<<",";
+			}
+			{
+				auto eov = lbm_csv_writer<T>::apply(csv_file,field.at({{i}}));
+				if(eov.is_error()) return eov;
+			}
+		}
+		return saw::make_void();
+	}
+};
+
+template<typename T>
+struct lbm_csv_writer<sch::Scalar<T>> {
+	static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::Scalar<T>>& field){
+		return lbm_csv_writer<T>::apply(csv_file,field.at({}));
+	}
+};
+
+template<typename... MemberT, saw::string_literal... Keys, uint64_t... Ghost, uint64_t... Dims>
+struct lbm_csv_writer<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>> final  {
+	template<uint64_t i>
+	static saw::error_or<void> iterate_i(
+			const std::filesystem::path& csv_dir, const std::string_view& file_base_name, uint64_t d_t, 
+			const saw::data<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>>& field){
+
+		if constexpr ( i < sizeof...(MemberT) ) {
+			using MT = typename saw::parameter_pack_type<i,sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>::type;
+			{
+				std::stringstream sstr;
+				sstr
+					<<file_base_name
+					<<"_"
+					<<MT::Key.view()
+					<<"_"
+					<<d_t
+					<<".csv"
+				;
+				std::ofstream csv_file{csv_dir / sstr.str() };
+
+				if( not csv_file.is_open() ){
+					return saw::make_error<saw::err::critical>("Could not open file.");
+				}
+				// 
+				auto eov = lbm_csv_writer<typename MT::ValueType>::apply(csv_file,field.template get<MT::KeyLiteral>(), MT::KeyLiteral.view());
+				if(eov.is_error()){
+					return eov;
+				}
+			}
+
+			return iterate_i<i+1u>(csv_dir, file_base_name, d_t,field);
+		}
+
+		return saw::make_void();
+	}
+
+
+	static saw::error_or<void> apply(
+			const std::filesystem::path& csv_dir, const std::string_view& file_base_name, uint64_t d_t, 
+			const saw::data<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>>& field){
+		
+		auto& field_0 = field.template get<saw::parameter_key_pack_type<0u,Keys...>::literal>();	
+		auto meta = field_0.get_dims();
+		
+		return iterate_i<0u>(csv_dir,file_base_name, d_t, field);
+	}
+};
+
+}
+
+template<typename Sch>
+saw::error_or<void> write_csv_file(const std::filesystem::path& out_dir, const std::string_view& file_name, uint64_t d_t, const saw::data<Sch>& field){
+
+	auto csv_dir = out_dir / "csv";
+	{
+		std::error_code ec;
+		std::filesystem::create_directories(csv_dir,ec);
+		if(ec != std::errc{}){
+			return saw::make_error<saw::err::critical>("Could not create directory for write_csv_file function");
+		}
+	}
+	auto eov = impl::lbm_csv_writer<Sch>::apply(csv_dir, file_name, d_t, field);
+	return eov;
+}
+
+}
+}
diff --git a/modules/core/c++/write_json.hpp b/modules/core/c++/write_json.hpp
new file mode 100644
index 0000000..19bbfa6
--- /dev/null
+++ b/modules/core/c++/write_json.hpp
@@ -0,0 +1,27 @@
+#pragma once
+
+#include <forstio/error.hpp>
+
+#include <forstio/codec/data.hpp>
+#include <forstio/codec/data_math.hpp>
+
+#include <forstio/codec/json/json.hpp>
+
+
+namespace kel {
+namespace lbm {
+
+template<typename Schema>
+saw::error_or<void> write_json(const std::filesystem::path& file_name, const saw::data<Schema>& field){
+	saw::data<Schema, saw::encode::Json> j_data;
+	{
+		saw::codec<Schema, saw::encode::Json> j_codec;
+		auto eov = j_codec.encode(field, j_data);
+		if(eov.is_error()){
+			return eov;
+		}
+	}
+	return saw::make_void();
+}
+}
+}
diff --git a/modules/core/c++/write_vtk.hpp b/modules/core/c++/write_vtk.hpp
new file mode 100644
index 0000000..e852172
--- /dev/null
+++ b/modules/core/c++/write_vtk.hpp
@@ -0,0 +1,256 @@
+#pragma once
+
+#include <forstio/error.hpp>
+
+#include <forstio/codec/data.hpp>
+#include <forstio/codec/data_math.hpp>
+
+#include "descriptor.hpp"
+#include "flatten.hpp"
+#include "chunk.hpp"
+
+#include <fstream>
+#include <filesystem>
+
+namespace kel {
+namespace lbm {
+namespace impl {
+
+template<typename CellFieldSchema>
+struct lbm_vtk_writer {
+};
+
+template<typename T, uint64_t D>
+struct lbm_vtk_writer<sch::Primitive<T,D>> {
+	static saw::error_or<void> apply_header(std::ostream& vtk_file, std::string_view name){
+		vtk_file<<"SCALARS "<<name<<" float 1\n";
+		vtk_file<<"LOOKUP_TABLE default\n";
+		return saw::make_void();
+	}
+	static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::Primitive<T,D>>& field){
+		if constexpr (std::is_same_v<T,sch::UnsignedInteger> and D == 1u) {
+			vtk_file<<field.template cast_to<sch::UInt16>().get()<<"\n";
+		}else{
+			vtk_file<<field.get()<<"\n";
+		}
+		return saw::make_void();
+	}
+};
+
+template<typename T, uint64_t D>
+struct lbm_vtk_writer<sch::FixedArray<T,D>> {
+	static saw::error_or<void> apply_header(std::ostream& vtk_file, std::string_view name){
+		vtk_file<<"VECTORS "<<name<<" float\n";
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::FixedArray<T,D>>& field){
+		saw::data<sch::FixedArray<sch::UInt64,D>> index;
+		for(saw::data<sch::UInt64> it{0}; it.get() < D; ++it){
+			index.at({0u}).set(0u);
+		}
+
+		// vtk_file<<"VECTORS "<<name<<" float\n";
+		for(uint64_t i = 0u; i < D; ++i){
+			if(i > 0){
+				vtk_file<<" ";
+			}
+			vtk_file<<field.at({i}).get();
+		}
+		for(uint64_t i = D; i < 3; ++i){
+			vtk_file<<" 0";
+		}
+		vtk_file<<"\n";
+		return saw::make_void();
+	}
+};
+
+template<typename T, uint64_t Ghost, uint64_t... D>
+struct lbm_vtk_writer<sch::Chunk<T,Ghost,D...>> {
+
+	template<uint64_t d>
+	static saw::error_or<void> apply_d(std::ostream& vtk_file, const saw::data<sch::Chunk<T,Ghost,D...>>& field, saw::data<sch::FixedArray<sch::UInt64,sizeof...(D)>>& index){
+		// VTK wants to iterate over z,y,x instead of x,y,z
+		// We could reorder the dimensions, but eh
+		if constexpr ( d > 0u){
+			for(index.at({d-1u}) = 0u; index.at({d-1u}) < field.get_dims().at({d-1u}); ++index.at({d-1u})){
+				auto eov = apply_d<d-1u>(vtk_file, field, index);
+			}
+		}else{
+			auto eov = lbm_vtk_writer<T>::apply(vtk_file, field.at(index));
+			if(eov.is_error()) return eov;
+		}
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::Chunk<T,Ghost,D...>>& field, std::string_view name){
+		{
+			auto eov = lbm_vtk_writer<T>::apply_header(vtk_file,name);
+			if(eov.is_error()){
+				return eov;
+			}
+		}
+		saw::data<sch::FixedArray<sch::UInt64,sizeof...(D)>> index;
+		for(saw::data<sch::UInt64> it{0}; it.get() < sizeof...(D); ++it){
+			index.at({0u}).set(0u);
+		}
+
+		{
+			auto eov = apply_d<sizeof...(D)>(vtk_file, field, index);
+			if(eov.is_error()){
+				return eov;
+			}
+		}
+
+		vtk_file<<"\n";
+
+		return saw::make_void();
+	}
+};
+
+template<typename T, uint64_t D>
+struct lbm_vtk_writer<sch::Vector<T,D>> {
+	static saw::error_or<void> apply_header(std::ostream& vtk_file, std::string_view name){
+		vtk_file<<"VECTORS "<<name<<" float\n";
+		return saw::make_void();
+	}
+	static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::Vector<T,D>>& field){
+		static_assert(D > 0, "Non-dimensionality is bad for velocity.");
+		static_assert(D <= 3, "4th dimension as well. Mostly due to vtk.");
+
+		// vtk_file<<"VECTORS "<<name<<" float\n";
+		for(uint64_t i = 0u; i < D; ++i){
+			if(i > 0){
+				vtk_file<<" ";
+			}
+			vtk_file<<field.at({{i}}).get();
+		}
+		for(uint64_t i = D; i < 3; ++i){
+			vtk_file<<" 0";
+		}
+		vtk_file<<"\n";
+		return saw::make_void();
+	}
+};
+
+template<typename T>
+struct lbm_vtk_writer<sch::Scalar<T>> {
+	static saw::error_or<void> apply_header(std::ostream& vtk_file, std::string_view name){
+		vtk_file<<"SCALARS "<<name<<" float 1\n";
+		vtk_file<<"LOOKUP_TABLE default\n";
+		return saw::make_void();
+	}
+	static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::Scalar<T>>& field){
+		// vtk_file<<"VECTORS "<<name<<" float\n";
+		vtk_file<<field.at({}).get();
+		vtk_file<<"\n";
+		return saw::make_void();
+	}
+};
+
+template<typename... MemberT, saw::string_literal... Keys, uint64_t... Ghost, uint64_t... Dims>
+struct lbm_vtk_writer<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>> final  {
+	template<uint64_t i>
+	static saw::error_or<void> iterate_i(std::ostream& vtk_file,
+		const saw::data<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>>& field){
+
+		if constexpr ( i < sizeof...(MemberT) ) {
+			using MT = typename saw::parameter_pack_type<i,sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>::type;
+			{
+				auto eov = lbm_vtk_writer<typename MT::ValueType>::apply(vtk_file,field.template get<MT::KeyLiteral>(), MT::KeyLiteral.view());
+				if(eov.is_error()){
+					return eov;
+				}
+			}
+
+			return iterate_i<i+1u>(vtk_file, field);
+		}
+
+		return saw::make_void();
+	}
+
+
+	static saw::error_or<void> apply(std::ostream& vtk_file,
+		const saw::data<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>>& field){
+		
+		vtk_file
+			<<"# vtk DataFile Version 3.0\n"
+			<<"LBM File\n"
+			<<"ASCII\n"
+			<<"DATASET STRUCTURED_POINTS\n"
+			<<"SPACING 1.0 1.0 1.0\n"
+			<<"ORIGIN 0.0 0.0 0.0\n"
+		;
+
+		auto& field_0 = field.template get<saw::parameter_key_pack_type<0u,Keys...>::literal>();	
+		auto meta = field_0.get_dims();
+		
+		saw::data<sch::UInt64> pd_size{1u};
+		// DIMENSIONS
+
+		{
+			vtk_file << "DIMENSIONS";
+			pd_size.set(saw::ct_multiply<uint64_t,Dims...>::value);
+			
+			static_assert(saw::ct_multiply<uint64_t,Dims...>::value > 0u, "Invalid Dim size resulting in length 0u");
+
+			for(saw::data<sch::UInt64> i{0u}; i.get() < sizeof...(Dims); ++i){
+				vtk_file << " " << meta.at(i).get(); 
+			}
+			for(saw::data<sch::UInt64> i{sizeof...(Dims)}; i.get() < 3u; ++i){
+				vtk_file << " 1";
+			}
+
+			vtk_file << "\n";
+		}
+		if constexpr (sizeof...(MemberT) > 0u){
+			// POINT DATA
+			{
+				vtk_file << "\nPOINT_DATA " << pd_size.get() <<"\n";
+			}
+
+			// HEADER TO BODY
+			{
+				vtk_file << "\n";
+			}
+		}
+		
+		return iterate_i<0u>(vtk_file, field);
+	}
+};
+
+}
+
+template<typename Sch>
+saw::error_or<void> write_vtk_file(const std::filesystem::path& out_dir, const std::string_view& file_name, uint64_t d_t, const saw::data<Sch>& field){
+
+	auto vtk_dir = out_dir / "vtk";
+	{
+		std::error_code ec;
+		std::filesystem::create_directories(vtk_dir,ec);
+		if(ec != std::errc{}){
+			return saw::make_error<saw::err::critical>("Could not create directory for write_vtk_file function");
+		}
+	}
+	std::string ft_name{file_name};
+
+	std::stringstream sstr;
+	sstr
+		<<file_name
+		<<"_"
+		<<d_t
+		<<".vtk"
+	;
+
+	std::ofstream vtk_file{vtk_dir / sstr.str() };
+
+	if( not vtk_file.is_open() ){
+		return saw::make_error<saw::err::critical>("Could not open file.");
+	}
+
+	auto eov = impl::lbm_vtk_writer<Sch>::apply(vtk_file, field);
+	return eov;
+}
+
+}
+}
diff --git a/modules/core/tests/SConscript b/modules/core/tests/SConscript
new file mode 100644
index 0000000..d1b381e
--- /dev/null
+++ b/modules/core/tests/SConscript
@@ -0,0 +1,32 @@
+#!/bin/false
+
+import os
+import os.path
+import glob
+
+
+Import('env')
+
+dir_path = Dir('.').abspath
+
+# Environment for base library
+test_cases_env = env.Clone();
+
+test_cases_env.Append(LIBS=['forstio-test']);
+
+test_cases_env.sources = sorted(glob.glob(dir_path + "/*.cpp"))
+test_cases_env.headers = sorted(glob.glob(dir_path + "/*.hpp"))
+
+env.sources += test_cases_env.sources;
+env.headers += test_cases_env.headers;
+
+objects_static = []
+test_cases_env.add_source_files(objects_static, test_cases_env.sources, shared=False);
+test_cases_env.program = test_cases_env.Program('#bin/tests', [objects_static]);
+# , env.library_static]);
+
+# Set Alias
+env.Alias('test', test_cases_env.program);
+env.Alias('check', test_cases_env.program);
+
+env.targets += ['test','check'];
diff --git a/modules/core/tests/chunk.cpp b/modules/core/tests/chunk.cpp
new file mode 100644
index 0000000..008d4c6
--- /dev/null
+++ b/modules/core/tests/chunk.cpp
@@ -0,0 +1,47 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/chunk.hpp"
+
+namespace {
+namespace sch {
+using namespace saw::schema;
+}
+
+SAW_TEST("Chunk Ghost size 0"){
+	using namespace kel;
+
+	using TestChunk = lbm::sch::Chunk<sch::UInt32, 0u, 2u, 3u, 4u>;
+	using TestArray = sch::FixedArray<sch::UInt32, 2u, 3u, 4u>;
+
+	SAW_EXPECT((std::is_same_v<TestChunk::InnerSchema,TestArray>), "Types are not identical");
+}
+
+SAW_TEST("Chunk Ghost size 1"){
+	using namespace kel;
+
+	using TestChunk = lbm::sch::Chunk<sch::UInt32,1u,2u,3u,4u>;
+	using TestArray = sch::FixedArray<sch::UInt32,4u,5u,6u>;
+
+	SAW_EXPECT((std::is_same_v<TestChunk::InnerSchema,TestArray>), "Types are not identical");
+}
+
+SAW_TEST("Chunk Ghost size 3"){
+	using namespace kel;
+
+	using TestChunk = lbm::sch::Chunk<sch::UInt32,3u,2u,3u,4u,10u,6u>;
+	using TestArray = sch::FixedArray<sch::UInt32,8u,9u,10u,16u,12u>;
+
+	SAW_EXPECT((std::is_same_v<TestChunk::InnerSchema,TestArray>), "Types are not identical");
+}
+
+SAW_TEST("Chunk 2D setup"){
+	using namespace kel;
+
+	using TestChunk = lbm::sch::Chunk<sch::Float32, 1u, 2u, 2u>;
+
+	auto chunk = saw::heap<saw::data<TestChunk>>();
+	chunk->at({{0u,0u}}).set(1.0f);
+	SAW_EXPECT((chunk->ghost_at({{1u,1u}}).get() == 1.0f),"Check if ghost is at shifted spot.");
+}
+
+}
diff --git a/modules/core/tests/collision.cpp b/modules/core/tests/collision.cpp
new file mode 100644
index 0000000..cd53336
--- /dev/null
+++ b/modules/core/tests/collision.cpp
@@ -0,0 +1,6 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/collision.hpp"
+
+namespace {
+}
diff --git a/modules/core/tests/converter.cpp b/modules/core/tests/converter.cpp
new file mode 100644
index 0000000..4fc536f
--- /dev/null
+++ b/modules/core/tests/converter.cpp
@@ -0,0 +1,26 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/converter.hpp"
+
+
+namespace {
+namespace sch {
+using namespace kel::lbm::sch;
+
+using T = Float64;
+}
+
+SAW_TEST("Si Meter to Lbm Meter"){
+	using namespace kel;
+
+	lbm::converter<sch::T> converter{
+		{0.1},
+		{0.1}
+	};
+
+	saw::data<sch::SiMeter<sch::T>> si_m{1.0};
+
+	auto lbm_m = converter.meter_si_to_lbm(si_m);
+	SAW_EXPECT(lbm_m.handle().get() == 10.0, "Correct si to lbm conversion");
+}
+}
diff --git a/modules/core/tests/descriptor.cpp b/modules/core/tests/descriptor.cpp
new file mode 100644
index 0000000..7f743ce
--- /dev/null
+++ b/modules/core/tests/descriptor.cpp
@@ -0,0 +1,39 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/descriptor.hpp"
+
+
+namespace {
+template<typename Descriptor>
+void check_opposite_dirs(){
+	using namespace kel;
+
+	using dfi = lbm::df_info<lbm::sch::Float32, Descriptor>;
+
+	for(uint64_t k = 0u; k < Descriptor::Q; ++k){
+		auto k_inv = dfi::opposite_index[k];
+
+		for(uint64_t i = 0u; i < Descriptor::D; ++i){
+			SAW_EXPECT(dfi::directions[k][i] == (-1*dfi::directions[k_inv][i]), "Opposites are inconsistent");
+		}
+	}
+}
+
+SAW_TEST("Opposites and Dirs D1Q3"){
+	using namespace kel;
+	check_opposite_dirs<lbm::sch::Descriptor<1,3>>();
+}
+
+SAW_TEST("Opposites and Dirs D2Q5"){
+	using namespace kel;
+	check_opposite_dirs<lbm::sch::Descriptor<2,5>>();
+}
+SAW_TEST("Opposites and Dirs D2Q9"){
+	using namespace kel;
+	check_opposite_dirs<lbm::sch::Descriptor<2,9>>();
+}
+SAW_TEST("Opposites and Dirs D3Q27"){
+	using namespace kel;
+	check_opposite_dirs<lbm::sch::Descriptor<3,27>>();
+}
+}
diff --git a/modules/core/tests/equilibrium.cpp b/modules/core/tests/equilibrium.cpp
new file mode 100644
index 0000000..20a1f08
--- /dev/null
+++ b/modules/core/tests/equilibrium.cpp
@@ -0,0 +1,41 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/equilibrium.hpp"
+
+
+namespace {
+
+template<typename Descriptor>
+void check_equilibrium(){
+	using namespace kel;
+
+	using dfi = lbm::df_info<lbm::sch::Float64,Descriptor>;
+
+	saw::data<lbm::sch::Scalar<lbm::sch::Float64>> rho;
+	rho.at({}).set(1.0);
+	saw::data<lbm::sch::Vector<lbm::sch::Float64,Descriptor::D>> vel;
+	for(saw::data<lbm::sch::UInt64> i{0u}; i.get() < Descriptor::D; ++i){
+		vel.at({{i}}) = {0.0};
+	}
+	auto eq = lbm::equilibrium<lbm::sch::Float64,Descriptor>(rho,vel);
+
+	for(saw::data<lbm::sch::UInt64> i{0u}; i.get() < Descriptor::Q; ++i){
+		SAW_EXPECT(eq.at(i).get() == dfi::weights[i.get()], std::string{"No velocity and normalized rho should be exactly the weights: "} + std::to_string(eq.at(i).get()) + std::string{" "} + std::to_string(dfi::weights[i.get()]));
+	}
+}
+
+SAW_TEST("Equilibrium at rest D1Q3"){
+	using namespace kel;
+	check_equilibrium<lbm::sch::Descriptor<1,3>>();
+}
+
+SAW_TEST("Equilibrium at rest D2Q5"){
+	using namespace kel;
+	check_equilibrium<lbm::sch::Descriptor<2,5>>();
+}
+
+SAW_TEST("Equilibrium at rest D2Q9"){
+	using namespace kel;
+	check_equilibrium<lbm::sch::Descriptor<2,9>>();
+}
+}
diff --git a/modules/core/tests/flatten.cpp b/modules/core/tests/flatten.cpp
new file mode 100644
index 0000000..84b3fa4
--- /dev/null
+++ b/modules/core/tests/flatten.cpp
@@ -0,0 +1,22 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/flatten.hpp"
+
+namespace {
+namespace sch {
+using namespace saw::schema;
+}
+
+SAW_TEST("Flatten Index Stride"){
+	using namespace kel;
+
+	constexpr saw::data<sch::UInt64> zero = lbm::flatten_index<sch::UInt64,3u>::stride<0u>({{2u,4u,3u}});
+	constexpr saw::data<sch::UInt64> one  = lbm::flatten_index<sch::UInt64,3u>::stride<1u>({{2u,4u,3u}});
+	constexpr saw::data<sch::UInt64> two  = lbm::flatten_index<sch::UInt64,3u>::stride<2u>({{2u,4u,3u}});
+
+	SAW_EXPECT(zero.get() == 1u, "Zero is correct");
+	SAW_EXPECT(one.get() == 2u, "Zero is correct");
+	SAW_EXPECT(two.get() == 8u, "Zero is correct");
+}
+
+}
diff --git a/modules/core/tests/iterator.cpp b/modules/core/tests/iterator.cpp
new file mode 100644
index 0000000..919c12c
--- /dev/null
+++ b/modules/core/tests/iterator.cpp
@@ -0,0 +1,72 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/iterator.hpp"
+
+#include <iostream>
+
+namespace {
+namespace sch {
+using namespace kel::lbm::sch;
+}
+
+SAW_TEST("Old Iterate"){
+	using namespace kel;
+
+	saw::data<sch::FixedArray<sch::UInt64,2u>> start{{0u,0u}};
+	saw::data<sch::FixedArray<sch::UInt64,2u>> end{{3u,3u}};
+
+	lbm::iterate_over([](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+		std::cout<<"Index: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl;
+	}, start, end);
+}
+
+SAW_TEST("CT Tuple Iterator"){
+	using namespace kel;
+
+	saw::data<
+		sch::Tuple<
+			sch::Float32,
+			sch::Int16,
+			sch::Float64,
+			sch::UInt32
+		>
+	> tup;
+
+	lbm::ct_tuple_iterator<std::decay_t<decltype(tup)>::Schema>::apply([](auto& val) -> saw::error_or<void> {
+		return saw::make_void();
+	},tup);
+}
+
+SAW_TEST("Old Iterate with Distance 1"){
+	using namespace kel;
+
+	saw::data<sch::FixedArray<sch::UInt64,2u>> start{{0u,0u}};
+	saw::data<sch::FixedArray<sch::UInt64,2u>> end{{4u,4u}};
+
+	lbm::iterate_over([](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+		std::cout<<"Index: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl;
+	}, start, end, {{1u,1u}});
+}
+
+SAW_TEST("Iterate - 2D"){
+	using namespace kel;
+
+	saw::data<sch::FixedArray<sch::UInt64,2u>> start{{0u,0u}};
+	saw::data<sch::FixedArray<sch::UInt64,2u>> end{{3u,3u}};
+
+	lbm::iterator<2u>::apply([](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+		std::cout<<"Index: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl;
+	}, start, end);
+}
+
+SAW_TEST("Iterate with Distance 1 - 2D"){
+	using namespace kel;
+
+	saw::data<sch::FixedArray<sch::UInt64,2u>> start{{0u,0u}};
+	saw::data<sch::FixedArray<sch::UInt64,2u>> end{{4u,4u}};
+
+	lbm::iterator<2u>::apply([](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+		std::cout<<"Index: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl;
+	}, start, end, {{1u,1u}});
+}
+}
diff --git a/modules/core/tests/math.cpp b/modules/core/tests/math.cpp
new file mode 100644
index 0000000..d456ce8
--- /dev/null
+++ b/modules/core/tests/math.cpp
@@ -0,0 +1,79 @@
+#include <forstio/test/suite.hpp>
+
+#include <iostream>
+#include "../c++/math/n_linear.hpp"
+
+namespace {
+namespace sch {
+using namespace saw::schema;
+}
+
+SAW_TEST("Math 1-Linear"){
+	using namespace kel;
+	
+	saw::data<sch::FixedArray<sch::Vector<sch::Float64,1u>,2u>> field;
+	{
+		field.at({{0u}}).at({{0u}}).set(-1.0);
+		field.at({{1u}}).at({{0u}}).set(2.0);
+	}
+	saw::data<sch::Vector<sch::Float64,1u>> pos;
+	pos.at({{0u}}).set(0.3);
+}
+
+SAW_TEST("Math/Floor Index and Fraction from Position"){
+	using namespace kel;
+
+	saw::data<sch::Vector<sch::Float64,2u>> pos;
+
+	{
+		pos.at({{0u}}) = 43.999;
+		pos.at({{1u}}) = -50.0;
+	}
+
+	auto ind_frac = lbm::position_to_index_and_fraction(pos);
+	auto& ind = ind_frac.template get<0u>();
+	for(uint64_t i = 0u; i < 2u; ++i){
+		std::cout<<ind.at({{i}}).get()<<" ";
+	}
+	std::cout<<std::endl;
+	auto& frac = ind_frac.template get<1u>();
+	for(uint64_t i = 0u; i < 2u; ++i){
+		std::cout<<frac.at({{i}}).get()<<" ";
+	}
+	std::cout<<std::endl;
+
+	SAW_EXPECT(true, "Default true check");
+}
+
+SAW_TEST("Math/Floor Index and Fraction from Position and Bounded"){
+	using namespace kel;
+
+	saw::data<sch::Vector<sch::Float64,2u>> pos;
+
+	{
+		pos.at({{0u}}) = 43.999;
+		pos.at({{1u}}) = -50.0;
+	}
+
+	saw::data<sch::Vector<sch::UInt64,2U>> bound;
+	{
+		bound.at({{0u}}) = 32u;
+		bound.at({{1u}}) = 16u;
+	}
+
+	auto ind_frac = lbm::position_to_index_and_fraction_bounded(pos,bound);
+	auto& ind = ind_frac.template get<0u>();
+	for(uint64_t i = 0u; i < 2u; ++i){
+		std::cout<<ind.at({{i}}).get()<<" ";
+	}
+	std::cout<<std::endl;
+	auto& frac = ind_frac.template get<1u>();
+	for(uint64_t i = 0u; i < 2u; ++i){
+		std::cout<<frac.at({{i}}).get()<<" ";
+	}
+	std::cout<<std::endl;
+
+	SAW_EXPECT(true, "Default true check");
+}
+
+}
diff --git a/modules/core/tests/memory.cpp b/modules/core/tests/memory.cpp
new file mode 100644
index 0000000..d59d212
--- /dev/null
+++ b/modules/core/tests/memory.cpp
@@ -0,0 +1,58 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/memory.hpp"
+
+namespace {
+namespace sch {
+using namespace saw::schema;
+
+using TStruct = Struct<
+	Member<Float32,"a">,
+	Member<Int8, "b">,
+	Member<UInt64, "c">
+>;
+
+using TChunk = kel::lbm::sch::Chunk<Float32,1u,2u,3u,4u>;
+using TBChunk = kel::lbm::sch::Chunk<UInt32,1u,2u,3u,4u>;
+
+using TChunkStruct = Struct<
+	Member<TChunk, "a">,
+	Member<TBChunk, "b">
+>;
+}
+
+SAW_TEST("Memory Estimate"){
+	using namespace kel::lbm;
+
+	SAW_EXPECT((memory_estimate<sch::Float32>().get() == 4u), "Float32 isn't 4 bytes" );
+	SAW_EXPECT((memory_estimate<sch::FixedArray<sch::Float32,5u,3u>>().get() == 60u), "FixedArray<Float32,5u,3u> isn't 60 bytes");
+	SAW_EXPECT((memory_estimate<sch::FixedArray<sch::Float32,5u,3u>, sch::UInt8>().get() == 61u), "FixedArray<Float32,5u,3u> + UInt8 isn't 61 bytes");
+}
+
+SAW_TEST("Memory Estimate Struct"){
+	using namespace kel::lbm;
+	
+	SAW_EXPECT((memory_estimate<sch::TStruct>().get() == 13u), "TStruct isn't 13 bytes" );
+	// SAW_EXPECT((memory_estimate<sch::TStruct>().get() == 13u), "TStruct isn't 13 bytes" );
+}
+
+SAW_TEST("Memory Estimate Scalar"){
+	using namespace kel::lbm;
+	
+	SAW_EXPECT((memory_estimate<sch::Scalar<sch::UInt8>>().get() == 1u), "Scalar of UInt8 isn't 1 bytes" );
+	// SAW_EXPECT((memory_estimate<sch::TStruct>().get() == 13u), "TStruct isn't 13 bytes" );
+}
+
+SAW_TEST("Memory Estimate Chunk"){
+	using namespace kel::lbm;
+	
+	SAW_EXPECT((memory_estimate<sch::TChunk>().get() == 480u), std::string{"TChunk isn't 480 bytes, but is "} + std::to_string(memory_estimate<sch::TChunk>().get()) );
+}
+
+SAW_TEST("Memory Estimate Struct of Chunk"){
+	using namespace kel::lbm;
+	
+	SAW_EXPECT((memory_estimate<sch::TChunkStruct>().get() == 960u), std::string{"TChunkStruct isn't 960 bytes, but is "} + std::to_string(memory_estimate<sch::TChunk>().get()) );
+}
+
+}
diff --git a/modules/core/tests/particles.cpp b/modules/core/tests/particles.cpp
new file mode 100644
index 0000000..133c343
--- /dev/null
+++ b/modules/core/tests/particles.cpp
@@ -0,0 +1,281 @@
+#include <forstio/test/suite.hpp>
+
+#include <iostream>
+#include "../c++/particle/particle.hpp"
+
+namespace {
+namespace sch {
+using namespace kel::lbm::sch;
+
+using T = Float64;
+}
+SAW_TEST("Verlet step 2D - Planar"){
+	using namespace kel;
+
+	saw::data<sch::Particle<sch::T,2u>> particle;
+	auto& body = particle.template get<"rigid_body">();
+	auto& pos = body.template get<"position">();
+	auto& pos_old = body.template get<"position_old">();
+
+	// auto& rot = body.template get<"rotation">();
+	auto& acc = body.template get<"acceleration">();
+
+	acc.at({{0}}).set({1.0});
+
+	saw::data<sch::Scalar<sch::T>> dt;
+	dt.at({}).set(0.5);
+	lbm::verlet_step_lambda<sch::T,2u>(particle,dt);
+
+	SAW_EXPECT(pos.at({{0}}).get() == 0.25, std::string{"Incorrect Pos X: "} + std::to_string(pos.at({{0}}).get()));
+	SAW_EXPECT(pos.at({{1}}).get() == 0.0, std::string{"Incorrect Pos Y: "} + std::to_string(pos.at({{1}}).get()));
+}
+/*
+SAW_TEST("No Collision Spheroid 2D"){
+	using namespace kel;
+
+	saw::data<sch::Particle<sch::T,2u>> part_a;
+	{
+		auto& body = part_a.template get<"rigid_body">();
+		auto& pos = body.template get<"position">();
+
+		pos.at({{0u}}) = 0.1;
+		pos.at({{1u}}) = 0.2;
+
+	}
+	saw::data<sch::Particle<sch::T,2u>> part_b;
+	{
+		auto& body = part_b.template get<"rigid_body">();
+		auto& pos = body.template get<"position">();
+
+		pos.at({{0u}}) = -2.1;
+		pos.at({{1u}}) = 0.0;
+	}
+
+	bool have_collided = lbm::broadphase_collision_check<sch::T,2u>(part_a,part_b);
+	SAW_EXPECT(not have_collided, "Particles shouldn't collide");
+}
+*/
+/*
+SAW_TEST("Collision Spheroid 2D"){
+	using namespace kel;
+
+	saw::data<sch::Particle<sch::T,2u>> part_a;
+	{
+		auto& body = part_a.template get<"rigid_body">();
+		auto& pos = body.template get<"position">();
+		auto& coll = part_a.template get<"collision">();
+		auto& radius = coll.template get<"radius">();
+
+		radius.at({}).set(1.0);
+
+		pos.at({{0u}}) = 0.1;
+		pos.at({{1u}}) = 0.2;
+
+	}
+	saw::data<sch::Particle<sch::T,2u>> part_b;
+	{
+		auto& body = part_b.template get<"rigid_body">();
+		auto& pos = body.template get<"position">();
+		auto& coll = part_b.template get<"collision">();
+		auto& radius = coll.template get<"radius">();
+
+		radius.at({}).set(1.0);
+		
+		pos.at({{0u}}) = -1.5;
+		pos.at({{1u}}) = 0.0;
+	}
+
+	bool have_collided = lbm::broadphase_collision_check<sch::T,2u>(part_a,part_b);
+	SAW_EXPECT(have_collided, "Particles should collide");
+}
+*/
+/*
+SAW_TEST("Moving particles 2D"){
+	using namespace kel;
+
+	saw::data<sch::Particle<sch::T,2u>> part_a;
+	{
+		auto& body = part_a.template get<"rigid_body">();
+		auto& pos = body.template get<"position">();
+		auto& pos_old = body.template get<"position_old">();
+		auto& coll = part_a.template get<"collision">();
+		auto& radius = coll.template get<"radius">();
+		auto& acc = body.template get<"acceleration">();
+
+		radius.at({}).set(1.0);
+
+		pos.at({{0u}}) = -5.0;
+		pos.at({{1u}}) = 0.2;
+		pos_old = pos;
+		acc.at({{0u}}) = 0.1;
+
+	}
+	saw::data<sch::Particle<sch::T,2u>> part_b;
+	{
+		auto& body = part_b.template get<"rigid_body">();
+		auto& pos = body.template get<"position">();
+		auto& pos_old = body.template get<"position_old">();
+		auto& coll = part_b.template get<"collision">();
+		auto& radius = coll.template get<"radius">();
+		auto& acc = body.template get<"acceleration">();
+
+		radius.at({}).set(1.0);
+		
+		pos.at({{0u}}) = 5.0;
+		pos.at({{1u}}) = 0.0;
+		pos_old = pos;
+		acc.at({{0u}}) = -0.1;
+	}
+
+	saw::data<sch::Scalar<sch::T>> dt;
+	dt.at({}).set(0.5);
+	bool has_collided = false;
+	for(uint64_t i = 0u; i < 32u; ++i){
+		lbm::verlet_step_lambda<sch::T,2u>(part_a,dt);
+		lbm::verlet_step_lambda<sch::T,2u>(part_b,dt);
+
+		has_collided = lbm::broadphase_collision_check<sch::T,2u>(part_a,part_b);
+		if(has_collided){
+			std::cout<<"Collided at: "<<i<<std::endl;
+			break;
+		}
+	}
+
+	SAW_EXPECT(has_collided, "Expecting collision within those steps");
+}
+*/
+SAW_TEST("Particle Matrix Rotation"){
+	using namespace kel;
+
+	saw::data<sch::Particle<sch::T,2u>> part_a;
+	{
+		auto& body = part_a.template get<"rigid_body">();
+		auto& pos = body.template get<"position">();
+		auto& pos_old = body.template get<"position_old">();
+		auto& acc = body.template get<"acceleration">();
+
+		pos.at({{0u}}) = -5.0;
+		pos.at({{1u}}) = 0.2;
+		pos_old = pos;
+	}
+}
+/*
+
+SAW_TEST("Particle Collision Impulse"){
+	using namespace kel;
+
+	saw::data<sch::Particle<sch::T,2u>> part_a;
+	{
+		auto& body = part_a.template get<"rigid_body">();
+		auto& pos = body.template get<"position">();
+		auto& pos_old = body.template get<"position_old">();
+		auto& coll = part_a.template get<"collision">();
+		auto& radius = coll.template get<"radius">();
+		auto& acc = body.template get<"acceleration">();
+
+		radius.at({}).set(1.0);
+
+		pos.at({{0u}}) = -5.0;
+		pos.at({{1u}}) = 0.2;
+		pos_old = pos;
+		acc.at({{0u}}) = 0.1;
+
+	}
+	saw::data<sch::Particle<sch::T,2u>> part_b;
+	{
+		auto& body = part_b.template get<"rigid_body">();
+		auto& pos = body.template get<"position">();
+		auto& pos_old = body.template get<"position_old">();
+		auto& coll = part_b.template get<"collision">();
+		auto& radius = coll.template get<"radius">();
+		auto& acc = body.template get<"acceleration">();
+
+		radius.at({}).set(1.5);
+		
+		pos.at({{0u}}) = 5.0;
+		pos.at({{1u}}) = 0.0;
+		pos_old = pos;
+		acc.at({{0u}}) = -0.1;
+	}
+}
+*/
+/*
+
+SAW_TEST("Minor Test for mask"){
+	using namespace kel;
+
+	lbm::particle_circle_geometry<sch::T> geo;
+
+	auto mask = geo.generate_mask<sch::T>(9u,1u);
+
+	auto& grid = mask.template get<"grid">();
+
+	for(saw::data<sch::UInt64> i{0u}; i < grid.template get_dim_size<0>(); ++i){
+		for(saw::data<sch::UInt64> j{0u}; j < grid.template get_dim_size<1>(); ++j){
+			std::cout<<grid.at({{i,j}}).get()<<" ";
+		}
+		std::cout<<"\n";
+	}
+	std::cout<<std::endl;
+
+	//saw::data<sch::Array<sch::T,2>> reference_mask{{{4+2,4+2}}};
+	//reference_mask.at({{0,0}});
+}
+*/
+/*
+SAW_TEST("Verlet integration test 2D"){
+	using namespace kel;
+	lbm::particle_system<sch::T,2,sch::Particle<sch::T,2>> system;
+
+	{
+		saw::data<sch::Particle<sch::T,2>> particle;
+		auto& rb = particle.template get<"rigid_body">();
+		auto& acc = rb.template get<"acceleration">();
+		auto& pos = rb.template get<"position">();
+		auto& pos_old = rb.template get<"position_old">();
+		pos = {{1e-1,1e-1}};
+		pos_old = {{0.0, 0.0}};
+		acc = {{0.0,-1e1}};
+
+		auto eov = system.add_particle(std::move(particle));
+		SAW_EXPECT(eov.is_value(), "Expected no error :)");
+	}
+	{
+		auto& p = system.at({0u});
+		auto& rb = p.template get<"rigid_body">();
+		auto& pos = rb.template get<"position">();
+
+		for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{2u}; ++i){
+			std::cout<<pos.at(i).get()<<" ";
+		}
+		std::cout<<std::endl;
+	}
+
+	for(uint64_t i = 0u; i < 36u; ++i){
+		system.step(saw::data<sch::T>{1e-1});
+
+		{
+			auto& p = system.at({0u});
+			auto& rb = p.template get<"rigid_body">();
+			auto& pos = rb.template get<"position">();
+
+			for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{2u}; ++i){
+				std::cout<<pos.at(i).get()<<" ";
+			}
+			std::cout<<"\n";
+
+			if(pos.at({1u}).get() < 0.0){
+				break;
+			}
+		}
+
+	}
+}
+*/
+
+SAW_TEST("Spheroid Particle / AABB"){
+	using namespace kel;
+
+	
+}
+}
diff --git a/modules/core/tests/vtk_write.cpp b/modules/core/tests/vtk_write.cpp
new file mode 100644
index 0000000..c99d752
--- /dev/null
+++ b/modules/core/tests/vtk_write.cpp
@@ -0,0 +1,93 @@
+#include <forstio/test/suite.hpp>
+
+#include <iostream>
+#include "../c++/write_vtk.hpp"
+
+#include <sstream>
+
+namespace {
+namespace sch {
+using namespace kel::lbm::sch;
+
+using T = Float64;
+using D2Q5 = Descriptor<2,5>;
+using D2Q9 = Descriptor<2,9>;
+
+template<typename Desc>
+using DfCell = Cell<T, Desc, 0, 0, 1>;
+
+template<typename Desc>
+using CellInfo = Cell<UInt8, D2Q9, 1, 0, 0>;
+
+/**
+ * Basic type for simulation
+ */
+template<typename Desc>
+using CellStruct = Struct<
+	Member<DfCell<Desc>, "dfs">,
+	Member<DfCell<Desc>, "dfs_old">,
+	Member<CellInfo<Desc>, "info">
+>;
+
+template<typename T, uint64_t D>
+using MacroStruct = Struct<
+	Member<FixedArray<T,D>, "velocity">,
+	Member<T, "pressure">
+>;
+
+}
+/*
+SAW_TEST("VTK Write test example"){
+	using namespace kel;
+
+	// write_vtk();
+
+	std::stringstream sstream;
+
+	saw::data<sch::Array<sch::MacroStruct<sch::T,2>, 2>> cells{{{2u,2u}}};
+
+	auto& cell_0 = cells.at({{{0,0}}});
+	cell_0.template get<"velocity">()= {{0.5,-0.1}};
+	cell_0.template get<"pressure">().set(1.1);
+
+	auto eov = lbm::impl::lbm_vtk_writer<sch::Array<sch::MacroStruct<sch::T,2>, 2>>::apply(sstream, cells);
+	SAW_EXPECT(eov.is_value(), "vtk writer failed to write");
+
+	// I want to print it to see it for myself. For now I have no tooling to more easily view associated and potentially generated files
+	std::cout<<sstream.str()<<std::endl;
+
+	static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS velocity float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n";
+	SAW_EXPECT(sstream.str() == comparison_str, "Expected different encoding");
+
+	// using Type = typename parameter_pack_type<i,T...>::type;
+}
+
+SAW_TEST("VTK Write raw test example"){
+	using namespace kel;
+
+	// write_vtk();
+
+	std::stringstream sstream;
+
+	saw::data<sch::Array<sch::MacroStruct<sch::T,2>, 2>> cells{{{2u,2u}}};
+
+	auto meta = cells.dims();
+
+	auto& cell_0 = cells.at({{{0,0}}});
+	cell_0.template get<"velocity">()= {{0.5,-0.1}};
+	cell_0.template get<"pressure">().set(1.1);
+
+	auto eov = lbm::impl::lbm_vtk_writer_raw<sch::MacroStruct<sch::T,2u>,2u>::apply(sstream, &cell_0, meta);
+	SAW_EXPECT(eov.is_value(), "vtk writer failed to write");
+
+	// I want to print it to see it for myself. For now I have no tooling to more easily view associated and potentially generated files
+	std::cout<<sstream.str()<<std::endl;
+
+	static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS velocity float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n";
+	SAW_EXPECT(sstream.str() == comparison_str, "Expected different encoding");
+
+	// using Type = typename parameter_pack_type<i,T...>::type;
+}
+*/
+
+}
diff --git a/modules/heterogenous/.nix/derivation.nix b/modules/heterogenous/.nix/derivation.nix
new file mode 100644
index 0000000..3686e21
--- /dev/null
+++ b/modules/heterogenous/.nix/derivation.nix
@@ -0,0 +1,41 @@
+{ lib
+, stdenv
+, scons
+, clang-tools
+, pname
+, version
+, forstio
+, kel-lbm
+, adaptive-cpp
+}:
+
+stdenv.mkDerivation {
+	pname = "${pname}-heterogenous";
+	inherit version;
+	src = ./..;
+
+	nativeBuildInputs = [
+		scons
+		clang-tools
+	];
+
+	buildInputs = [
+		forstio.core
+		forstio.async
+		forstio.codec
+		forstio.codec-unit
+		forstio.codec-json
+		forstio.remote
+		kel-lbm.core
+	];
+
+	doCheck = true;
+	checkPhase = ''
+		scons test
+		./bin/tests
+	'';
+
+	preferLocalBuild = true;
+
+	outputs = [ "out" "dev" ];
+}
diff --git a/modules/heterogenous/SConstruct b/modules/heterogenous/SConstruct
new file mode 100644
index 0000000..8b3ab01
--- /dev/null
+++ b/modules/heterogenous/SConstruct
@@ -0,0 +1,75 @@
+#!/usr/bin/env python3
+
+import sys
+import os
+import os.path
+import glob
+import re
+
+
+if sys.version_info < (3,):
+    def isbasestring(s):
+        return isinstance(s,basestring)
+else:
+    def isbasestring(s):
+        return isinstance(s, (str,bytes))
+
+def add_kel_source_files(self, sources, filetype, lib_env=None, shared=False, target_post=""):
+
+    if isbasestring(filetype):
+        dir_path = self.Dir('.').abspath
+        filetype = sorted(glob.glob(dir_path+"/"+filetype))
+
+    for path in filetype:
+        target_name = re.sub( r'(.*?)(\.cpp|\.c\+\+)', r'\1' + target_post, path )
+        if shared:
+            target_name+='.os'
+            sources.append( self.SharedObject( target=target_name, source=path ) )
+        else:
+            target_name+='.o'
+            sources.append( self.StaticObject( target=target_name, source=path ) )
+    pass
+
+def isAbsolutePath(key, dirname, env):
+	assert os.path.isabs(dirname), "%r must have absolute path syntax" % (key,)
+
+env_vars = Variables(
+	args=ARGUMENTS
+)
+
+env_vars.Add('prefix',
+	help='Installation target location of build results and headers',
+	default='/usr/local/',
+	validator=isAbsolutePath
+)
+
+env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[],
+    CPPDEFINES=['SAW_UNIX'],
+    CXXFLAGS=[
+        '-std=c++20',
+        '-g',
+        '-Wall',
+        '-Wextra'
+    ],
+    LIBS=[
+        'forstio-core'
+    ]
+);
+env.__class__.add_source_files = add_kel_source_files
+env.Tool('compilation_db');
+env.cdb = env.CompilationDatabase('compile_commands.json');
+
+env.objects = [];
+env.sources = [];
+env.headers = [];
+env.targets = [];
+
+Export('env')
+SConscript('c++/SConscript')
+SConscript('tests/SConscript')
+
+env.Alias('cdb', env.cdb);
+env.Alias('all', [env.targets]);
+env.Default('all');
+
+env.Alias('install', '$prefix')
diff --git a/modules/heterogenous/c++/SConscript b/modules/heterogenous/c++/SConscript
new file mode 100644
index 0000000..b66e979
--- /dev/null
+++ b/modules/heterogenous/c++/SConscript
@@ -0,0 +1,27 @@
+#!/bin/false
+
+import os
+import os.path
+import glob
+
+
+Import('env')
+
+dir_path = Dir('.').abspath
+
+# Environment for base library
+hetero_env = env.Clone();
+
+hetero_env.sources = sorted(glob.glob(dir_path + "/*.cpp"));
+hetero_env.headers = sorted(glob.glob(dir_path + "/*.hpp"));
+
+env.sources += hetero_env.sources;
+env.headers += hetero_env.headers;
+
+## Static lib
+objects = []
+hetero_env.add_source_files(objects, hetero_env.sources, shared=False);
+env.library_static = hetero_env.StaticLibrary('#build/kel-lbm-hetero', [objects]);
+
+env.Install('$prefix/lib/', env.library_static);
+env.Install('$prefix/include/kel/lbm/hetero/', hetero_env.headers);
diff --git a/modules/heterogenous/c++/common.hpp b/modules/heterogenous/c++/common.hpp
new file mode 100644
index 0000000..789c876
--- /dev/null
+++ b/modules/heterogenous/c++/common.hpp
@@ -0,0 +1,12 @@
+#pragma once
+
+#include <forstio/codec/data.hpp>
+
+namespace kel {
+namespace lbm {
+template<typename T, uint64_t D>
+class super_block {
+};
+
+}
+}
diff --git a/modules/heterogenous/tests/SConscript b/modules/heterogenous/tests/SConscript
new file mode 100644
index 0000000..d1b381e
--- /dev/null
+++ b/modules/heterogenous/tests/SConscript
@@ -0,0 +1,32 @@
+#!/bin/false
+
+import os
+import os.path
+import glob
+
+
+Import('env')
+
+dir_path = Dir('.').abspath
+
+# Environment for base library
+test_cases_env = env.Clone();
+
+test_cases_env.Append(LIBS=['forstio-test']);
+
+test_cases_env.sources = sorted(glob.glob(dir_path + "/*.cpp"))
+test_cases_env.headers = sorted(glob.glob(dir_path + "/*.hpp"))
+
+env.sources += test_cases_env.sources;
+env.headers += test_cases_env.headers;
+
+objects_static = []
+test_cases_env.add_source_files(objects_static, test_cases_env.sources, shared=False);
+test_cases_env.program = test_cases_env.Program('#bin/tests', [objects_static]);
+# , env.library_static]);
+
+# Set Alias
+env.Alias('test', test_cases_env.program);
+env.Alias('check', test_cases_env.program);
+
+env.targets += ['test','check'];
diff --git a/modules/png/.nix/derivation.nix b/modules/png/.nix/derivation.nix
new file mode 100644
index 0000000..6e52083
--- /dev/null
+++ b/modules/png/.nix/derivation.nix
@@ -0,0 +1,38 @@
+{ lib
+, stdenv
+, scons
+, clang-tools
+, pname
+, version
+, forstio
+, kel-lbm
+, adaptive-cpp
+}:
+
+stdenv.mkDerivation {
+	pname = "${pname}-png";
+	inherit version;
+	src = ./..;
+
+	nativeBuildInputs = [
+		scons
+		clang-tools
+	];
+
+	buildInputs = [
+		forstio.core
+		forstio.async
+		forstio.codec
+		kel-lbm.core
+	];
+
+	doCheck = true;
+	checkPhase = ''
+		scons test
+		./bin/tests
+	'';
+
+	preferLocalBuild = true;
+
+	outputs = [ "out" "dev" ];
+}
diff --git a/modules/png/SConstruct b/modules/png/SConstruct
new file mode 100644
index 0000000..8b3ab01
--- /dev/null
+++ b/modules/png/SConstruct
@@ -0,0 +1,75 @@
+#!/usr/bin/env python3
+
+import sys
+import os
+import os.path
+import glob
+import re
+
+
+if sys.version_info < (3,):
+    def isbasestring(s):
+        return isinstance(s,basestring)
+else:
+    def isbasestring(s):
+        return isinstance(s, (str,bytes))
+
+def add_kel_source_files(self, sources, filetype, lib_env=None, shared=False, target_post=""):
+
+    if isbasestring(filetype):
+        dir_path = self.Dir('.').abspath
+        filetype = sorted(glob.glob(dir_path+"/"+filetype))
+
+    for path in filetype:
+        target_name = re.sub( r'(.*?)(\.cpp|\.c\+\+)', r'\1' + target_post, path )
+        if shared:
+            target_name+='.os'
+            sources.append( self.SharedObject( target=target_name, source=path ) )
+        else:
+            target_name+='.o'
+            sources.append( self.StaticObject( target=target_name, source=path ) )
+    pass
+
+def isAbsolutePath(key, dirname, env):
+	assert os.path.isabs(dirname), "%r must have absolute path syntax" % (key,)
+
+env_vars = Variables(
+	args=ARGUMENTS
+)
+
+env_vars.Add('prefix',
+	help='Installation target location of build results and headers',
+	default='/usr/local/',
+	validator=isAbsolutePath
+)
+
+env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[],
+    CPPDEFINES=['SAW_UNIX'],
+    CXXFLAGS=[
+        '-std=c++20',
+        '-g',
+        '-Wall',
+        '-Wextra'
+    ],
+    LIBS=[
+        'forstio-core'
+    ]
+);
+env.__class__.add_source_files = add_kel_source_files
+env.Tool('compilation_db');
+env.cdb = env.CompilationDatabase('compile_commands.json');
+
+env.objects = [];
+env.sources = [];
+env.headers = [];
+env.targets = [];
+
+Export('env')
+SConscript('c++/SConscript')
+SConscript('tests/SConscript')
+
+env.Alias('cdb', env.cdb);
+env.Alias('all', [env.targets]);
+env.Default('all');
+
+env.Alias('install', '$prefix')
diff --git a/modules/png/c++/SConscript b/modules/png/c++/SConscript
new file mode 100644
index 0000000..75b8645
--- /dev/null
+++ b/modules/png/c++/SConscript
@@ -0,0 +1,27 @@
+#!/bin/false
+
+import os
+import os.path
+import glob
+
+
+Import('env')
+
+dir_path = Dir('.').abspath
+
+# Environment for base library
+png_env = env.Clone();
+
+png_env.sources = sorted(glob.glob(dir_path + "/*.cpp"));
+png_env.headers = sorted(glob.glob(dir_path + "/*.hpp"));
+
+env.sources += png_env.sources;
+env.headers += png_env.headers;
+
+## Static lib
+objects = []
+png_env.add_source_files(objects, png_env.sources, shared=False);
+env.library_static = png_env.StaticLibrary('#build/kel-lbm-png', [objects]);
+
+env.Install('$prefix/lib/', env.library_static);
+env.Install('$prefix/include/kel/lbm/png/', png_env.headers);
diff --git a/modules/png/c++/png.cpp b/modules/png/c++/png.cpp
new file mode 100644
index 0000000..248a8a5
--- /dev/null
+++ b/modules/png/c++/png.cpp
@@ -0,0 +1 @@
+#include "png.hpp"
diff --git a/modules/png/c++/png.hpp b/modules/png/c++/png.hpp
new file mode 100644
index 0000000..f675a99
--- /dev/null
+++ b/modules/png/c++/png.hpp
@@ -0,0 +1,7 @@
+#pragma once
+
+namespace kel {
+namespace lbm {
+
+}
+}
diff --git a/modules/sycl/.nix/derivation.nix b/modules/sycl/.nix/derivation.nix
new file mode 100644
index 0000000..4422739
--- /dev/null
+++ b/modules/sycl/.nix/derivation.nix
@@ -0,0 +1,43 @@
+{ lib
+, stdenv
+, scons
+, clang-tools
+, pname
+, version
+, forstio
+, kel
+, adaptive-cpp
+}:
+
+stdenv.mkDerivation {
+	pname = "${pname}-sycl";
+	inherit version;
+	src = ./..;
+
+	nativeBuildInputs = [
+		scons
+		clang-tools
+	];
+
+	buildInputs = [
+		forstio.core
+		forstio.async
+		forstio.codec
+		forstio.codec-unit
+		forstio.codec-json
+		forstio.remote
+		#forstio.remote-sycl
+		kel.lbm.core
+		adaptive-cpp
+	];
+
+	doCheck = true;
+	checkPhase = ''
+		scons test
+		./bin/tests
+	'';
+
+	preferLocalBuild = true;
+
+	outputs = [ "out" "dev" ];
+}
diff --git a/modules/sycl/SConstruct b/modules/sycl/SConstruct
new file mode 100644
index 0000000..f8bfbcb
--- /dev/null
+++ b/modules/sycl/SConstruct
@@ -0,0 +1,80 @@
+#!/usr/bin/env python3
+
+import sys
+import os
+import os.path
+import glob
+import re
+
+
+if sys.version_info < (3,):
+    def isbasestring(s):
+        return isinstance(s,basestring)
+else:
+    def isbasestring(s):
+        return isinstance(s, (str,bytes))
+
+def add_kel_source_files(self, sources, filetype, lib_env=None, shared=False, target_post=""):
+
+    if isbasestring(filetype):
+        dir_path = self.Dir('.').abspath
+        filetype = sorted(glob.glob(dir_path+"/"+filetype))
+
+    for path in filetype:
+        target_name = re.sub( r'(.*?)(\.cpp|\.c\+\+)', r'\1' + target_post, path )
+        if shared:
+            target_name+='.os'
+            sources.append( self.SharedObject( target=target_name, source=path ) )
+        else:
+            target_name+='.o'
+            sources.append( self.StaticObject( target=target_name, source=path ) )
+    pass
+
+def isAbsolutePath(key, dirname, env):
+	assert os.path.isabs(dirname), "%r must have absolute path syntax" % (key,)
+
+env_vars = Variables(
+	args=ARGUMENTS
+)
+
+env_vars.Add('prefix',
+	help='Installation target location of build results and headers',
+	default='/usr/local/',
+	validator=isAbsolutePath
+)
+
+env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[],
+    CPPDEFINES=['SAW_UNIX'],
+    CXX = ['syclcc-clang'],
+    CXXFLAGS=[
+        '-std=c++20',
+        '-g',
+        '-O3',
+        '-Wall',
+        '-Wextra',
+        '-isystem', 'AdaptiveCpp'
+    ],
+    LIBS=[
+        'forstio-core',
+        'acpp-rt',
+        'omp'
+    ]
+);
+env.__class__.add_source_files = add_kel_source_files
+env.Tool('compilation_db');
+env.cdb = env.CompilationDatabase('compile_commands.json');
+
+env.objects = [];
+env.sources = [];
+env.headers = [];
+env.targets = [];
+
+Export('env')
+SConscript('c++/SConscript')
+SConscript('tests/SConscript')
+
+env.Alias('cdb', env.cdb);
+env.Alias('all', [env.targets]);
+env.Default('all');
+
+env.Alias('install', '$prefix')
diff --git a/modules/sycl/c++/SConscript b/modules/sycl/c++/SConscript
new file mode 100644
index 0000000..2ed63ba
--- /dev/null
+++ b/modules/sycl/c++/SConscript
@@ -0,0 +1,27 @@
+#!/bin/false
+
+import os
+import os.path
+import glob
+
+
+Import('env')
+
+dir_path = Dir('.').abspath
+
+# Environment for base library
+core_env = env.Clone();
+
+core_env.sources = sorted(glob.glob(dir_path + "/*.cpp"));
+core_env.headers = sorted(glob.glob(dir_path + "/*.hpp"));
+
+env.sources += core_env.sources;
+env.headers += core_env.headers;
+
+## Static lib
+objects = []
+core_env.add_source_files(objects, core_env.sources, shared=False);
+env.library_static = core_env.StaticLibrary('#build/kel-lbm-sycl', [objects]);
+
+env.Install('$prefix/lib/', env.library_static);
+env.Install('$prefix/include/kel/lbm/sycl/', core_env.headers);
diff --git a/modules/sycl/c++/common.hpp b/modules/sycl/c++/common.hpp
new file mode 100644
index 0000000..8ff76fc
--- /dev/null
+++ b/modules/sycl/c++/common.hpp
@@ -0,0 +1,12 @@
+#pragma once
+
+#include <kel/lbm/common.hpp>
+#include <AdaptiveCpp/sycl/sycl.hpp>
+
+namespace kel {
+namespace lbm {
+namespace sycl {
+using namespace acpp::sycl;
+}
+}
+}
diff --git a/modules/sycl/c++/data.hpp b/modules/sycl/c++/data.hpp
new file mode 100644
index 0000000..9f43848
--- /dev/null
+++ b/modules/sycl/c++/data.hpp
@@ -0,0 +1,961 @@
+#include "common.hpp"
+#include <kel/lbm/lbm.hpp>
+
+namespace kel {
+namespace lbm {
+namespace encode {
+template<typename Encode>
+struct Sycl {
+};
+}
+
+/*
+namespace impl {
+template<typename Schema>
+struct struct_has_only_equal_dimension_array{};
+}
+*/
+}
+}
+
+namespace saw {
+template<typename Sch, uint64_t... Dims, typename Encode>
+class data<schema::FixedArray<Sch,Dims...>, kel::lbm::encode::Sycl<Encode>> final {
+public:
+	using Schema = schema::FixedArray<Sch,Dims...>;
+private:
+	acpp::sycl::queue* q_;
+	data<Sch,Encode>* values_;
+
+	SAW_FORBID_COPY(data);
+	SAW_FORBID_MOVE(data);
+public:
+	data(const data<typename meta_schema<Schema>::MetaSchema>& meta__, acpp::sycl::queue& q__):
+		q_{&q__},
+		values_{nullptr}
+	{
+		(void) meta__;
+		SAW_ASSERT(q_);
+		values_ = acpp::sycl::malloc_device<data<Sch,Encode>>(ct_multiply<uint64_t,Dims...>::value,*q_);
+		SAW_ASSERT(values_);
+	}
+
+	data(acpp::sycl::queue& q__):
+		q_{&q__},
+		values_{nullptr}
+	{
+		SAW_ASSERT(q_);
+		values_ = acpp::sycl::malloc_device<data<Sch,Encode>>(ct_multiply<uint64_t,Dims...>::value,*q_);
+		SAW_ASSERT(values_);
+	}
+
+	~data(){
+		if(not values_){
+			return;
+		}
+		SAW_ASSERT(q_);
+
+		acpp::sycl::free(values_,*q_);
+		values_ = nullptr;
+	}
+
+	static constexpr data<schema::FixedArray<schema::UInt64, sizeof...(Dims)>> get_dims() {
+		return saw::data<schema::FixedArray<schema::UInt64, sizeof...(Dims)>>{{Dims...}};
+	}
+
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index){
+		return values_[kel::lbm::flatten_index<schema::UInt64,sizeof...(Dims)>::apply(index,get_dims()).get()];
+	}
+	
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index) const{
+		return values_[kel::lbm::flatten_index<schema::UInt64,sizeof...(Dims)>::apply(index,get_dims()).get()];
+	}
+
+	constexpr data<Sch,Encode>* flat_data() const {
+		return values_;
+	}
+};
+
+template<typename Sch, uint64_t... Dims, typename Encode>
+class data<schema::Ptr<schema::FixedArray<Sch,Dims...>>, kel::lbm::encode::Sycl<Encode>> final {
+public:
+	using Schema = schema::Ptr<schema::FixedArray<Sch,Dims...>>;
+private:
+	data<Sch,Encode>* values_;
+
+public:
+	SAW_DEFAULT_COPY(data);
+	SAW_DEFAULT_MOVE(data);
+	
+	data():
+		values_{nullptr}
+	{}
+
+	data(data<schema::FixedArray<Sch,Dims...>, kel::lbm::encode::Sycl<Encode>>& values__):
+		values_{values__.flat_data()}
+	{}
+
+	data(data<Sch,Encode>* values__):
+		values_{values__}
+	{}
+
+	static constexpr data<schema::FixedArray<schema::UInt64, sizeof...(Dims)>> get_dims() {
+		return saw::data<schema::FixedArray<schema::UInt64, sizeof...(Dims)>>{{Dims...}};
+	}
+
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index){
+		return values_[kel::lbm::flatten_index<schema::UInt64,sizeof...(Dims)>::apply(index,get_dims()).get()];
+	}
+	
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index) const{
+		return values_[kel::lbm::flatten_index<schema::UInt64,sizeof...(Dims)>::apply(index,get_dims()).get()];
+	}
+
+	constexpr data<Sch,Encode>* flat_data() const {
+		return values_;
+	}
+};
+
+template<typename Sch, uint64_t Dims, typename Encode>
+class data<schema::Array<Sch,Dims>, kel::lbm::encode::Sycl<Encode>> final {
+public:
+	using Schema = schema::Array<Sch,Dims>;
+private:
+	static_assert(Dims > 0u, "Zero Dim Arrays make no sense here. If you meant to use this for math style approaches then use Tensor instead of Array");
+
+	data<Sch,Encode>* values_;
+	data<schema::FixedArray<schema::UInt64,Dims>,Encode> meta_;
+	acpp::sycl::queue* q_;
+
+	SAW_FORBID_COPY(data);
+	SAW_FORBID_MOVE(data);
+public:
+	data(const data<typename meta_schema<Schema>::MetaSchema>& meta__, acpp::sycl::queue& q__):
+		values_{nullptr},
+		meta_{meta__},
+		q_{&q__}
+	{
+		SAW_ASSERT(q_);
+		values_ = acpp::sycl::malloc_device<data<Sch,Encode>>(flat_size().get(),*q_);
+		SAW_ASSERT(values_);
+	}
+
+	data(acpp::sycl::queue& q__):
+		values_{nullptr},
+		meta_{},
+		q_{&q__}
+	{
+		SAW_ASSERT(q_);
+	}
+
+	~data(){
+		if(not values_){
+			return;
+		}
+		SAW_ASSERT(q_);
+
+		acpp::sycl::free(values_,*q_);
+		values_ = nullptr;
+	}
+
+	constexpr data<schema::FixedArray<schema::UInt64, Dims>, Encode> meta() const {
+		return meta_;
+	}
+
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,Dims>, Encode>& index){
+		return values_[kel::lbm::flatten_index<schema::UInt64,Dims>::apply(index,meta()).get()];
+	}
+
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,Dims>, Encode>& index) const{
+		return values_[kel::lbm::flatten_index<schema::UInt64,Dims>::apply(index,meta()).get()];
+	}
+
+	constexpr error_or<void> reset_to(const data<typename meta_schema<Schema>::MetaSchema>& meta_arg){
+		SAW_ASSERT(q_);
+		meta_ = meta_arg;
+
+		if(values_){
+			acpp::sycl::free(values_,*q_);
+		}
+		values_ = acpp::sycl::malloc_device<data<Sch,Encode>>(flat_size().get(),*q_);
+		SAW_ASSERT(q_);
+
+		return make_void();
+	}
+
+	constexpr data<Sch,Encode>* flat_data() const {
+		return values_;
+	}
+
+	constexpr data<schema::UInt64,Encode> flat_size() const {
+		data<schema::UInt64> mult{1u};
+
+		for(uint64_t i{0u}; i < Dims; ++i){
+			mult = mult * meta_.at({i});
+		}
+
+		return mult;
+	}
+};
+
+template<typename Sch, uint64_t Dims, typename Encode>
+class data<schema::Ptr<schema::Array<Sch,Dims>>, kel::lbm::encode::Sycl<Encode>> final {
+public:
+	using Schema = schema::Ptr<schema::Array<Sch,Dims>>;
+private:
+	static_assert(Dims > 0u, "Zero Dim Arrays make no sense here. If you meant to use this for math style approaches then use Tensor instead of Array");
+
+	data<Sch,Encode>* values_;
+	data<schema::FixedArray<schema::UInt64,Dims>,Encode> meta_;
+
+public:
+	SAW_DEFAULT_COPY(data);
+	SAW_DEFAULT_MOVE(data);
+
+	data():
+		values_{nullptr}
+	{}
+
+	data(const data<schema::Array<Sch,Dims>, kel::lbm::encode::Sycl<Encode>>& values__):
+		values_{values__.flat_data()},
+		meta_{values__.meta()}
+	{}
+
+	constexpr data<schema::FixedArray<schema::UInt64, Dims>, Encode> meta() const {
+		return meta_;
+	}
+
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,Dims>, Encode>& index){
+		return values_[kel::lbm::flatten_index<schema::UInt64,Dims>::apply(index,meta()).get()];
+	}
+
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,Dims>, Encode>& index) const{
+		return values_[kel::lbm::flatten_index<schema::UInt64,Dims>::apply(index,meta()).get()];
+	}
+
+	constexpr data<Sch,Encode>* flat_data() const {
+		return values_;
+	}
+
+	constexpr data<schema::UInt64,Encode> flat_size() const {
+		data<schema::UInt64> mult{1u};
+
+		for(uint64_t i{0u}; i < Dims; ++i){
+			mult = mult * meta_.at({i});
+		}
+
+		return mult;
+	}
+};
+
+template<typename Sch, uint64_t Ghost, uint64_t... Sides, typename Encode>
+class data<kel::lbm::sch::Chunk<Sch,Ghost,Sides...>,kel::lbm::encode::Sycl<Encode>> final {
+public:
+	using Schema = kel::lbm::sch::Chunk<Sch,Ghost,Sides...>;
+private:
+	using InnerSchema = typename Schema::InnerSchema;
+	using ValueSchema = typename InnerSchema::ValueType;
+
+	data<InnerSchema, kel::lbm::encode::Sycl<Encode>> values_;
+public:
+	data(const data<typename meta_schema<Schema>::MetaSchema>& meta__, acpp::sycl::queue& q__):
+		values_{meta__,q__}
+	{}
+
+	data(acpp::sycl::queue& q__):
+		values_{q__}
+	{}
+
+	constexpr data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
+		return values_.at(index);
+	}
+
+	constexpr data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
+		return values_.at(index);
+	}
+
+	static constexpr auto get_ghost_dims() {
+		return data<InnerSchema,kel::lbm::encode::Sycl<Encode>>::get_dims();
+	}
+
+	static constexpr auto ghost_meta() {
+		return data<InnerSchema,kel::lbm::encode::Sycl<Encode>>::meta();
+	}
+
+	data<ValueSchema, Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
+		std::decay_t<decltype(index)> ind;
+		for(uint64_t i = 0u; i < sizeof...(Sides); ++i){
+			ind.at({i}) = index.at({i}) + Ghost;
+		}
+		return values_.at(ind);
+	}
+
+	data<ValueSchema, Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
+		std::decay_t<decltype(index)> ind;
+		for(uint64_t i = 0u; i < sizeof...(Sides); ++i){
+			ind.at({i}) = index.at({i}) + Ghost;
+		}
+		return values_.at(ind);
+	}
+
+	static constexpr auto get_dims(){
+		return data<schema::FixedArray<schema::UInt64, sizeof...(Sides)>,Encode>{{Sides...}};
+	}
+
+	static constexpr auto meta(){
+		return data<schema::FixedArray<schema::UInt64, sizeof...(Sides)>,Encode>{{Sides...}};
+	}
+
+	auto flat_data() const {
+		return values_.flat_data();
+	}
+
+	static constexpr auto flat_size() {
+		return data<InnerSchema,kel::lbm::encode::Sycl<Encode>>::flat_size();
+	}
+};
+
+template<typename Sch, uint64_t Ghost, uint64_t... Sides, typename Encode>
+class data<schema::Ptr<kel::lbm::sch::Chunk<Sch,Ghost,Sides...>>,kel::lbm::encode::Sycl<Encode>> final {
+public:
+	using Schema = schema::Ptr<kel::lbm::sch::Chunk<Sch,Ghost,Sides...>>;
+private:
+	using InnerSchema = typename kel::lbm::sch::Chunk<Sch,Ghost,Sides...>::InnerSchema;
+	using ValueSchema = typename InnerSchema::ValueType;
+	
+	data<schema::Ptr<InnerSchema>,kel::lbm::encode::Sycl<Encode>> values_;
+
+public:
+	SAW_DEFAULT_MOVE(data);
+	SAW_DEFAULT_COPY(data);
+
+	data():
+		values_{nullptr}
+	{}
+
+	data(const data<kel::lbm::sch::Chunk<Sch,Ghost,Sides...>, kel::lbm::encode::Sycl<Encode>>& values__):
+		values_{values__.flat_data()}
+	{}
+
+	data(data<Sch,Encode>* values__):
+		values_{values__}
+	{}
+
+	constexpr data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
+		return values_.at(index);
+	}
+
+	constexpr data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
+		return values_.at(index);
+	}
+
+	static constexpr auto get_ghost_dims() {
+		return data<InnerSchema,kel::lbm::encode::Sycl<Encode>>::get_dims();
+	}
+
+	static constexpr auto ghost_meta() {
+		return data<InnerSchema,kel::lbm::encode::Sycl<Encode>>::meta();
+	}
+
+	data<ValueSchema, Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
+		std::decay_t<decltype(index)> ind;
+		for(uint64_t i = 0u; i < sizeof...(Sides); ++i){
+			ind.at({i}) = index.at({i}) + Ghost;
+		}
+		return values_.at(ind);
+	}
+
+	data<ValueSchema, Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
+		std::decay_t<decltype(index)> ind;
+		for(uint64_t i = 0u; i < sizeof...(Sides); ++i){
+			ind.at({i}) = index.at({i}) + Ghost;
+		}
+		return values_.at(ind);
+	}
+
+	static constexpr auto meta(){
+		return data<schema::FixedArray<schema::UInt64, sizeof...(Sides)>,Encode>{{Sides...}};
+	}
+
+	static constexpr auto get_dims(){
+		return data<schema::FixedArray<schema::UInt64, sizeof...(Sides)>,Encode>{{Sides...}};
+	}
+
+	auto flat_data() const {
+		return values_.flat_data();
+	}
+
+	static constexpr auto flat_size() {
+		return data<InnerSchema,kel::lbm::encode::Sycl<Encode>>::flat_size();
+	}
+};
+
+template<typename... Members, typename Encode>
+struct data<schema::Tuple<Members...>, kel::lbm::encode::Sycl<Encode>> final {
+public:
+	using StorageT = std::tuple<data<Members,kel::lbm::encode::Sycl<Encode>>...>;
+	using Schema = schema::Tuple<Members...>;
+private:
+	StorageT members_;
+
+	/**
+	 * A helper constructor to forward the sycl queue to the inner "default" constructors
+	 */
+	template<std::size_t... Is>
+	constexpr data(acpp::sycl::queue& q, std::index_sequence<Is...>):
+		members_{(static_cast<void>(Is), q)...}
+	{}
+public:
+	/*
+	data(data<typename meta_schema<Schema>::MetaSchema>& meta__, acpp::sycl::queue& q__):
+		data{q__, std::make_index_sequence<sizeof...(Members)>{}}
+	{}
+	*/
+
+	data(acpp::sycl::queue& q__):
+		data{q__, std::make_index_sequence<sizeof...(Members)>{}}
+	{
+		q__.wait();
+	}
+
+	template<size_t i>
+	auto& get(){
+		return std::get<i>(members_);
+	}
+
+	template<size_t i>
+	auto& get() const {
+		return std::get<i>(members_);
+	}
+};
+
+template<typename... Members, typename Encode>
+struct data<schema::Ptr<schema::Tuple<Members...>>, kel::lbm::encode::Sycl<Encode>> final {
+public:
+	using StorageT = std::tuple<data<schema::Ptr<Members>,kel::lbm::encode::Sycl<Encode>>...>;
+	using Schema = schema::Tuple<Members...>;
+private:
+	StorageT members_;
+
+	/**
+	 * A helper constructor to forward the sycl queue to the inner "default" constructors
+	 */
+public:
+	data() = default;
+
+	template<size_t i>
+	auto& get(){
+		return std::get<i>(members_);
+	}
+
+	template<size_t i>
+	const auto& get() const {
+		return std::get<i>(members_);
+	}
+};
+
+template<typename... Members, typename Encode>
+class data<schema::Struct<Members...>, kel::lbm::encode::Sycl<Encode> > final {
+public:
+	using StorageT = std::tuple<data<typename Members::ValueType,kel::lbm::encode::Sycl<Encode>>...>;
+	using Schema = schema::Struct<Members...>;
+private:
+	/**
+	 * @todo Check by static assert that the members all have the same dimensions. Alternatively
+	 * Do it here by specializing.
+	 */
+	StorageT members_;
+
+	/**
+	 * A helper constructor to forward the sycl queue to the inner "default" constructors
+	 */
+	template<std::size_t... Is>
+	constexpr data(acpp::sycl::queue& q, std::index_sequence<Is...>):
+		members_{(static_cast<void>(Is), q)...}
+	{}
+public:
+	data(acpp::sycl::queue& q__):
+		data{q__, std::make_index_sequence<sizeof...(Members)>{}}
+	{
+		q__.wait();
+	}
+
+	template<size_t i>
+	auto& get(){
+		return std::get<i>(members_);
+	}
+
+	template<size_t i>
+	auto& get() const {
+		return std::get<i>(members_);
+	}
+
+	template<saw::string_literal K>
+	auto& get(){
+		return std::get<parameter_key_pack_index<K, Members::KeyLiteral...>::value>(members_);
+	}
+	
+	template<saw::string_literal K>
+	auto& get() const {
+		return std::get<parameter_key_pack_index<K, Members::KeyLiteral...>::value>(members_);
+	}
+};
+
+template<typename... Sch, saw::string_literal... Keys, typename Encode>
+class data<schema::Ptr<schema::Struct<schema::Member<Sch,Keys>...>>, kel::lbm::encode::Sycl<Encode> > final {
+public:
+	using StorageT = std::tuple<data<schema::Ptr<Sch>,kel::lbm::encode::Sycl<Encode>>...>;
+	using Schema = schema::Struct<schema::Member<schema::Ptr<Sch>,Keys>...>;
+private:
+	/**
+	 * @todo Check by static assert that the members all have the same dimensions. Alternatively
+	 * Do it here by specializing.
+	 */
+	StorageT members_;
+public:
+	data() = default;
+
+	template<size_t i>
+	auto& get(){
+		return std::get<i>(members_);
+	}
+
+	template<size_t i>
+	auto& get() const {
+		return std::get<i>(members_);
+	}
+
+	template<saw::string_literal K>
+	auto& get(){
+		return std::get<parameter_key_pack_index<K, Keys...>::value>(members_);
+	}
+	
+	template<saw::string_literal K>
+	auto& get() const {
+		return std::get<parameter_key_pack_index<K, Keys...>::value>(members_);
+	}
+};
+}
+
+namespace kel {
+namespace lbm {
+namespace impl {
+template<typename Sch, typename Encode>
+struct sycl_copy_helper;
+
+template<typename... Members, typename Encode>
+struct sycl_copy_helper<sch::Struct<Members...>, Encode> final {
+	using Schema = sch::Struct<Members...>;
+
+	template<uint64_t i>
+	static saw::error_or<void> copy_to_device_member(saw::data<Schema,Encode>& host_data, saw::data<Schema,encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+		if constexpr (i < sizeof...(Members)){
+			using M = typename saw::parameter_pack_type<i,Members...>::type;
+			auto& host_member_data = host_data.template get<M::KeyLiteral>();
+			auto& sycl_member_data = sycl_data.template get<M::KeyLiteral>();
+
+			auto eov = sycl_copy_helper<typename M::ValueType,Encode>::copy_to_device(host_member_data,sycl_member_data,q);
+			if(eov.is_error()){
+				return eov;
+			}
+
+			return copy_to_device_member<i+1u>(host_data,sycl_data,q);
+		}
+		
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> copy_to_device(saw::data<Schema,Encode>& host_data, saw::data<Schema, encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+
+		return copy_to_device_member<0u>(host_data, sycl_data, q);
+	}
+
+	template<uint64_t i>
+	static saw::error_or<void> copy_to_host_member(saw::data<Schema,encode::Sycl<Encode>>& sycl_data, saw::data<Schema,Encode>& host_data, sycl::queue& q){
+		if constexpr (i < sizeof...(Members)){
+			using M = typename saw::parameter_pack_type<i,Members...>::type;
+			auto& host_member_data = host_data.template get<M::KeyLiteral>();
+			auto& sycl_member_data = sycl_data.template get<M::KeyLiteral>();
+
+			auto eov = sycl_copy_helper<typename M::ValueType,Encode>::copy_to_host(sycl_member_data,host_member_data,q);
+			if(eov.is_error()){
+				return eov;
+			}
+
+			return copy_to_host_member<i+1u>(sycl_data,host_data,q);
+		}
+		
+		return saw::make_void();
+	}
+	
+
+	static saw::error_or<void> copy_to_host(saw::data<Schema,encode::Sycl<Encode>>& sycl_data, saw::data<Schema,Encode>& host_data, sycl::queue& q){
+		return copy_to_host_member<0u>(sycl_data, host_data, q);
+	}
+
+	template<uint64_t i>
+	static saw::error_or<void> malloc_on_device_member(saw::data<Schema,Encode>& host_data, saw::data<Schema,encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+
+		if constexpr (i < sizeof...(Members)){
+			using M = typename saw::parameter_pack_type<i,Members...>::type;
+			auto& host_member_data = host_data.template get<i>();
+			auto& sycl_member_data = sycl_data.template get<i>();
+
+			auto eov = sycl_copy_helper<typename M::ValueType,Encode>::malloc_on_device(host_member_data,sycl_member_data,q);
+			if(eov.is_error()){
+				return eov;
+			}
+
+			return malloc_on_device_member<i+1u>(host_data,sycl_data,q);
+		}
+
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> malloc_on_device(saw::data<Schema,Encode>& host_data, saw::data<Schema,encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+		return malloc_on_device_member<0u>(host_data,sycl_data,q);
+	}
+};
+
+template<typename... Members, typename Encode>
+struct sycl_copy_helper<sch::Tuple<Members...>, Encode> final {
+	using Schema = sch::Tuple<Members...>;
+
+	template<uint64_t i>
+	static saw::error_or<void> copy_to_device_member(saw::data<Schema,Encode>& host_data, saw::data<Schema,encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+		if constexpr (i < sizeof...(Members)){
+			using M = typename saw::parameter_pack_type<i,Members...>::type;
+			auto& host_member_data = host_data.template get<i>();
+			auto& sycl_member_data = sycl_data.template get<i>();
+
+			auto eov = sycl_copy_helper<M,Encode>::copy_to_device(host_member_data,sycl_member_data,q);
+			if(eov.is_error()){
+				return eov;
+			}
+
+			return copy_to_device_member<i+1u>(host_data,sycl_data,q);
+		}
+
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> copy_to_device(saw::data<Schema,Encode>& host_data, saw::data<Schema, encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+		return copy_to_device_member<0u>(host_data, sycl_data, q);
+	}
+
+	template<uint64_t i>
+	static saw::error_or<void> copy_to_host_member(saw::data<Schema,encode::Sycl<Encode>>& sycl_data, saw::data<Schema,Encode>& host_data, sycl::queue& q){
+		if constexpr (i < sizeof...(Members)){
+			using M = typename saw::parameter_pack_type<i,Members...>::type;
+
+			auto& host_member_data = host_data.template get<i>();
+			auto& sycl_member_data = sycl_data.template get<i>();
+
+			auto eov = sycl_copy_helper<M,Encode>::copy_to_host(sycl_member_data,host_member_data,q);
+			if(eov.is_error()){
+				return eov;
+			}
+
+			return copy_to_host_member<i+1u>(sycl_data,host_data,q);
+		}
+
+		return saw::make_void();
+	}
+
+
+	static saw::error_or<void> copy_to_host(
+		saw::data<Schema,Encode>& host_data,
+		saw::data<Schema,encode::Sycl<Encode>>& sycl_data,
+		sycl::queue& q
+	){
+		return copy_to_host_member<0u>(sycl_data, host_data, q);
+	}
+
+	template<uint64_t i>
+	static saw::error_or<void> malloc_on_device_member(saw::data<Schema,Encode>& host_data, saw::data<Schema,encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+
+		if constexpr (i < sizeof...(Members)){
+			using M = typename saw::parameter_pack_type<i,Members...>::type;
+			auto& host_member_data = host_data.template get<i>();
+			auto& sycl_member_data = sycl_data.template get<i>();
+
+			auto eov = sycl_copy_helper<M,Encode>::malloc_on_device(host_member_data,sycl_member_data,q);
+			if(eov.is_error()){
+				return eov;
+			}
+
+			return malloc_on_device_member<i+1u>(host_data,sycl_data,q);
+		}
+
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> malloc_on_device(saw::data<Schema,Encode>& host_data, saw::data<Schema,encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+		return malloc_on_device_member<0u>(host_data,sycl_data,q);
+	}
+};
+
+template<typename Sch, uint64_t... Dims, typename Encode>
+struct sycl_copy_helper<sch::FixedArray<Sch,Dims...>, Encode> final {
+	using Schema = sch::FixedArray<Sch,Dims...>;
+	
+	static saw::error_or<void> copy_to_host(saw::data<Schema,encode::Sycl<Encode>>& sycl_data, saw::data<Schema,Encode>& host_data, sycl::queue& q){
+		auto host_ptr = host_data.flat_data();
+		auto sycl_ptr = sycl_data.flat_data();
+
+		static_assert(sizeof(std::decay_t<decltype(sycl_ptr)>) == sizeof(std::decay_t<decltype(host_ptr)>), "Unequal size");
+
+		q.submit([&](acpp::sycl::handler& h){
+			h.copy(sycl_ptr,host_ptr, saw::ct_multiply<uint64_t,Dims...>::value);
+		}).wait();
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> copy_to_device(saw::data<Schema,Encode>& host_data, saw::data<Schema,encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+		auto host_ptr = host_data.flat_data();
+		auto sycl_ptr = sycl_data.flat_data();
+
+		static_assert(sizeof(std::decay_t<decltype(sycl_ptr)>) == sizeof(std::decay_t<decltype(host_ptr)>), "Unequal size");
+
+		q.submit([&](acpp::sycl::handler& h){
+			h.copy(host_ptr,sycl_ptr, saw::ct_multiply<uint64_t,Dims...>::value);
+		}).wait();
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> malloc_on_device(saw::data<Schema,Encode>& host_data, saw::data<Schema,encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+		(void) host_data;
+		(void) sycl_data;
+		(void) q;
+		return saw::make_void();
+	}
+};
+
+template<typename Sch, uint64_t Ghost, uint64_t... Dims, typename Encode>
+struct sycl_copy_helper<sch::Chunk<Sch,Ghost,Dims...>, Encode> final {
+	using Schema = sch::Chunk<Sch,Ghost,Dims...>;
+
+	static saw::error_or<void> copy_to_host(saw::data<Schema,encode::Sycl<Encode>>& sycl_data, saw::data<Schema,Encode>& host_data, sycl::queue& q){
+		auto host_ptr = host_data.flat_data();
+		auto sycl_ptr = sycl_data.flat_data();
+
+		static_assert(sizeof(std::decay_t<decltype(sycl_ptr)>) == sizeof(std::decay_t<decltype(host_ptr)>), "Unequal size");
+
+		auto flat_size = host_data.flat_size();
+
+		q.submit([&](acpp::sycl::handler& h){
+			h.copy(sycl_ptr,host_ptr, flat_size.get());
+		}).wait();
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> copy_to_device(saw::data<Schema,Encode>& host_data, saw::data<Schema,encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+		auto host_ptr = host_data.flat_data();
+		auto sycl_ptr = sycl_data.flat_data();
+
+		static_assert(sizeof(std::decay_t<decltype(sycl_ptr)>) == sizeof(std::decay_t<decltype(host_ptr)>), "Unequal size");
+
+		auto flat_size = host_data.flat_size();
+
+		q.submit([&](acpp::sycl::handler& h){
+			h.copy(host_ptr,sycl_ptr, flat_size.get());
+		}).wait();
+		return saw::make_void();
+	}
+};
+
+template<typename Sch, uint64_t Dims, typename Encode>
+struct sycl_copy_helper<sch::Array<Sch,Dims>, Encode> final {
+	using Schema = sch::Array<Sch,Dims>;
+
+	static saw::error_or<void> copy_to_host(saw::data<Schema,encode::Sycl<Encode>>& sycl_data, saw::data<Schema,Encode>& host_data, sycl::queue& q){
+		auto host_ptr = host_data.flat_data();
+		auto sycl_ptr = sycl_data.flat_data();
+
+		static_assert(sizeof(std::decay_t<decltype(sycl_ptr)>) == sizeof(std::decay_t<decltype(host_ptr)>), "Unequal size");
+
+		SAW_ASSERT(host_data.flat_size() == sycl_data.flat_size());
+		q.submit([&](acpp::sycl::handler& h){
+			h.copy(sycl_ptr,host_ptr, host_data.flat_size().get());
+		}).wait();
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> copy_to_device(saw::data<Schema,Encode>& host_data, saw::data<Schema,encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+
+		{
+			auto hm = host_data.meta();
+			auto sm = sycl_data.meta();
+			bool equ{true};
+			for(uint64_t i{0u}; i < Dims; ++i){
+				equ &= (hm.at({i}).get() == sm.at({i}).get());
+			}
+			if(not equ){
+				sycl_data.reset_to(hm);
+			}
+		}
+
+		auto host_ptr = host_data.flat_data();
+		auto sycl_ptr = sycl_data.flat_data();
+		static_assert(sizeof(std::decay_t<decltype(sycl_ptr)>) == sizeof(std::decay_t<decltype(host_ptr)>), "Unequal size");
+
+		q.submit([&](acpp::sycl::handler& h){
+			h.copy(host_ptr,sycl_ptr, host_data.flat_size().get());
+		}).wait();
+
+		return saw::make_void();
+	}
+
+	static saw::error_or<void> malloc_on_device(saw::data<Schema,Encode>& host_data, saw::data<Schema,encode::Sycl<Encode>>& sycl_data, sycl::queue& q){
+		sycl_data = {host_data.meta(),q};
+		return saw::make_void();
+	}
+};
+
+
+template<typename Schema, typename Encode>
+struct make_view_helper;
+
+template<typename Sch, uint64_t Dims, typename Encode>
+struct make_view_helper<sch::Array<Sch,Dims>, encode::Sycl<Encode>> {
+public:
+static saw::error_or<void> apply(
+	saw::data<sch::Array<Sch,Dims>, encode::Sycl<Encode>>& dat,
+	saw::data<sch::Ptr<sch::Array<Sch,Dims>>, encode::Sycl<Encode>>& dat_view
+){
+	dat_view = {dat};
+	return saw::make_void();
+}
+};
+
+template<typename Sch, uint64_t... Dims, typename Encode>
+struct make_view_helper<sch::FixedArray<Sch,Dims...>, encode::Sycl<Encode>> {
+public:
+	static saw::error_or<void> apply(
+		saw::data<sch::FixedArray<Sch,Dims...>, encode::Sycl<Encode>>& dat,
+		saw::data<sch::Ptr<sch::FixedArray<Sch,Dims...>>, encode::Sycl<Encode>>& dat_view
+	){
+		dat_view = {dat};
+		return saw::make_void();
+	}
+};
+
+template<typename Sch, uint64_t Ghost, uint64_t... Dims, typename Encode>
+struct make_view_helper<sch::Chunk<Sch,Ghost,Dims...>, encode::Sycl<Encode>> {
+public:
+	static saw::error_or<void> apply(
+		saw::data<sch::Chunk<Sch,Ghost,Dims...>, encode::Sycl<Encode>>& dat,
+		saw::data<sch::Ptr<sch::Chunk<Sch,Ghost,Dims...>>, encode::Sycl<Encode>>& dat_view
+	){
+		dat_view = {dat};
+		return saw::make_void();
+	}
+};
+
+template<typename... Sch, saw::string_literal... Keys, typename Encode>
+struct make_view_helper<sch::Struct<sch::Member<Sch,Keys>...>, encode::Sycl<Encode>> {
+private:
+template<uint64_t i>
+static saw::error_or<void> apply_i(
+		saw::data<sch::Struct<sch::Member<Sch,Keys>...>, encode::Sycl<Encode>>& dat,
+		saw::data<sch::Ptr<sch::Struct<sch::Member<Sch,Keys>...>>, encode::Sycl<Encode>>& dat_view
+	){
+		if constexpr (i < sizeof...(Sch)){
+			using M = typename saw::parameter_pack_type<i,sch::Member<Sch,Keys>...>::type;
+
+			auto eov = make_view_helper<typename M::ValueType,encode::Sycl<Encode>>::apply(dat.template get<M::KeyLiteral>(),dat_view.template get<M::KeyLiteral>());
+			if(eov.is_error()){
+				return eov;
+			}
+
+			return apply_i<i+1u>(dat,dat_view);
+		}
+
+		return saw::make_void();
+	}
+public:
+static saw::error_or<void> apply(
+	saw::data<sch::Struct<sch::Member<Sch,Keys>...>, encode::Sycl<Encode>>& dat,
+	saw::data<sch::Ptr<sch::Struct<sch::Member<Sch,Keys>...>>, encode::Sycl<Encode>>& dat_view
+){
+	return apply_i<0u>(dat,dat_view);
+}
+};
+
+template<typename... Sch, typename Encode>
+struct make_view_helper<sch::Tuple<Sch...>, encode::Sycl<Encode>> {
+private:
+	template<uint64_t i>
+	static saw::error_or<void> apply_i(
+		saw::data<sch::Tuple<Sch...>, encode::Sycl<Encode>>& dat,
+		saw::data<sch::Ptr<sch::Tuple<Sch...>>, encode::Sycl<Encode>>& dat_view
+	){
+		if constexpr (i < sizeof...(Sch)){
+			using M = typename saw::parameter_pack_type<i,Sch...>::type;
+
+			auto eov = make_view_helper<M,encode::Sycl<Encode>>::apply(dat.template get<i>(), dat_view.template get<i>());
+			if(eov.is_error()){
+				return eov;
+			}
+
+			return apply_i<i+1u>(dat,dat_view);
+		}
+
+		return saw::make_void();
+	}
+public:
+	static saw::error_or<void> apply(
+		saw::data<sch::Tuple<Sch...>, encode::Sycl<Encode>>& dat,
+		saw::data<sch::Ptr<sch::Tuple<Sch...>>, encode::Sycl<Encode>> dat_view
+	){
+		return apply_i<0u>(dat,dat_view);
+	}
+};
+}
+
+// Ptr => Ptr<Chunk<T,Ghost,Dims...>>> => Ptr<FixedArray<T,Dims+Ghost...>
+template<typename Schema,typename Encode>
+auto make_view(saw::data<Schema,Encode>& dat){
+	saw::data<sch::Ptr<Schema>,Encode> dat_view;
+	auto eov = impl::make_view_helper<Schema,Encode>::apply(dat,dat_view);
+	(void) eov;
+
+	return dat_view;
+}
+
+class device final {
+private:
+	sycl::queue q_;
+
+	SAW_FORBID_COPY(device);
+	SAW_FORBID_MOVE(device);
+public:
+	device() = default;
+	~device() = default;
+
+	template<typename Sch, typename Encode>
+	saw::error_or<void> copy_to_device(saw::data<Sch,Encode>& host_data, saw::data<Sch,encode::Sycl<Encode>>& sycl_data){
+		return impl::sycl_copy_helper<Sch,Encode>::copy_to_device(host_data, sycl_data, q_);
+	}
+
+	template<typename Sch, typename Encode>
+	saw::error_or<void> copy_to_host(saw::data<Sch,encode::Sycl<Encode>>& sycl_data, saw::data<Sch,Encode>& host_data){
+		return impl::sycl_copy_helper<Sch,Encode>::copy_to_host(sycl_data, host_data, q_);
+	}
+
+	template<typename Sch, typename Encode>
+	saw::error_or<void> malloc_on_device(
+			saw::data<Sch,Encode>& host_data,
+			saw::data<Sch,encode::Sycl<Encode>>& sycl_data
+			){
+		auto eov = impl::sycl_copy_helper<Sch,Encode>::malloc_on_device(host_data, sycl_data, q_);
+		q_.wait();
+		return eov;
+	}
+
+	auto& get_handle(){
+		return q_;
+	}
+};
+}
+}
diff --git a/modules/sycl/c++/index.hpp b/modules/sycl/c++/index.hpp
new file mode 100644
index 0000000..0d2c035
--- /dev/null
+++ b/modules/sycl/c++/index.hpp
@@ -0,0 +1,18 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+template<typename D>
+saw::data<sch::FixedArray<sch::UInt64,D>> sycl_to_saw_index(const acpp::sycl::id<D>& idx){
+	saw::data<sch::FixedArray<sch::UInt64,D>> trans_index;
+
+	for(uint64_t i{0u}; i < D; ++i){
+		trans_index.at({{i}}).set(idx[i]);
+	}
+
+	return trans_index;
+}
+}
+}
diff --git a/modules/sycl/c++/lbm.hpp b/modules/sycl/c++/lbm.hpp
new file mode 100644
index 0000000..3e67ae6
--- /dev/null
+++ b/modules/sycl/c++/lbm.hpp
@@ -0,0 +1,3 @@
+#pragma once
+
+#include "data.hpp"
diff --git a/modules/sycl/c++/sycl.hpp b/modules/sycl/c++/sycl.hpp
new file mode 100644
index 0000000..973b511
--- /dev/null
+++ b/modules/sycl/c++/sycl.hpp
@@ -0,0 +1,5 @@
+#pragma once
+
+#include "common.hpp"
+#include "data.hpp"
+#include "index.hpp"
diff --git a/modules/sycl/tests/SConscript b/modules/sycl/tests/SConscript
new file mode 100644
index 0000000..d1b381e
--- /dev/null
+++ b/modules/sycl/tests/SConscript
@@ -0,0 +1,32 @@
+#!/bin/false
+
+import os
+import os.path
+import glob
+
+
+Import('env')
+
+dir_path = Dir('.').abspath
+
+# Environment for base library
+test_cases_env = env.Clone();
+
+test_cases_env.Append(LIBS=['forstio-test']);
+
+test_cases_env.sources = sorted(glob.glob(dir_path + "/*.cpp"))
+test_cases_env.headers = sorted(glob.glob(dir_path + "/*.hpp"))
+
+env.sources += test_cases_env.sources;
+env.headers += test_cases_env.headers;
+
+objects_static = []
+test_cases_env.add_source_files(objects_static, test_cases_env.sources, shared=False);
+test_cases_env.program = test_cases_env.Program('#bin/tests', [objects_static]);
+# , env.library_static]);
+
+# Set Alias
+env.Alias('test', test_cases_env.program);
+env.Alias('check', test_cases_env.program);
+
+env.targets += ['test','check'];
diff --git a/modules/sycl/tests/data.cpp b/modules/sycl/tests/data.cpp
new file mode 100644
index 0000000..4321a0d
--- /dev/null
+++ b/modules/sycl/tests/data.cpp
@@ -0,0 +1,56 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/lbm.hpp"
+
+namespace {
+
+namespace sch {
+using namespace kel::lbm::sch;
+using TestObjSchema = Tuple<
+	Member<FixedArray<UInt64,2u,2u>, "foo">,
+	Member<Array<Float32>, "bar">,
+	Member<
+		Array<
+			Struct<
+				Member<FixedArray<Float32,2u>,"pos">
+			>
+		>,
+		"baz"
+	>
+>;
+}
+
+SAW_TEST("Sycl Data Compilation"){
+	acpp::sycl::queue q;
+	saw::data<
+		saw::schema::Struct<
+			saw::schema::Member<
+				kel::lbm::sch::Chunk<saw::schema::UInt8,1u,1u,1u>,
+				"test"
+			>
+		>,
+		kel::lbm::encode::Sycl<saw::encode::Native>
+	> dat{q};
+
+	auto& test_f = dat.template get<"test">();
+
+	// test_f.at({}).set(1);
+	// SAW_EXPECT(test_f.at({}).get() == 1, "Value check failed");
+}
+
+SAW_TEST("Sycl Data Array of Struct"){
+	acpp::sycl::queue q;
+
+	saw::data<sch::Array<sch::Float64>, kel::lbm::encode::Sycl<saw::encode::Native>> a{{{2u}},q};
+}
+
+/*
+SAW_TEST("Sycl Data Compilation for Particle Similacrum"){
+	acpp::sycl::queue q;
+
+	saw::data<
+		sch::TestObjSchema
+	> a;
+}
+*/
+}