Moving project structure for more less compilation

42 files changed, 503 insertions, 114 deletions

diff --git a/.nix/derivation.nix b/.nix/derivation.nix
index 917c21c..5d77b56 100644
--- a/.nix/derivation.nix
+++ b/.nix/derivation.nix
@@ -7,8 +7,8 @@
 
 stdenv.mkDerivation {
   pname = "kel-lbm";
-  version = "0.0.1";
-  src = ./..;
+  version = "0.0.2";
+  src = ../lib;
 
   nativeBuildInputs = [
     scons
diff --git a/default.nix b/default.nix
index 52458aa..b6f8f5b 100644
--- a/default.nix
+++ b/default.nix
@@ -37,6 +37,11 @@ in rec {
 			inherit pname version stdenv forstio;
 			kel-lbm = lbm;
 		};
+
+		meta_2d = pkgs.callPackage ./examples/meta_2d/.nix/derivation.nix {
+			inherit pname version stdenv forstio;
+			kel-lbm = lbm;
+		};
 	};
 
 	debug = {
diff --git a/examples/cavity_2d_gpu/.nix/derivation.nix b/examples/cavity_2d_gpu/.nix/derivation.nix
index f8d7314..73b81b7 100644
--- a/examples/cavity_2d_gpu/.nix/derivation.nix
+++ b/examples/cavity_2d_gpu/.nix/derivation.nix
@@ -24,6 +24,8 @@ stdenv.mkDerivation {
     forstio.async
     forstio.codec
     forstio.codec-unit
+    forstio.remote
+    forstio.remote-sycl
 		forstio.codec-json
 		adaptive-cpp
 		kel-lbm
diff --git a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
index cc0a811..c4836d0 100644
--- a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
+++ b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
@@ -1,7 +1,6 @@
 #include <kel/lbm/lbm.hpp>
 
 #include <forstio/codec/data.hpp>
-// #include <forstio/remote/
 
 #include <AdaptiveCpp/sycl/sycl.hpp>
 
@@ -9,6 +8,83 @@
 #include <fstream>
 #include <cmath>
 
+namespace saw {
+namespace encode {
+template<typename InnerEnc>
+struct Sycl {
+};
+}
+template<typename Schema>
+class data<Schema, encode::Sycl<encode::Native>> {
+private:
+	cl::sycl::buffer<data<Schema, encode::Native>> data_;
+	data<schema::UInt64, encode::Native> size_;
+public:
+	data(const data<Schema, encode::Native>& data__):
+		data_{&data__, 1u},
+		size_{data__.size()}
+	{}
+
+	auto& get_handle() {
+		return data_;
+	}
+
+	const auto& get_handle() const {
+		return data_;
+	}
+
+	data<schema::UInt64, encode::Native> size() const {
+		return size_;
+	}
+
+	template<cl::sycl::access::mode AccessMode>
+	auto access(cl::sycl::handler& h){
+		return data_.template get_access<AccessMode>(h);
+	}
+
+	template<cl::sycl::access::mode AccessMode>
+	auto access(cl::sycl::handler& h) const {
+		return data_.template get_access<AccessMode>(h);
+	}
+};
+
+template<typename Sch, uint64_t Dim>
+class data<schema::Array<Sch, Dim>, encode::Sycl<encode::Native>> {
+public:
+	using Schema = schema::Array<Sch,Dim>;
+private:
+	cl::sycl::buffer<data<Sch, encode::Native>> data_;
+	data<schema::UInt64, encode::Native> size_;
+public:
+	data(const data<Schema, encode::Native>& host_data__):
+		data_{&host_data__.at({0u}),host_data__.size().get()},
+		size_{host_data__.size()}
+	{}
+
+	auto& get_handle() {
+		return data_;
+	}
+
+	const auto& get_handle() const {
+		return data_;
+	}
+
+	data<schema::UInt64, encode::Native> size() const {
+		return size_;
+	}
+
+	template<cl::sycl::access::mode AccessMode>
+	auto access(cl::sycl::handler& h){
+		return data_.template get_access<AccessMode>(h);
+	}
+
+	template<cl::sycl::access::mode AccessMode>
+	auto access(cl::sycl::handler& h) const {
+		return data_.template get_access<AccessMode>(h);
+	}
+};
+}
+
 namespace kel {
 namespace lbm {
 namespace sch {
@@ -36,10 +112,13 @@ using CellInfo = Cell<UInt8, Desc, 1u, 0u, 0u>;
  * Basic type for simulation
  */
 template<typename Desc>
-using CellStruct = Struct<
-	Member<DfCell<Desc>, "dfs">,
-	Member<DfCell<Desc>, "dfs_old">,
-	Member<CellInfo<Desc>, "info">
+using CellStruct = CellFieldStruct<
+	Desc,
+	Struct<
+		Member<CellField<Desc,DfCell<Desc>>, "dfs">,
+		Member<CellField<Desc,DfCell<Desc>>, "dfs_old">,
+		Member<CellField<Desc,CellInfo<Desc>>, "info">
+	>
 >;
 
 template<typename T, uint64_t D>
@@ -48,7 +127,7 @@ using MacroStruct = Struct<
 	Member<T, "pressure">
 >;
 
-using CavityFieldD2Q9 = CellField<D2Q9, CellStruct<D2Q9>>;
+using CavityFieldD2Q9 = CellStruct<D2Q9>;
 }
 
 
@@ -104,22 +183,20 @@ public:
 		*/
 	}
 };
-}
-}
 
 constexpr size_t dim_size = 2;
 constexpr size_t dim_x = 128;
 constexpr size_t dim_y = 128;
 
-void set_geometry(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
+void set_geometry(saw::data<sch::CavityFieldD2Q9>& lattice){
 	using namespace kel::lbm;
-	auto meta = latt.meta();
+	auto meta = lattice.meta();
+	auto& info_field = lattice.template get<"info">();
 	/**
 	 * Set ghost
 	 */
 	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
+		auto& info = info_field.at(index);
 
 		info({0u}).set(0u);
 
@@ -129,8 +206,7 @@ void set_geometry(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 	 * Set wall
 	 */
 	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
+		auto& info = info_field.at(index);
 
 		info({0u}).set(2u);
 
@@ -140,8 +216,7 @@ void set_geometry(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 	 * Set fluid
 	 */
 	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
+		auto& info = info_field.at(index);
 
 		info({0u}).set(1u);
 
@@ -151,30 +226,32 @@ void set_geometry(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 	 * Set top lid
 	 */
 	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
+		auto& info = info_field.at(index);
 
 		info({0u}).set(3u);
 
 	}, {{0u,1u}}, {{meta.at({0u}), 2u}}, {{2u,0u}});
 }
 
-void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
+void set_initial_conditions(
+	saw::data<sch::CavityFieldD2Q9>& lattice
+){
 	using namespace kel::lbm;
 
 	saw::data<sch::T> rho{1.0};
 	saw::data<sch::FixedArray<sch::T,sch::D2Q9::D>> vel{{0.0,0.0}};
 	auto eq = equilibrium<sch::T,sch::D2Q9>(rho, vel);
 
-	auto meta = latt.meta();
+	auto meta = lattice.meta();
+	auto& dfs_field = lattice.template get<"dfs">();
+	auto& dfs_old_field = lattice.template get<"dfs_old">();
 
 	/**
 	 * Set distribution
 	 */
 	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& dfs = cell.template get<"dfs">();
-		auto& dfs_old = cell.template get<"dfs_old">();
+		auto& dfs = dfs_field.at(index);
+		auto& dfs_old = dfs_old_field.at(index);
 
 		for(saw::data<sch::UInt64> k = 0; k < saw::data<sch::UInt64>{sch::D2Q9::Q}; ++k){
 			dfs(k) = eq.at(k);
@@ -184,84 +261,150 @@ void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 	}, {{0u,0u}}, meta);
 }
 
+struct sycl_component_context {
+	acpp::sycl::handler* handler = nullptr;
+};
+
 void lbm_step(
-	acpp::sycl::buffer<saw::data<sch::CavityFieldD2Q9>>& latt,
+	saw::data<sch::CellField<sch::D2Q9,sch::CellInfo<sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>>& info,
+	saw::data<sch::CellField<sch::D2Q9,kel::lbm::sch::DfCell  <sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>>& dfs,
+	saw::data<sch::CellField<sch::D2Q9,kel::lbm::sch::DfCell  <sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>>& dfs_old,
 	bool even_step,
 	uint64_t time_step,
 	acpp::sycl::queue& sycl_q
 ){
-	using namespace kel::lbm;
-	using dfi = df_info<sch::T,sch::D2Q9>;
+  using namespace kel::lbm;
+  using namespace acpp;
+
+  using dfi = df_info<sch::T, sch::D2Q9>;
+
+  //component<sch::T, sch::D2Q9, cmpt::BGK> coll{0.59};
+  //component<sch::T, sch::D2Q9, cmpt::BounceBack> bb;
+  component<sch::T, sch::D2Q9, cmpt::MovingWall> bb_lid;
+  bb_lid.lid_vel = {0.1, 0.0};
+
+  const size_t Nx = latt.template get_dim_size<0>().get();
+  const size_t Ny = latt.template get_dim_size<1>().get();
+
+  // Submit collision kernel
+  sycl_q.submit([&](sycl::handler& cgh) {
+    // Accessor for latt with read/write
+    auto info_acc = info.template get_access<sycl::access::mode::read>(cgh);
+    auto dfs_acc = dfs.template get_access<sycl::access::mode::read_write>(cgh);
+    auto dfs_old_acc = dfs_old.template get_access<sycl::access::mode::read_write>(cgh);
+
+    cgh.parallel_for(
+      sycl::range<2>{Nx, Ny}, [=](sycl::id<2> idx) {
+        size_t i = idx[0];
+        size_t j = idx[1];
+
+        // Read cell info
+        auto& info = info_acc[idx];
+
+        switch (info({0u}).get()) {
+          case 1u: {
+            // bb.apply(latt_acc, {i, j}, time_step);
+				// bool is_even = ((time_step % 2) == 0);
+
+				auto& dfs_old = (is_even) ? dfs_old_acc[idx] : dfs_acc[idx];
+				// 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()));
+				}
+            break;
+          }
+          case 2u: {
+            // coll.apply(latt_acc, {i, j}, time_step);
+				// bool is_even = ((time_step % 2) == 0);
 
-  /**
-	* 1. Relaxation parameter \tau
-	*/
-	component<sch::T, sch::D2Q9, cmpt::BGK> coll{0.59};
-	component<sch::T, sch::D2Q9, cmpt::BounceBack> bb;
-	component<sch::T, sch::D2Q9, cmpt::MovingWall> bb_lid;
-	bb_lid.lid_vel = {0.1,0.0};
+				auto& dfs_old = (is_even) ? dfs_old_acc[idx] : dfs_acc[idx];
+				// auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
 
-	auto meta = latt.meta();
+				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);
 
-	/**
-	 * Collision
-	 */
-	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
-
-		auto& dfs = cell.template get<"dfs">();
-		auto& dfs_old = cell.template get<"dfs_old">();
-
-		switch(info({0u}).get()){
-			case 1u: {
-				bb.apply(latt, index, time_step);
-				break;
-			}
-			case 2u: {
-				coll.apply(latt, index, time_step);
-				break;
-			}
-			default:
-				break;
-		}
+				using dfi = df_info<T,Descriptor>;
 
-	}, {{0u,0u}}, meta);
+				auto& force = cell.template get<"force">();
 
-	// Stream
-	for(uint64_t i = 1u; (i+1u) < latt.template get_dim_size<0>().get(); ++i){
-		for(uint64_t j = 1u; (j+1u) < latt.template get_dim_size<1>().get(); ++j){
-			auto& cell = latt({{i,j}});
-			auto& df_new = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">();
-			auto& info_new = cell.template get<"info">();
-
-			if(info_new({0u}).get() > 0u && info_new({0u}).get() != 3u){
-				for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){
-					auto dir = dfi::directions[dfi::opposite_index[k]];
-					auto& cell_dir_old = latt({{i+dir[0],j+dir[1]}});
-
-					auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">();
-					auto& info_old = cell_dir_old.template get<"info">();
-
-					if( info_old({0}).get() == 3u ){
-						auto& df_old_loc = even_step ? latt({{i,j}}).template get<"dfs_old">() : latt({{i,j}}).template get<"dfs">();
-						df_new({k}) = df_old_loc({dfi::opposite_index.at(k)}) - 2.0 * dfi::inv_cs2 * dfi::weights.at(k) * 1.0 * ( bb_lid.lid_vel[0] * dir[0] + bb_lid.lid_vel[1] * dir[1]);
-						// dfs({dfi::opposite_index.at(i)}) = dfs_cpy({i}) - 2.0 * dfi::weights[i] * 1.0 * ( lid_vel[0] * dfi::directions[i][0] + lid_vel[1] * dfi::directions[i][1]) * dfi::inv_cs2;
-					} else {
-						df_new({k}) = df_old({k});
+				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_);
 				}
-			}
-		}
-	}
+            break;
+          }
+          default:
+            // Do nothing
+            break;
+        }
+      });
+  });
+
+  // Submit streaming kernel
+  sycl_q.submit([&](sycl::handler& cgh) {
+    auto info_acc = info.template get_access<sycl::access::mode::read_write>(cgh);
+
+    cgh.parallel_for(
+      sycl::range<2>{Nx - 2, Ny - 2}, [=](sycl::id<2> idx) {
+        size_t i = idx[0] + 1; // avoid boundary
+        size_t j = idx[1] + 1;
+
+        auto& df_new = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">();
+        auto& info_new = cell.template get<"info">();
+
+        if (info_new({0u}).get() > 0u && info_new({0u}).get() != 3u) {
+          for (uint64_t k = 0u; k < sch::D2Q9::Q; ++k) {
+            auto dir = dfi::directions[dfi::opposite_index[k]];
+            auto& cell_dir_old = latt_acc({i + dir[0], j + dir[1]});
+
+            auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">();
+            auto& info_old = cell_dir_old.template get<"info">();
+
+            if (info_old({0}).get() == 3u) {
+              auto& df_old_loc = even_step ? latt_acc({i, j}).template get<"dfs_old">() : latt_acc({i, j}).template get<"dfs">();
+              df_new({k}) = df_old_loc({dfi::opposite_index.at(k)}) - 2.0 * dfi::inv_cs2 * dfi::weights.at(k) * 1.0 * (bb_lid.lid_vel[0] * dir[0] + bb_lid.lid_vel[1] * dir[1]);
+            } else {
+              df_new({k}) = df_old({k});
+            }
+          }
+        }
+      });
+  });
+
+  sycl_q.wait();
+}
+}
 }
 
 int main(){
 	using namespace kel::lbm;
 
-	saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>> dim{{dim_x, dim_y}};
+	saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>> dim{{dim_x,info_field.at dim_y}};
+	const dim_size = dim_x * dim_y;
 
-	saw::data<sch::CavityFieldD2Q9, saw::encode::Native> lattice{dim};
+	saw::data<sch::CavityFieldD2Q9>> lattice{dim};
 
 	converter<sch::T> conv{
 		{0.1},
@@ -287,13 +430,18 @@ int main(){
 	 */
 	set_initial_conditions(lattice);
 
+	auto& info_field = lattice.template get<"info">();
+	auto& dfs_field = lattice.template get<"dfs">();
+	auto& dfs_old_field = lattice.template get<"dfs_old">();
+
 	acpp::sycl::queue sycl_q{acpp::sycl::default_selector_v, acpp::sycl::property::queue::in_order{}};
-	acpp::sycl::buffer<saw::data<sch::CavityFieldD2Q9>> df_sycl{&lattice, 1u};
+	saw::data<sch::CellField<Desc,sch::CellInfo<sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>> sycl_info_field{info_field};
+	saw::data<sch::CellField<Desc,sch::DfCell<sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>> sycl_dfs_field{dfs_field};
+	saw::data<sch::CellField<Desc,sch::DfCell<sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>> sycl_dfs_old_field{dfs_old_field};
 
 	/**
 	 * Timeloop
 	 */
-
 	uint64_t lattice_steps = 512000u;
 	bool even_step = true;
 
@@ -310,7 +458,7 @@ int main(){
 			write_vtk_file(vtk_f_name, macros);
 		}
 
-		lbm_step(df_sycl, (i%2u == 0u), i, sycl_q);
+		lbm_step(sycl_info_field, sycl_dfs_field, sycl_dfs_old_field, (i%2u == 0u), i, sycl_q);
 	}
 	return 0;
 }
diff --git a/examples/config.json b/examples/config.json
index 7c6d3a2..ba5d447 100644
--- a/examples/config.json
+++ b/examples/config.json
@@ -1,5 +1,5 @@
 {
-	"delta_x" : 1.0,
-	"delta_t" : 1.0,
-	"kinematic_viscosity" : 3e-2
+	"delta_x" : 10.0,
+	"delta_t" : 1e5,
+	"kinematic_viscosity" : 1.5e-5
 }
diff --git a/examples/meta_2d/.nix/derivation.nix b/examples/meta_2d/.nix/derivation.nix
new file mode 100644
index 0000000..6dcc00e
--- /dev/null
+++ b/examples/meta_2d/.nix/derivation.nix
@@ -0,0 +1,33 @@
+{ lib
+, stdenv
+, scons
+, clang-tools
+, forstio
+, pname
+, version
+, kel-lbm
+}:
+
+stdenv.mkDerivation {
+	pname = pname + "-examples-" + "meta_2d";
+	inherit version;
+  src = ./..;
+
+  nativeBuildInputs = [
+    scons
+    clang-tools
+  ];
+
+  buildInputs = [
+    forstio.core
+    forstio.async
+    forstio.codec
+    forstio.codec-unit
+		forstio.codec-json
+		kel-lbm
+  ];
+
+	preferLocalBuild = true;
+
+  outputs = [ "out" "dev" ];
+}
diff --git a/examples/meta_2d/SConscript b/examples/meta_2d/SConscript
new file mode 100644
index 0000000..f14fd77
--- /dev/null
+++ b/examples/meta_2d/SConscript
@@ -0,0 +1,32 @@
+#!/bin/false
+
+import os
+import os.path
+import glob
+
+
+Import('env')
+
+dir_path = Dir('.').abspath
+
+# Environment for base library
+examples_env = env.Clone();
+
+examples_env.sources = sorted(glob.glob(dir_path + "/*.cpp"))
+examples_env.headers = sorted(glob.glob(dir_path + "/*.hpp"))
+
+env.sources += examples_env.sources;
+env.headers += examples_env.headers;
+
+# Cavity2D
+examples_objects = [];
+examples_env.add_source_files(examples_objects, ['meta_2d.cpp'], shared=False);
+examples_env.meta_2d = examples_env.Program('#bin/meta_2d', [examples_objects]);
+
+# Set Alias
+env.examples = [
+	examples_env.meta_2d
+];
+env.Alias('examples', env.examples);
+env.targets += ['examples'];
+env.Install('$prefix/bin/', env.examples);
diff --git a/SConstruct b/examples/meta_2d/SConstruct
index 72f7c7f..a7201cb 100644
--- a/SConstruct
+++ b/examples/meta_2d/SConstruct
@@ -57,8 +57,7 @@ env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[],
         '-Wextra'
     ],
     LIBS=[
-        'forstio-core',
-        'forstio-codec'
+        'forstio-core'
     ]
 );
 env.__class__.add_source_files = add_kel_source_files
