Writing was flipped on non square, cubic domains

7 files changed, 543 insertions, 21 deletions

diff --git a/c++/descriptor.hpp b/c++/descriptor.hpp
index 23c82f2..0d05a04 100644
--- a/c++/descriptor.hpp
+++ b/c++/descriptor.hpp
@@ -202,6 +202,7 @@ public:
 	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);
 	}
diff --git a/c++/geometry.hpp b/c++/geometry.hpp
index fe0fe7e..9802feb 100644
--- a/c++/geometry.hpp
+++ b/c++/geometry.hpp
@@ -2,11 +2,13 @@
 
 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++/iterator.hpp b/c++/iterator.hpp
index fcc50bc..78babff 100644
--- a/c++/iterator.hpp
+++ b/c++/iterator.hpp
@@ -14,5 +14,19 @@ void iterate_over(Func&& func, const saw::data<sch::FixedArray<sch::UInt64,2u>>&
 	}
 	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
index d331c6a..6bcd1e7 100644
--- a/c++/lbm.hpp
+++ b/c++/lbm.hpp
@@ -5,6 +5,7 @@
 #include "config.hpp"
 #include "component.hpp"
 #include "equilibrium.hpp"
+#include "macroscopic.hpp"
 #include "write_vtk.hpp"
 
 #include <forstio/codec/unit/unit_print.hpp>
diff --git a/c++/write_vtk.hpp b/c++/write_vtk.hpp
index 5cbc6c0..40597fd 100644
--- a/c++/write_vtk.hpp
+++ b/c++/write_vtk.hpp
@@ -62,13 +62,13 @@ struct lbm_vtk_writer<sch::Array<sch::Struct<sch::Member<StructT,StructN>...> ,
 		constexpr auto Lit = saw::parameter_key_pack_type<i, StructN...>::literal;
 		using Type = typename saw::parameter_pack_type<i,StructT...>::type;
 
-		if constexpr (Dep == Dim){
+		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 < Dim, "Don't fall into this case");
-			for(index.at({Dep}) = 0; index.at({Dep}) < field.get_dims().at({Dep}); ++index.at({Dep})){
-				auto eov = write_i_iterate_d<i,Dep+1>(vtk_file, field, index);
+			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;
 				}
@@ -95,7 +95,7 @@ struct lbm_vtk_writer<sch::Array<sch::Struct<sch::Member<StructT,StructN>...> ,
 			}
 
 		}
-		return write_i_iterate_d<i,0u>(vtk_file, field, index);
+		return write_i_iterate_d<i,Dim>(vtk_file, field, index);
 	}
 
 	template<uint64_t i>
