Moving project structure for more less compilation

23 files changed, 1722 insertions, 0 deletions

diff --git a/lib/c++/SConscript b/lib/c++/SConscript
new file mode 100644
index 0000000..85a078f
--- /dev/null
+++ b/lib/c++/SConscript
@@ -0,0 +1,32 @@
+#!/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"));
+
+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);
diff --git a/lib/c++/args.hpp b/lib/c++/args.hpp
new file mode 100644
index 0000000..99c6172
--- /dev/null
+++ b/lib/c++/args.hpp
@@ -0,0 +1,12 @@
+#pragma once
+
+#include <forstio/codec/data.hpp>
+#include <forstio/codec/args.hpp>
+
+namespace kel {
+namespace lbm {
+namespace sch {
+using namespace saw::schema;
+}
+}
+}
diff --git a/lib/c++/boundary.hpp b/lib/c++/boundary.hpp
new file mode 100644
index 0000000..5cc7705
--- /dev/null
+++ b/lib/c++/boundary.hpp
@@ -0,0 +1,154 @@
+#pragma once
+
+#include "component.hpp"
+#include "equilibrium.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct BounceBack {};
+
+template<bool East>
+struct ZouHeHorizontal{};
+
+template<bool North>
+struct ZouHeVertical{};
+}
+
+/**
+ * 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(saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, uint64_t time_step){
+		bool is_even = ((time_step % 2u) == 0u);
+
+		// This is a ref
+		auto& cell = field(index);
+
+		auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
+		// auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
+
+		using dfi = df_info<T,Descriptor>;
+
+		// Technically use .copy()
+		auto df_cpy = dfs_old.copy();
+
+		for(uint64_t i = 1u; i < Descriptor::Q; ++i){
+			dfs_old({i}) = df_cpy({dfi::opposite_index.at(i)});
+		}
+	}
+};
+
+/**
+ * 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> pressure_setting_;
+	saw::data<FP> rho_setting_;
+public:
+	component(const saw::data<FP>& pressure_setting__):
+		pressure_setting_{pressure_setting__},
+		rho_setting_{pressure_setting__ * df_info<FP,Descriptor>::inv_cs2}
+	{}
+
+	template<typename CellFieldSchema>
+	void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, uint64_t time_step){
+		using dfi = df_info<FP,Descriptor>;
+
+		bool is_even = ((time_step % 2) == 0);
+		auto& cell = field(index);
+
+		auto& info = cell.template get<"info">();
+		if(info({0u}).get() == 0u){
+			return;
+		}
+		auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
+		// auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
+
+		/**
+		 * Sum all known DFs
+		 */
+		saw::data<FP> sum_df{0};
+		for(saw::data<sch::UInt64> k{0u}; k < saw::data<sch::UInt64>{Descriptor::Q}; ++k){
+			auto c_k = dfi::directions[k.get()];
+			auto& cell_n = field({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}});
+			auto& info_n = cell_n.template get<"info">();
+			auto info_n_val = info_n({0u});
+			auto k_opp = dfi::opposite_index[k.get()];
+
+			if(info_n_val.get() > 0u){
+				sum_df += dfs_old({k_opp});
+			}
+		}
+
+		/**
+		 * Get the sum of the unknown dfs and precalculate the direction
+		 */
+		constexpr int known_dir = Dir ? 1 : -1;
+		auto sum_unknown_dfs = (rho_setting_ - sum_df) * saw::data<FP>{known_dir};
+
+		for(saw::data<sch::UInt64> k{0u}; k < saw::data<sch::UInt64>{Descriptor::Q}; ++k){
+			auto c_k = dfi::directions[k.get()];
+			auto& cell_n = field({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}});
+			auto& info_n = cell_n.template get<"info">();
+			auto info_n_val = info_n({0u});
+			// auto k_opp = dfi::opposite_index[k.get()];
+
+			if(info_n_val.get() > 0u){
+				sum_unknown_dfs += dfs_old({k}) * c_k[0u];
+			}
+		}
+
+		auto vel_x = sum_unknown_dfs / rho_setting_;
+
+		static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly");
+
+		if constexpr (Dir) {
+			dfs_old({2u}) = dfs_old({1u}) + saw::data<FP>{2.0 / 3.0} * rho_setting_ * vel_x;
+			dfs_old({6u}) = dfs_old({5u}) + saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old({4u}) - dfs_old({3u}));
+			dfs_old({8u}) = dfs_old({7u}) + saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old({3u}) - dfs_old({4u}));
+		}else if constexpr (not Dir){
+			dfs_old({1u}) = dfs_old({2u}) - saw::data<FP>{2.0 / 3.0} * rho_setting_ * vel_x;
+			dfs_old({5u}) = dfs_old({6u}) - saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old({3u}) - dfs_old({4u}));
+			dfs_old({7u}) = dfs_old({8u}) - saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old({4u}) - dfs_old({3u}));
+		}
+	}
+};
+}
+}
diff --git a/lib/c++/collision.hpp b/lib/c++/collision.hpp
new file mode 100644
index 0000000..9ab542b
--- /dev/null
+++ b/lib/c++/collision.hpp
@@ -0,0 +1,129 @@
+#pragma once
+
+#include "macroscopic.hpp"
+#include "component.hpp"
+#include "equilibrium.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:
+	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_}
+	{}
+
+	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(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, uint64_t time_step){
+		bool is_even = ((time_step % 2) == 0);
+		auto& cell = field(index);
+
+		auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
+		// auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
+
+		saw::data<T> rho;
+		saw::data<sch::FixedArray<T,Descriptor::D>> vel;
+		compute_rho_u<T,Descriptor>(dfs_old,rho,vel);
+		auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			dfs_old({i}) = dfs_old({i}) + frequency_ * (eq.at(i) - dfs_old({i}));
+			// dfs_old({i}).set(dfs_old({i}).get() + (1.0 / relaxation_) * (eq.at(i).get() - dfs_old({i}).get()));
+		}
+	}
+};
+
+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_;
+public:
+	component(typename saw::native_data_type<T>::type relaxation__):
+		relaxation_{relaxation__},
+		frequency_{typename saw::native_data_type<T>::type(1) / relaxation_}
+	{}
+
+	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>
+	void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, uint64_t time_step){
+		bool is_even = ((time_step % 2) == 0);
+		auto& cell = field(index);
+
+		auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
+		// auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
+
+		saw::data<T> rho;
+		saw::data<sch::FixedArray<T,Descriptor::D>> vel;
+		compute_rho_u<T,Descriptor>(dfs_old,rho,vel);
+		auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+		using dfi = df_info<T,Descriptor>;
+
+		auto& force = cell.template get<"force">();
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			// saw::data<T> ci_min_u{0};
+			saw::data<T> ci_dot_u{0};
+			for(uint64_t d = 0u; d < Descriptor::D; ++d){
+				ci_dot_u = ci_dot_u + vel.at({d}) * saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])};
+			}
+			saw::data<T> F_i{0};// = saw::data<T>{dfi::weights[i]} * (ci_dot_F * dfi::inv_cs2 + ci_dot_u * ci_dot_F * (dfi::inv_cs2 * dfi::inv_cs2));
+
+			for(uint64_t d = 0u; d < Descriptor::D; ++d){
+				F_i = F_i +
+					saw::data<T>{dfi::weights[i]} * ((saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])}
+					- vel.at({d}) ) * dfi::inv_cs2 + ci_dot_u * saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])} * dfi::inv_cs2  * dfi::inv_cs2 ) * force({d});
+			}
+
+
+			dfs_old({i}) = dfs_old({i}) + frequency_ * (eq.at(i) - dfs_old({i}) ) + F_i * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_);
+		}
+	}
+};
+}
+}
diff --git a/lib/c++/component.hpp b/lib/c++/component.hpp
new file mode 100644
index 0000000..1e5dbbf
--- /dev/null
+++ b/lib/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/lib/c++/config.hpp b/lib/c++/config.hpp
new file mode 100644
index 0000000..64f7a0f
--- /dev/null
+++ b/lib/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/lib/c++/converter.hpp b/lib/c++/converter.hpp
new file mode 100644
index 0000000..5c19c68
--- /dev/null
+++ b/lib/c++/converter.hpp
@@ -0,0 +1,79 @@
+#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};
+	}
+
+	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/lib/c++/descriptor.hpp b/lib/c++/descriptor.hpp
new file mode 100644
index 0000000..51d5814
--- /dev/null
+++ b/lib/c++/descriptor.hpp
@@ -0,0 +1,365 @@
+#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;
+};
+
+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 Schema = CellFieldStruct<Desc, Struct<CellFieldMembers...>>;
+	using InnerSchema = Struct<CellFieldMembers...>;
+	using MetaSchema = FixedArray<UInt64, Desc::D>;
+};
+
+}
+
+template<typename T, typename Desc>
+class df_info{};
+
+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},
+		{-1, 0},
+		{ 1, 0},
+		{ 0,-1},
+		{ 0, 1},
+		{-1,-1},
+		{ 1, 1},
+		{-1, 1},
+		{ 1,-1}
+	}};
+
+	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},
+		{-1, 0, 0},
+		{ 1, 0, 0},
+		{ 0,-1, 0},
+		{ 0, 1, 0},
+		{-1,-1, 0},
+		{ 1, 1, 0},
+		{-1, 1, 0},
+		{ 1,-1, 0}
+	}};
+
+	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 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 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;
+	}
+};
+
+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_.get_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);
+	}
+};
+
+/**
+ * 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::Name>() = {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::Name>().dims();
+	}
+
+	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/lib/c++/environment.hpp b/lib/c++/environment.hpp
new file mode 100644
index 0000000..4915e3a
--- /dev/null
+++ b/lib/c++/environment.hpp
@@ -0,0 +1,24 @@
+#pragma once
+
+#include <forstio/error.hpp>
+#include <filesystem>
+#include <cstdlib>
+
+namespace kel {
+namespace lbm {
+
+/**
+ * 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(){
+	const char* home_dir = std::getenv("HOME");
+	if(not home_dir){
+		return saw::make_error<saw::err::not_found>("Couldn't find home dir");
+	}
+
+	return std::filesystem::path{home_dir} / ".lbm";
+}
+}
+}
diff --git a/lib/c++/equilibrium.hpp b/lib/c++/equilibrium.hpp
new file mode 100644
index 0000000..bb55d00
--- /dev/null
+++ b/lib/c++/equilibrium.hpp
@@ -0,0 +1,49 @@
+#pragma once
+
+#include "descriptor.hpp"
+
+namespace kel {
+namespace lbm {
+template<typename T, typename Descriptor>
+saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<T> rho, saw::data<sch::FixedArray<T,Descriptor::D>> vel){
+	using dfi = df_info<T, Descriptor>;
+
+	saw::data<sch::FixedArray<T,Descriptor::Q>> eq;
+	//   ^
+	//  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 < eq.template get_dim_size<0u>(); ++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.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/lib/c++/geometry.hpp b/lib/c++/geometry.hpp
new file mode 100644
index 0000000..9802feb
--- /dev/null
+++ b/lib/c++/geometry.hpp
@@ -0,0 +1,14 @@
+#pragma once
+
+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){
+
+	}
+};
+*/
+}
+}
diff --git a/lib/c++/hlbm.hpp b/lib/c++/hlbm.hpp
new file mode 100644
index 0000000..1c665ce
--- /dev/null
+++ b/lib/c++/hlbm.hpp
@@ -0,0 +1,24 @@
+#pragma once
+
+#include "macroscopic.hpp"
+#include "component.hpp"
+#include "equilibrium.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct HLBM {};
+}
+
+/**
+ * HLBM collision operator for LBM
+ */
+template<typename T, typename Descriptor>
+class component<T, Descriptor, cmpt::HLBM> {
+private:
+public:
+	component() = default;
+};
+
+}
+}
diff --git a/lib/c++/iterator.hpp b/lib/c++/iterator.hpp
new file mode 100644
index 0000000..78babff
--- /dev/null
+++ b/lib/c++/iterator.hpp
@@ -0,0 +1,32 @@
+#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;
+}
+/* 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/lib/c++/lbm.hpp b/lib/c++/lbm.hpp
new file mode 100644
index 0000000..39c42af
--- /dev/null
+++ b/lib/c++/lbm.hpp
@@ -0,0 +1,34 @@
+#pragma once
+
+#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 "iterator.hpp"
+#include "macroscopic.hpp"
+#include "write_vtk.hpp"
+#include "util.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){
+	std::cout
+		<<"[LBM Meta]\n"
+		<<"==========\n"
+		<<"\n"
+		<<"Δx:     "<<conv.delta_x()<<"\n"
+		<<"Δt:     "<<conv.delta_t()<<"\n"
+		<<"KinVis: "<<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"
+	;
+}
+}
+}
diff --git a/lib/c++/lbm_unit.hpp b/lib/c++/lbm_unit.hpp
new file mode 100644
index 0000000..2d90652
--- /dev/null
+++ b/lib/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/lib/c++/macroscopic.hpp b/lib/c++/macroscopic.hpp
new file mode 100644
index 0000000..51368e9
--- /dev/null
+++ b/lib/c++/macroscopic.hpp
@@ -0,0 +1,92 @@
+#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::Cell<T, Desc, 0, 0, 1>>& dfs,
+		typename saw::native_data_type<T>::type& rho,
+		std::array<typename saw::native_data_type<T>::type, 2>& vel
+	)
+{
+	using dfi = df_info<T, Desc>;
+
+	rho = 0;
+	std::fill(vel.begin(), vel.end(), 0);
+
+	for(size_t j = 0; j < Desc::Q; ++j){
+		rho += dfs(j).get();
+		for(size_t i = 0; i < Desc::D; ++i){
+			vel[i] += dfi::directions[j][i] * dfs(j).get();
+		}
+	}
+
+	for(size_t i = 0; i < Desc::D; ++i){
+		vel[i] /= rho;
+	}
+}
+
+/**
+ * Calculate the macroscopic variables rho and u in Lattice Units.
+ */
+template<typename T, typename Desc>
+void compute_rho_u (
+		const saw::data<sch::Cell<T, Desc, 0, 0, 1>>& dfs,
+		saw::ref<saw::data<T>> rho,
+		saw::ref<saw::data<sch::FixedArray<T,Desc::D>>> vel
+	)
+{
+	using dfi = df_info<T, Desc>;
+
+	rho().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() = rho() + dfs(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(j);
+		}
+	}
+
+	for(size_t i = 0; i < Desc::D; ++i){
+		vel().at({i}) = vel().at({i}) / rho();
+	}
+}
+
+/**
+ * Calculate the macroscopic variables rho and u in Lattice Units.
+ */
+template<typename T, typename Desc>
+void compute_rho_u (
+		const saw::data<sch::Cell<T, Desc, 0, 0, 1>>& dfs,
+		saw::ref<saw::data<T>> rho,
+		saw::ref<saw::data<sch::Vector<T,Desc::D>>> vel
+	)
+{
+	using dfi = df_info<T, Desc>;
+
+	rho().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() = rho() + dfs(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(j);
+		}
+	}
+
+	for(size_t i = 0; i < Desc::D; ++i){
+		vel().at({{i}}) = vel().at({{i}}) / rho();
+	}
+}
+}
+}
diff --git a/lib/c++/particle/geometry/circle.hpp b/lib/c++/particle/geometry/circle.hpp
new file mode 100644
index 0000000..331f06d
--- /dev/null
+++ b/lib/c++/particle/geometry/circle.hpp
@@ -0,0 +1,66 @@
+#pragma once
+
+#include "../particle.hpp"
+
+namespace kel {
+namespace lbm {
+
+template<typename T>
+class particle_circle_geometry {
+private:
+public:
+	particle_circle_geometry()
+	{}
+
+	template<typename MT = T>
+	saw::data<sch::ParticleMask<MT,2>> generate_mask(uint64_t resolution, uint64_t boundary_nodes = 0) const {
+		saw::data<sch::ParticleMask<MT,2>> 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});
+
+		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);
+					}
+				}
+			}
+		}
+
+		saw::data<T> total_mass{};
+		iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+			auto ind_vec = saw::math::vectorize_data(index).template cast_to<T>();
+			for(uint64_t i{0u}; i < 2u; ++i){
+				ind_vec.at({{i}}) = ind_vec.at({{i}}) * grid.at(index);
+			}
+			com = com + ind_vec;
+
+			total_mass = total_mass + grid.at(index);
+		},{{0u,0u}}, grid.dims());
+
+		for(uint64_t i{0u}; i < 2u; ++i){
+			com.at({{i}}) = com.at({{i}}) / total_mass;
+		}
+
+		return mask;
+	}
+};
+
+}
+}
diff --git a/lib/c++/particle/particle.hpp b/lib/c++/particle/particle.hpp
new file mode 100644
index 0000000..39aadfb
--- /dev/null
+++ b/lib/c++/particle/particle.hpp
@@ -0,0 +1,113 @@
+#pragma once
+
+#include <forstio/codec/data.hpp>
+#include <forstio/codec/data_math.hpp>
+#include <forstio/codec/math.hpp>
+
+namespace kel {
+namespace lbm {
+namespace sch {
+using namespace saw::schema;
+
+template<typename T, uint64_t D>
+using ParticleRigidBody = Struct<
+	Member<Vector<T,D>, "position">,
+	Member<Vector<T,D>, "position_old">,
+	Member<T, "rotation">,
+	Member<T, "rotation_old">,
+
+	Member<Vector<T,D>, "acceleration">,
+	Member<T, "rotational_acceleration">
+>;
+
+template<typename T, uint64_t D>
+using ParticleMask = Struct<
+	Member<Array<T,D>, "grid">,
+	Member<Vector<T,D>, "center_of_mass">
+>;
+
+template<typename T, uint64_t D>
+using Particle = Struct<
+	Member<ParticleRigidBody<T,D>, "rigid_body">,
+	Member<ParticleMask<T,D>, "mask">,
+	Member<T, "mass">,
+	Member<T, "size">
+>;
+}
+
+template<typename T, uint64_t D, typename ParticleCollision = sch::ParticleMask<T,D> >
+class particle_system {
+private:
+	saw::data<sch::Array<sch::Particle<T,D>>> particles_;
+
+	void verlet_step(saw::data<sch::Particle<T,D>>& particle, saw::data<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& rot = body.template get<"rotation">();
+		auto& acc = body.template get<"acceleration">();
+
+		auto tsd_squared = time_step_delta * time_step_delta;
+
+		saw::data<sch::Vector<T,D>> pos_new;
+		// Actual step
+		for(uint64_t i = 0u; i < D; ++i){
+			pos_new.at({{i}}) = saw::data<T>{2.0} * pos.at({{i}}) - pos_old.at({{i}}) + acc.at({{i}}) * tsd_squared;
+		}
+
+		pos_old = pos;
+		pos = pos_new;
+	}
+public:
+	/**
+	 *  Add particle to this class and return an id representing this particle
+	 */
+	saw::error_or<saw::data<sch::UInt64>> add_particle(saw::data<sch::Particle<T,D>> particle__){
+		auto size = particles_.size();
+		auto eov = particles_.add(std::move(particle__));
+		if(eov.is_error()){
+			return std::move(eov.get_error());
+		}
+
+		return size;
+	}
+
+	/*
+	saw::data<sch::Particle<T,D>>& get_particle(saw::data<sch::UInt64> id){
+	}
+	*/
+
+	void step(saw::data<T> time_step_delta){
+		for(saw::data<sch::UInt64> i{0u}; i < particles_.size(); ++i){
+			verlet_step(particles_.at(i), time_step_delta);
+		}
+	}
+
+	template<typename LbmLattice>
+	void update_particle_border(saw::data<LbmLattice>& latt){
+		(void) latt;
+		for(auto& iter : particles_){
+			auto& par = iter;
+
+			auto& body = par.template get<"rigid_body">();
+			auto& size = par.template get<"size">();
+
+
+		}
+	}
+
+	saw::data<sch::UInt64> size() const {
+		return particles_.size();
+	}
+
+	/**
+	 * Mostly meant for unforeseen use cases.
+	 */
+	saw::data<sch::Particle<T,D>>& at(saw::data<sch::UInt64> i){
+		return particles_.at(i);
+	}
+};
+}
+}
diff --git a/lib/c++/statistics.hpp b/lib/c++/statistics.hpp
new file mode 100644
index 0000000..c07ccb7
--- /dev/null
+++ b/lib/c++/statistics.hpp
@@ -0,0 +1,11 @@
+#pragma once
+
+namespace kel {
+namespace lbm {
+template<typename T>
+class statistics {
+private:
+public:
+};
+}
+}
diff --git a/lib/c++/term_renderer.hpp b/lib/c++/term_renderer.hpp
new file mode 100644
index 0000000..5cbb551
--- /dev/null
+++ b/lib/c++/term_renderer.hpp
@@ -0,0 +1,6 @@
+#pragma once
+
+namespace kel {
+namespace lbm {
+}
+}
diff --git a/lib/c++/util.hpp b/lib/c++/util.hpp
new file mode 100644
index 0000000..0bdebd1
--- /dev/null
+++ b/lib/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/lib/c++/write_json.hpp b/lib/c++/write_json.hpp
new file mode 100644
index 0000000..19bbfa6
--- /dev/null
+++ b/lib/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/lib/c++/write_vtk.hpp b/lib/c++/write_vtk.hpp
new file mode 100644
index 0000000..0647db5
--- /dev/null
+++ b/lib/c++/write_vtk.hpp
@@ -0,0 +1,205 @@
+#pragma once
+
+#include <forstio/error.hpp>
+
+#include <forstio/codec/data.hpp>
+#include <forstio/codec/data_math.hpp>
+
+#include "descriptor.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\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){
+		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){
+		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, 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<uint64_t Dim, typename... StructT, saw::string_literal... StructN>
+struct lbm_vtk_writer<sch::Array<sch::Struct<sch::Member<StructT,StructN>...> , Dim>> {
+
+	template<uint64_t i, uint64_t Dep>
+	static saw::error_or<void> write_i_iterate_d(std::ostream& vtk_file, const saw::data<sch::Array<sch::Struct<sch::Member<StructT,StructN>...>,Dim>>& field, saw::data<sch::FixedArray<sch::UInt64,Dim>>& index){
+		constexpr auto Lit = saw::parameter_key_pack_type<i, StructN...>::literal;
+		using Type = typename saw::parameter_pack_type<i,StructT...>::type;
+
+		if constexpr (Dep == 0u){
+			return lbm_vtk_writer<Type>::apply(vtk_file, field.at(index).template get<Lit>());
+		} else {
+			// Dep < Dim // I hope
+			static_assert(Dep > 0u, "Don't fall into this case");
+			for(index.at({Dep-1u}) = 0; index.at({Dep-1u}) < field.get_dims().at({Dep-1u}); ++index.at({Dep-1u})){
+				auto eov = write_i_iterate_d<i,Dep-1u>(vtk_file, field, index);
+				if(eov.is_error()){
+					return eov;
+				}
+			}
+		}
+		return saw::make_void();
+	}
+
+	template<uint64_t i>
+	static saw::error_or<void> write_i(std::ostream& vtk_file, const saw::data<sch::Array<sch::Struct<sch::Member<StructT,StructN>...>,Dim>>& field){
+
+		// auto meta = field.get_dims();
+		saw::data<sch::FixedArray<sch::UInt64,Dim>> index;
+		for(saw::data<sch::UInt64> it{0}; it.get() < Dim; ++it){
+			index.at({0u}).set(0u);
+		}
+		// vtk_file write?
+		// Data header
+		{
+
+			auto eov = lbm_vtk_writer<typename saw::parameter_pack_type<i,StructT...>::type>::apply_header(vtk_file, saw::parameter_key_pack_type<i, StructN...>::literal.view());
+			if(eov.is_error()){
+				return eov;
+			}
+
+		}
+		return write_i_iterate_d<i,Dim>(vtk_file, field, index);
+	}
+
+	template<uint64_t i>
+	static saw::error_or<void> iterate_i(std::ostream& vtk_file, const saw::data<sch::Array<sch::Struct<sch::Member<StructT,StructN>...>, Dim>>& field){
+		// constexpr auto Lit = saw::parameter_key_pack_type<i, StructN...>::literal;
+		{
+			auto eov = write_i<i>(vtk_file, field);
+			if(eov.is_error()){
+				return eov;
+			}
+		}
+		if constexpr ( (i+1u) < sizeof...(StructT) ){
+			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::Array<sch::Struct<sch::Member<StructT,StructN>...>, Dim>>& 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 meta = field.get_dims();
+		saw::data<sch::UInt64> pd_size{1u};
+		// DIMENSIONS
+
+		{
+			vtk_file << "DIMENSIONS";
+			for(saw::data<sch::UInt64> i{0u}; i.get() < Dim; ++i){
+				pd_size = pd_size * meta.at(i);
+				vtk_file << " " << meta.at(i).get(); 
+			}
+			for(saw::data<sch::UInt64> i{Dim}; i.get() < 3u; ++i){
+				vtk_file << " 1";
+			}
+
+			vtk_file << "\n";
+		}
+
+		if constexpr (sizeof...(StructT) > 0u){
+			// POINT DATA
+			{
+				vtk_file << "POINT_DATA " << pd_size.get() <<"\n";
+			}
+
+			// HEADER TO BODY
+			{
+				vtk_file << "\n";
+			}
+
+			return iterate_i<0u>(vtk_file, field);
+		}
+
+		return saw::make_void();
+	}
+};
+}
+
+template<typename Struct>
+saw::error_or<void> write_vtk_file(const std::filesystem::path& file_name, const saw::data<Struct>& field){
+
+	std::ofstream vtk_file{file_name};
+
+	if(!vtk_file.is_open()){
+		return saw::make_error<saw::err::critical>("Could not open file.");
+	}
+
+
+	auto eov = impl::lbm_vtk_writer<Struct>::apply(vtk_file, field);
+	return eov;
+}
+}
+}