@@ -71,9 +70,7 @@ env.headers = [];
 env.targets = [];
 
 Export('env')
-SConscript('c++/SConscript')
-SConscript('examples/SConscript')
-SConscript('tests/SConscript')
+SConscript('SConscript')
 
 env.Alias('cdb', env.cdb);
 env.Alias('all', [env.targets]);
diff --git a/examples/meta_2d.cpp b/examples/meta_2d/meta_2d.cpp
index ac0857b..6d19800 100644
--- a/examples/meta_2d.cpp
+++ b/examples/meta_2d/meta_2d.cpp
@@ -1,4 +1,4 @@
-#include "../c++/lbm.hpp"
+#include <kel/lbm/lbm.hpp>
 
 #include <iostream>
 
diff --git a/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp
index f3caf3a..e4c6de6 100644
--- a/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp
+++ b/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp
@@ -516,8 +516,9 @@ void couple_particles_to_lattice(
 		auto v_n = saw::math::dot(solid_normal, (p_vel + p_acc));
 
 		if(v_n.at({}).get() < 0.0){
-			solid_normal.at({{0u}}) = solid_normal.at({{0u}})*v_n.at({});
-			solid_normal.at({{1u}}) = solid_normal.at({{1u}})*v_n.at({});
+			auto v_n_2 = v_n + v_n;
+			solid_normal.at({{0u}}) = solid_normal.at({{0u}})*(v_n_2.at({}));
+			solid_normal.at({{1u}}) = solid_normal.at({{1u}})*(v_n_2.at({}));
 			p_acc = p_acc - solid_normal;
 		}
 
diff --git a/lib/SConstruct b/lib/SConstruct
new file mode 100644
index 0000000..8b3ab01
--- /dev/null
+++ b/lib/SConstruct
@@ -0,0 +1,75 @@
+#!/usr/bin/env python3
+
+import sys
+import os
+import os.path
+import glob
+import re
+
+
+if sys.version_info < (3,):
+    def isbasestring(s):
+        return isinstance(s,basestring)
+else:
+    def isbasestring(s):
+        return isinstance(s, (str,bytes))
+
+def add_kel_source_files(self, sources, filetype, lib_env=None, shared=False, target_post=""):
+
+    if isbasestring(filetype):
+        dir_path = self.Dir('.').abspath
+        filetype = sorted(glob.glob(dir_path+"/"+filetype))
+
+    for path in filetype:
+        target_name = re.sub( r'(.*?)(\.cpp|\.c\+\+)', r'\1' + target_post, path )
+        if shared:
+            target_name+='.os'
+            sources.append( self.SharedObject( target=target_name, source=path ) )
+        else:
+            target_name+='.o'
+            sources.append( self.StaticObject( target=target_name, source=path ) )
+    pass
+
+def isAbsolutePath(key, dirname, env):
+	assert os.path.isabs(dirname), "%r must have absolute path syntax" % (key,)
+
+env_vars = Variables(
+	args=ARGUMENTS
+)
+
+env_vars.Add('prefix',
+	help='Installation target location of build results and headers',
+	default='/usr/local/',
+	validator=isAbsolutePath
+)
+
+env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[],
+    CPPDEFINES=['SAW_UNIX'],
+    CXXFLAGS=[
+        '-std=c++20',
+        '-g',
+        '-Wall',
+        '-Wextra'
+    ],
+    LIBS=[
+        'forstio-core'
+    ]
+);
+env.__class__.add_source_files = add_kel_source_files
+env.Tool('compilation_db');
+env.cdb = env.CompilationDatabase('compile_commands.json');
+
+env.objects = [];
+env.sources = [];
+env.headers = [];
+env.targets = [];
+
+Export('env')
+SConscript('c++/SConscript')
+SConscript('tests/SConscript')
+
+env.Alias('cdb', env.cdb);
+env.Alias('all', [env.targets]);
+env.Default('all');
+
+env.Alias('install', '$prefix')
diff --git a/c++/SConscript b/lib/c++/SConscript
index 85a078f..85a078f 100644
--- a/c++/SConscript
+++ b/lib/c++/SConscript
diff --git a/c++/args.hpp b/lib/c++/args.hpp
index 99c6172..99c6172 100644
--- a/c++/args.hpp
+++ b/lib/c++/args.hpp
diff --git a/c++/boundary.hpp b/lib/c++/boundary.hpp
index 4e30f71..5cc7705 100644
--- a/c++/boundary.hpp
+++ b/lib/c++/boundary.hpp
@@ -18,8 +18,8 @@ struct ZouHeVertical{};
 /**
  * Full-Way BounceBack.
  */