diff --git a/examples/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu.cpp
new file mode 100644
index 0000000..e1e6333
--- /dev/null
+++ b/examples/cavity_2d_gpu.cpp
@@ -0,0 +1,314 @@
+#include "../c++/descriptor.hpp"
+#include "../c++/macroscopic.hpp"
+#include "../c++/lbm.hpp"
+#include "../c++/component.hpp"
+#include "../c++/collision.hpp"
+#include "../c++/boundary.hpp"
+
+#include <forstio/codec/data.hpp>
+// #include <forstio/remote/
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+
+namespace kel {
+namespace lbm {
+namespace sch {
+using namespace saw::schema;
+
+/**
+ * Basic distribution function
+ * Base type
+ * D
+ * Q
+ * Scalar factor
+ * D factor
+ * Q factor
+ */
+using T = Float32;
+using D2Q5 = Descriptor<2u,5u>;
+using D2Q9 = Descriptor<2u,9u>;
+
+template<typename Desc>
+using DfCell = Cell<T, Desc, 0u, 0u, 1u>;
+
+template<typename Desc>
+using CellInfo = Cell<UInt8, D2Q9, 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">
+>;
+
+template<typename T, uint64_t D>
+using MacroStruct = Struct<
+	Member<FixedArray<T,D>, "velocity">,
+	Member<T, "pressure">
+>;
+
+using CavityFieldD2Q9 = CellField<D2Q9, CellStruct<D2Q9>>;
+}
+
+
+/*
+template<typename T, typename Encode>
+class df_cell_view;
+*/
+/**
+ * Minor helper for the AA-Pull Pattern, so I can use only one lattice
+ *
+ * Am I sure I want to use AA this way?
+ * Esoteric Twist technically reduces the needed memory access footprint
+ */
+/*
+template<typename Desc, size_t SN, size_t DN, size_t QN, typename Encode>
+class df_cell_view<sch::Cell<sch::T, Desc, SN, DN, QN>, Encode> {
+public:
+	using Schema = sch::Cell<sch::T,Desc,SN,DN,QN>;
+private:
+		std::array<std::decay_t<typename saw::native_data_type<sch::T>::type>*, QN> view_;
+public:
+		df_cell_view(const std::array<std::decay_t<typename saw::native_data_type<sch::T>::type>*, QN>& view):
+				view_{view}
+		{}
+};
+*/
+namespace cmpt {
+struct MovingWall {};
+}
+
+/**
+ * Full-Way moving wall Bounce back, something is not right here.
+ * Technically it should reflect properly.
+ */
+template<typename Desc>
+class component<sch::T, Desc, cmpt::MovingWall> {
+public:
+	std::array<typename saw::native_data_type<sch::T>::type, Desc::D> lid_vel;
+
+public:
+	void apply(
+		saw::data<sch::DfCell<Desc>>& dfs
+	){
+		using dfi = df_info<sch::T,Desc>;
+
+		// Technically use .copy()
+		/*
+		auto dfs_cpy = dfs;
+
+		for(uint64_t i = 0u; i < Desc::Q; ++i){
+			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;
+		}
+		*/
+	}
+};
+}
+}
+
+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){
+	using namespace kel::lbm;
+	/**
+	 * Set ghost
+	 */
+	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+		auto& cell = latt(index);
+		auto& info = cell.template get<"info">();
+
+		info({0u}).set(0u);
+
+	}, {{0u,0u}}, meta);
+
+	/**
+	 * Set wall
+	 */
+	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+		auto& cell = latt(index);
+		auto& info = cell.template get<"info">();
+
+		info({0u}).set(2u);
+
+	}, {{0u,0u}}, meta, {{1u,1u}});
+
+	/**
+	 * Set fluid
+	 */
+	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+		auto& cell = latt(index);
+		auto& info = cell.template get<"info">();
+
+		info({0u}).set(1u);
+
+	}, {{0u,0u}}, meta, {{2u,2u}});
+
+	/**
+	 * Set top lid
+	 */
+	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+		auto& cell = latt(index);
+		auto& info = cell.template get<"info">();
+
+		info({0u}).set(3u);
+
+	}, {{0u,1u}}, {{meta.at({0u}), 2u}}, {{2u,0u}});
+}
+
+void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
+	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();
+
+	/**
+	 * 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">();
+
+		for(saw::data<sch::UInt64> k = 0; k < saw::data<sch::UInt64>{sch::D2Q9::Q}; ++k){
+			dfs(k) = eq.at(k);
+			dfs_old(k) = eq.at(k);
+		}
+
+	}, {{0u,0u}}, meta);
+}
+
+void lbm_step(
+	saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt,
+	uint64_t time_step
+){
+	using namespace kel::lbm;
+	using dfi = df_info<sch::T,sch::D2Q9>;
+
+  /**
+	* 1. Relaxation parameter \tau
+	*/
+	component<sch::T, sch::D2Q9, cmpt::BGK> coll{0.5384};
+	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 meta = latt.meta();
+
+	/**
+	 * 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: {
+				coll.apply(latt, index, time_step);
+				break;
+			}
+			case 2u: {
+				bb.apply(latt, index, time_step);
+				break;
+			}
+			default:
+				break;
+		}
+
+	}, {{0u,0u}}, meta);
+
+	// 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});
+					}
+				}
+			}
+		}
+	}
+}
+
+int main(){
+	using namespace kel::lbm;
+
+	saw::remote<saw::rmt::Sycl> sycl_rmt;
+
+	saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>> dim{{dim_x, dim_y}};
+
+	saw::data<sch::CavityFieldD2Q9, saw::encode::Native> lattice{dim};
+
+	converter<sch::T> conv{
+		{0.1},
+		{0.1}
+	};
+
+	print_lbm_meta<sch::T, sch::D2Q9>(conv, {1e-3});
+
+	//auto& df_field = lattices.at(0).template get<"dfs">();
+	//for(uint64_t i = 0; i < df_field.get_dim_size<0u>(); ++i){
+	//	lattices.at(i) = {dim_x, dim_y};
+	//}
+
+	/**
+	 * Set meta information describing what this cell is
+	 */
+	set_geometry(lattice);
+
+	/**
+	 *
+	 */
+	set_initial_conditions(lattice);
+
+	/**
+	 * Timeloop
+	 */
+
+	uint64_t lattice_steps = 512000u;
+	bool even_step = true;
+
+	uint64_t print_every = 256u;
+	uint64_t file_no = 0u;
+
+	saw::data<sch::Array<sch::MacroStruct<sch::T,sch::D2Q9::D>,sch::D2Q9::D>> macros{dim};
+
+	for(uint64_t i = 0u; i < 256u; ++i){
+		{
+			std::string vtk_f_name{"tmp/poiseulle_2d_"};
+			vtk_f_name += std::to_string(i) + ".vtk";
+			write_vtk_file(vtk_f_name, macros);
+		}
+
+		// lbm_step(lattice, i);
+	}
+	return 0;
+}
diff --git a/examples/poiseulle_2d.cpp b/examples/poiseulle_2d.cpp
index dcdc7c3..54623f7 100644
--- a/examples/poiseulle_2d.cpp
+++ b/examples/poiseulle_2d.cpp
@@ -45,8 +45,92 @@ using MacroStruct = Struct<
 	Member<T, "pressure">
 >;
 
