Moving project structure for more less compilation

22 files changed, 0 insertions, 1626 deletions

diff --git a/c++/SConscript b/c++/SConscript
deleted file mode 100644
index 85a078f..0000000
--- a/c++/SConscript
+++ /dev/null
@@ -1,32 +0,0 @@
-#!/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/c++/args.hpp b/c++/args.hpp
deleted file mode 100644
index 99c6172..0000000
--- a/c++/args.hpp
+++ /dev/null
@@ -1,12 +0,0 @@
-#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/c++/boundary.hpp b/c++/boundary.hpp
deleted file mode 100644
index 4e30f71..0000000
--- a/c++/boundary.hpp
+++ /dev/null
@@ -1,154 +0,0 @@
-#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>
-class component<T, Descriptor, cmpt::BounceBack> {
-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>& 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>
-class component<FP, Descriptor, cmpt::ZouHeHorizontal<Dir>> 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>& 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/c++/collision.hpp b/c++/collision.hpp
deleted file mode 100644
index 47870c1..0000000
--- a/c++/collision.hpp
+++ /dev/null
@@ -1,129 +0,0 @@
-#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>
-class component<T, Descriptor, cmpt::BGK> {
-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>& 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>
-class component<T, Descriptor, cmpt::BGKGuo> {
-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>& 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/c++/component.hpp b/c++/component.hpp
deleted file mode 100644
index 87ac31c..0000000
--- a/c++/component.hpp
+++ /dev/null
@@ -1,38 +0,0 @@
-#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>
-class component {};
-}
-}
diff --git a/c++/config.hpp b/c++/config.hpp
deleted file mode 100644
index 64f7a0f..0000000
--- a/c++/config.hpp
+++ /dev/null
@@ -1,53 +0,0 @@
-#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/c++/converter.hpp b/c++/converter.hpp
deleted file mode 100644
index 5c19c68..0000000
--- a/c++/converter.hpp
+++ /dev/null
@@ -1,79 +0,0 @@
-#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/c++/descriptor.hpp b/c++/descriptor.hpp
deleted file mode 100644
index d04d217..0000000
--- a/c++/descriptor.hpp
+++ /dev/null
@@ -1,296 +0,0 @@
-#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;
-
-template<typename Desc, typename... CellFieldTypes, saw::string_literal... CellFieldNames>
-struct CellField<
-	Desc,
-	Struct<
-		Member<CellFieldTypes, CellFieldNames>...
-	>
->;
-
-}
-
-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);
-	}
-};
-}
diff --git a/c++/environment.hpp b/c++/environment.hpp
deleted file mode 100644
index 4915e3a..0000000
--- a/c++/environment.hpp
+++ /dev/null
@@ -1,24 +0,0 @@
-#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/c++/equilibrium.hpp b/c++/equilibrium.hpp
deleted file mode 100644
index bb55d00..0000000
--- a/c++/equilibrium.hpp
+++ /dev/null
@@ -1,49 +0,0 @@
-#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/c++/geometry.hpp b/c++/geometry.hpp
deleted file mode 100644
index 9802feb..0000000
--- a/c++/geometry.hpp
+++ /dev/null
@@ -1,14 +0,0 @@
-#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/c++/hlbm.hpp b/c++/hlbm.hpp
deleted file mode 100644
index 1c665ce..0000000
--- a/c++/hlbm.hpp
+++ /dev/null
@@ -1,24 +0,0 @@
-#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/c++/iterator.hpp b/c++/iterator.hpp
deleted file mode 100644
index 78babff..0000000
--- a/c++/iterator.hpp
+++ /dev/null
@@ -1,32 +0,0 @@
-#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/c++/lbm.hpp b/c++/lbm.hpp
deleted file mode 100644
index 39c42af..0000000
--- a/c++/lbm.hpp
+++ /dev/null
@@ -1,34 +0,0 @@
-#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/c++/lbm_unit.hpp b/c++/lbm_unit.hpp
deleted file mode 100644
index 2d90652..0000000
--- a/c++/lbm_unit.hpp
+++ /dev/null
@@ -1,70 +0,0 @@
-#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/c++/macroscopic.hpp b/c++/macroscopic.hpp
deleted file mode 100644
index 51368e9..0000000
--- a/c++/macroscopic.hpp
+++ /dev/null
@@ -1,92 +0,0 @@
-#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/c++/particle/geometry/circle.hpp b/c++/particle/geometry/circle.hpp
deleted file mode 100644
index 331f06d..0000000
--- a/c++/particle/geometry/circle.hpp
+++ /dev/null
@@ -1,66 +0,0 @@
-#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/c++/particle/particle.hpp b/c++/particle/particle.hpp
deleted file mode 100644
index 39aadfb..0000000
--- a/c++/particle/particle.hpp
+++ /dev/null
@@ -1,113 +0,0 @@
-#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/c++/statistics.hpp b/c++/statistics.hpp
deleted file mode 100644
index c07ccb7..0000000
--- a/c++/statistics.hpp
+++ /dev/null
@@ -1,11 +0,0 @@
-#pragma once
-
-namespace kel {
-namespace lbm {
-template<typename T>
-class statistics {
-private:
-public:
-};
-}
-}
diff --git a/c++/term_renderer.hpp b/c++/term_renderer.hpp
deleted file mode 100644
index 5cbb551..0000000
--- a/c++/term_renderer.hpp
+++ /dev/null
@@ -1,6 +0,0 @@
-#pragma once
-
-namespace kel {
-namespace lbm {
-}
-}
diff --git a/c++/util.hpp b/c++/util.hpp
deleted file mode 100644
index 0bdebd1..0000000
--- a/c++/util.hpp
+++ /dev/null
@@ -1,93 +0,0 @@
-#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/c++/write_vtk.hpp b/c++/write_vtk.hpp
deleted file mode 100644
index 0647db5..0000000
--- a/c++/write_vtk.hpp
+++ /dev/null
@@ -1,205 +0,0 @@
-#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;
-}
-}
-}