-template<typename T, typename Descriptor>
-class component<T, Descriptor, cmpt::BounceBack> {
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::BounceBack, Encode> {
 private:
 public:
 	component() = default;
@@ -44,7 +44,7 @@ public:
 	 * 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){
+	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
@@ -76,8 +76,8 @@ public:
  * 0 - 2 - 2
  *
  */
-template<typename FP, typename Descriptor, bool Dir>
-class component<FP, Descriptor, cmpt::ZouHeHorizontal<Dir>> final {
+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_;
@@ -88,7 +88,7 @@ public:
 	{}
 
 	template<typename CellFieldSchema>
-	void apply(saw::data<CellFieldSchema>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, uint64_t time_step){
+	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);
diff --git a/c++/collision.hpp b/lib/c++/collision.hpp
index 47870c1..9ab542b 100644
--- a/c++/collision.hpp
+++ b/lib/c++/collision.hpp
@@ -14,8 +14,8 @@ struct BGKGuo {};
 /**
  * Standard BGK collision operator for LBM
  */
-template<typename T, typename Descriptor>
-class component<T, Descriptor, cmpt::BGK> {
+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_;
@@ -45,7 +45,7 @@ public:
 	 * Raw setup
 	 */
 	template<typename CellFieldSchema>
-	void apply(saw::data<CellFieldSchema>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, uint64_t time_step){
+	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);
 
@@ -64,8 +64,8 @@ public:
 	}
 };
 
-template<typename T, typename Descriptor>
-class component<T, Descriptor, cmpt::BGKGuo> {
+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_;
@@ -90,7 +90,7 @@ public:
 	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){
+	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);
 
diff --git a/c++/component.hpp b/lib/c++/component.hpp
index 87ac31c..1e5dbbf 100644
--- a/c++/component.hpp
+++ b/lib/c++/component.hpp
@@ -32,7 +32,7 @@ public:
 /**
  * Component class which forms the basis of each processing step
  */
-template<typename T, typename Descriptor, typename Cmpt>
+template<typename T, typename Descriptor, typename Cmpt, typename Encode = saw::encode::Native>
 class component {};
 }
 }