+template<typename T>
+using GeometryStruct = Struct<
+	Member<T, "info">
+>;
+
 using CavityFieldD2Q9 = CellField<D2Q9, CellStruct<D2Q9>>;
 }
+
+namespace cmpt {
+template<bool East>
+struct PressureBoundaryRestrictedVelocityTo {};
+}
+
+/**
+ * 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 East>
+struct component<FP,Descriptor, cmpt::PressureBoundaryRestrictedVelocityTo<East>> {
+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__}
+	{}
+
+	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& 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<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});
+			}
+		}
+		constexpr int known_dir = East ? 1 : 1;
+		auto sum_unknown_dfs = (pressure_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 / pressure_setting_;
+
+		if constexpr (East) {
+			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 East){
+			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}));
+		}
+	}
+};
 }
 }
 
@@ -114,8 +198,96 @@ void set_geometry(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 	using namespace kel::lbm;
 
+	saw::data<sch::T> rho{1.0};
+	saw::data<sch::FixedArray<sch::T,sch::D2Q9::D>> vel{{0.01,0.0}};
+	auto eq = equilibrium<sch::T,sch::D2Q9>(rho, vel);
+
 	auto meta = latt.meta();
 
+	/**
+	 * 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">();
+
+		for(saw::data<sch::UInt64> k = 0; k < saw::data<sch::UInt64>{sch::D2Q9::Q}; ++k){
+			dfs(k) = eq.at(k);
+			dfs_old(k) = eq.at(k);
+		}
+
+	}, {{0u,0u}}, meta);
+}
+
+void lbm_step(
+	saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt,
+	uint64_t time_step
+){
+	using namespace kel::lbm;
+	using dfi = df_info<sch::T,sch::D2Q9>;
+
+	bool even_step = ((time_step % 2u) == 0u);
+  /**
+	* 1. Relaxation parameter \tau
+	*/
+	component<sch::T, sch::D2Q9, cmpt::BGK> coll{0.5384};
+	component<sch::T, sch::D2Q9, cmpt::BounceBack> bb;
+	component<sch::T, sch::D2Q9, cmpt::PressureBoundaryRestrictedVelocityTo<true>> inlet{1.1};
+	component<sch::T, sch::D2Q9, cmpt::PressureBoundaryRestrictedVelocityTo<false>> outlet{0.9};
+
+	auto meta = latt.meta();
+
+	/**
+	 * 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: {
+				coll.apply(latt, index, time_step);
+				break;
+			}
+			case 2u: {
+				bb.apply(latt, index, time_step);
+				break;
+			}
+			case 3u: {
+				inlet.apply(latt, index, time_step);
+				break;
+			}
+			case 4u: {
+				outlet.apply(latt, index, time_step);
+				break;
+			}
+			default:
+				break;
+		}
+
+	}, {{0u,0u}}, meta);
+
+	// Stream
+	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+		auto& cell = latt(index);
+		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){
+			for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){
+				auto dir = dfi::directions[dfi::opposite_index[k]];
+				auto& cell_dir_old = latt({{index.at({0u})+dir[0],index.at({1u})+dir[1]}});
+
+				auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">();
+
+				df_new({k}) = df_old({k});
+			}
+		}
+	}, {{0u,0u}}, meta);
 }
 
 int main(int argc, char** argv){
@@ -148,39 +320,57 @@ int main(int argc, char** argv){
 
 	print_lbm_meta<sch::Float64,sch::Descriptor<2u,9u>>(conv, {conf.template get<"kinematic_viscosity">()});
 
-	saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>> dim{{32u, 8u}};
+	saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>> dim{{1024u, 128u}};
 	saw::data<sch::CavityFieldD2Q9, saw::encode::Native> lattice{dim};
-
-	set_geometry(lattice);
+	auto meta = lattice.meta();
 
 	/**
-	 * Print basic setup info
+	 * Setup geometry
 	 */