diff --git a/c++/config.hpp b/lib/c++/config.hpp
index 64f7a0f..64f7a0f 100644
--- a/c++/config.hpp
+++ b/lib/c++/config.hpp
diff --git a/c++/converter.hpp b/lib/c++/converter.hpp
index 5c19c68..5c19c68 100644
--- a/c++/converter.hpp
+++ b/lib/c++/converter.hpp
diff --git a/c++/descriptor.hpp b/lib/c++/descriptor.hpp
index d04d217..51d5814 100644
--- a/c++/descriptor.hpp
+++ b/lib/c++/descriptor.hpp
@@ -25,7 +25,9 @@ struct Cell {
 };
 
 template<typename Desc, typename Cell>
-struct CellField;
+struct CellField{
+	using Descriptor = Desc;
+};
 
 template<typename Desc, typename... CellFieldTypes, saw::string_literal... CellFieldNames>
 struct CellField<
@@ -33,7 +35,16 @@ struct CellField<
 	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>;
+};
 
 }
 
@@ -292,5 +303,63 @@ public:
 	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/c++/environment.hpp b/lib/c++/environment.hpp
index 4915e3a..4915e3a 100644
--- a/c++/environment.hpp
+++ b/lib/c++/environment.hpp
diff --git a/c++/equilibrium.hpp b/lib/c++/equilibrium.hpp
index bb55d00..bb55d00 100644
--- a/c++/equilibrium.hpp
+++ b/lib/c++/equilibrium.hpp
diff --git a/c++/geometry.hpp b/lib/c++/geometry.hpp
index 9802feb..9802feb 100644
--- a/c++/geometry.hpp
+++ b/lib/c++/geometry.hpp
diff --git a/c++/hlbm.hpp b/lib/c++/hlbm.hpp
index 1c665ce..1c665ce 100644
--- a/c++/hlbm.hpp
+++ b/lib/c++/hlbm.hpp
diff --git a/c++/iterator.hpp b/lib/c++/iterator.hpp
index 78babff..78babff 100644
--- a/c++/iterator.hpp
+++ b/lib/c++/iterator.hpp
diff --git a/c++/lbm.hpp b/lib/c++/lbm.hpp
index 39c42af..39c42af 100644
--- a/c++/lbm.hpp
+++ b/lib/c++/lbm.hpp
diff --git a/c++/lbm_unit.hpp b/lib/c++/lbm_unit.hpp
index 2d90652..2d90652 100644
--- a/c++/lbm_unit.hpp
+++ b/lib/c++/lbm_unit.hpp
diff --git a/c++/macroscopic.hpp b/lib/c++/macroscopic.hpp
index 51368e9..51368e9 100644
--- a/c++/macroscopic.hpp
+++ b/lib/c++/macroscopic.hpp
diff --git a/c++/particle/geometry/circle.hpp b/lib/c++/particle/geometry/circle.hpp
index 331f06d..331f06d 100644
--- a/c++/particle/geometry/circle.hpp
+++ b/lib/c++/particle/geometry/circle.hpp
diff --git a/c++/particle/particle.hpp b/lib/c++/particle/particle.hpp
index 39aadfb..39aadfb 100644
--- a/c++/particle/particle.hpp
+++ b/lib/c++/particle/particle.hpp
diff --git a/c++/statistics.hpp b/lib/c++/statistics.hpp
index c07ccb7..c07ccb7 100644
--- a/c++/statistics.hpp
+++ b/lib/c++/statistics.hpp
diff --git a/c++/term_renderer.hpp b/lib/c++/term_renderer.hpp
index 5cbb551..5cbb551 100644
--- a/c++/term_renderer.hpp
+++ b/lib/c++/term_renderer.hpp
diff --git a/c++/util.hpp b/lib/c++/util.hpp
index 0bdebd1..0bdebd1 100644
--- a/c++/util.hpp
+++ b/lib/c++/util.hpp
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/c++/write_vtk.hpp b/lib/c++/write_vtk.hpp
index 0647db5..0647db5 100644
--- a/c++/write_vtk.hpp
+++ b/lib/c++/write_vtk.hpp
diff --git a/tests/SConscript b/lib/tests/SConscript
index d1b381e..d1b381e 100644
--- a/tests/SConscript
+++ b/lib/tests/SConscript
diff --git a/tests/converter.cpp b/lib/tests/converter.cpp
index 4fc536f..4fc536f 100644
--- a/tests/converter.cpp
+++ b/lib/tests/converter.cpp
diff --git a/tests/descriptor.cpp b/lib/tests/descriptor.cpp
index a8337e6..a8337e6 100644
--- a/tests/descriptor.cpp
+++ b/lib/tests/descriptor.cpp
diff --git a/tests/equilibrium.cpp b/lib/tests/equilibrium.cpp
index 9201e55..9201e55 100644
--- a/tests/equilibrium.cpp
+++ b/lib/tests/equilibrium.cpp
diff --git a/tests/iterator.cpp b/lib/tests/iterator.cpp
index cd1cb7c..cd1cb7c 100644
--- a/tests/iterator.cpp
+++ b/lib/tests/iterator.cpp
diff --git a/tests/particle_flow_coupling.cpp b/lib/tests/particle_flow_coupling.cpp
index c3e3769..c3e3769 100644
--- a/tests/particle_flow_coupling.cpp
+++ b/lib/tests/particle_flow_coupling.cpp
diff --git a/tests/particles.cpp b/lib/tests/particles.cpp
index 277a8d0..277a8d0 100644
--- a/tests/particles.cpp
+++ b/lib/tests/particles.cpp
diff --git a/tests/vtk_write.cpp b/lib/tests/vtk_write.cpp
index 0df9998..0df9998 100644
--- a/tests/vtk_write.cpp
+++ b/lib/tests/vtk_write.cpp