-	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = lattice(index);
-		auto& info = cell.template get<"info">();
-		auto info_v = info({0u}).get();
-		std::cout<<(static_cast<uint32_t>(info_v));
+	set_geometry(lattice);
 
-		if( (index.at({1u}).get()+1u) < dim.at({1u}).get()){
-			std::cout<<" ";
-		}else{
-			std::cout<<"\n";
-		}
+	{
+
+		saw::data<sch::Array<sch::GeometryStruct<sch::T>, sch::D2Q9::D>> geo{dim};
+
+		iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+			auto& cell = lattice(index);
+			auto& info = cell.template get<"info">();
 
-	}, {{0u,0u}}, dim);
+			geo(index).template get<"info">().set(info({0u}).get());
+		}, {{0u,0u}}, dim);
+
+		std::string vtk_f_name{"tmp/geometry.vtk"};
+		write_vtk_file(vtk_f_name, geo);
+	}
+
+	/**
+	 * Setup DFs
+	 */
+	set_initial_conditions(lattice);
 
 	saw::data<sch::Array<sch::MacroStruct<sch::T,sch::D2Q9::D>,sch::D2Q9::D>> macros{dim};
 
 	for(uint64_t i = 0u; i < 256u; ++i){
+		bool even_step = ((i % 2u) == 0u);
 
 		{
+			// Stream
+			iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+				auto& cell = lattice(index);
+				auto& dfs = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">();
+
+				auto& rho = macros.at(index).template get<"pressure">();
+				auto& vel = macros.at(index).template get<"velocity">();
+				compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel);
+
+			}, {{0u,0u}}, meta);
 			std::string vtk_f_name{"tmp/poiseulle_2d_"};
 			vtk_f_name += std::to_string(i) + ".vtk";
 			write_vtk_file(vtk_f_name, macros);
 		}
 
-		// lbm_step(lattice, i);
+		lbm_step(lattice, i);
 	}
 
 	return 0;