Simulation works partially. Reconfirm with writing data out

5 files changed, 61 insertions, 69 deletions

diff --git a/lib/core/c++/boundary.hpp b/lib/core/c++/boundary.hpp
index 8ab2457..ec16edf 100644
--- a/lib/core/c++/boundary.hpp
+++ b/lib/core/c++/boundary.hpp
@@ -44,47 +44,22 @@ public:
 	 * Raw setup
 	 */
 	template<typename CellFieldSchema>
-	void apply(saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step){
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step)const{
 		bool is_even = ((time_step.get() % 2u) == 0u);
 
 		// This is a ref
-		auto& 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">();
+		// auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">();
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
 
 		using dfi = df_info<T,Descriptor>;
 
 		// Technically use .copy()
-		auto df_cpy = dfs_old.copy();
+		auto& dfs_old = dfs_old_f.at(index);
+		auto df_cpy = dfs_old;
 
-		for(uint64_t i = 1u; i < Descriptor::Q; ++i){
-			dfs_old({i}) = df_cpy({dfi::opposite_index.at(i)});
-		}
-	}
-	
-	template<typename CellStructSchema>
-	void apply(
-		saw::data<CellStructSchema, Encode>* field, 
-		const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& meta,
-		const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index,
-		saw::data<sch::UInt64> time_step
-	){
-		bool is_even = ((time_step.get() % 2u) == 0u);
-
-		// This is a ref
-		auto& cell = field[0u];
-
-		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)});
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			dfs_old.at({i}) = df_cpy.at({dfi::opposite_index.at(i)});
 		}
 	}
 };
@@ -104,59 +79,46 @@ public:
 template<typename FP, typename Descriptor, bool Dir, typename Encode>
 class component<FP, Descriptor, cmpt::ZouHeHorizontal<Dir>, Encode> final {
 private:
-	saw::data<FP> pressure_setting_;
 	saw::data<FP> rho_setting_;
 public:
-	component(const saw::data<FP>& pressure_setting__):
-		pressure_setting_{pressure_setting__},
-		rho_setting_{pressure_setting__ * df_info<FP,Descriptor>::inv_cs2}
+	component(const saw::data<FP>& density_setting__):
+		rho_setting_{density_setting__}
 	{}
 
 	template<typename CellFieldSchema>
-	void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step){
+	void apply(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
 		using dfi = df_info<FP,Descriptor>;
+		constexpr int known_dir = Dir ? 1 : -1;
 
 		bool is_even = ((time_step.get() % 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">();
+		auto& info_f = field.template get<"info">();
+		
+		// auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">();
+		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
 
+		auto& dfs_old = dfs_old_f.at(index);
 		/**
 		 * 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});
+			if(c_k[0u]*known_dir <= 0){
+				sum_df += dfs_old.at({k});
 			}
 		}
 
 		/**
 		 * 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];
+			if(c_k[0u]*known_dir > 0){
+				sum_unknown_dfs += dfs_old.at({k}) * c_k[0u];
 			}
 		}
 
@@ -165,13 +127,13 @@ public:
 		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}));
+			dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data<FP>{2.0 / 3.0} * rho_setting_ * vel_x;
+			dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u}));
+			dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({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}));
+			dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data<FP>{2.0 / 3.0} * rho_setting_ * vel_x;
+			dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u}));
+			dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u}));
 		}
 	}
 };
diff --git a/lib/core/c++/chunk.hpp b/lib/core/c++/chunk.hpp
index 5d20faa..a1f2451 100644
--- a/lib/core/c++/chunk.hpp
+++ b/lib/core/c++/chunk.hpp
@@ -35,6 +35,7 @@ using SuperChunk = Array<ChunkSchema,Dim>;
 }
 
 namespace saw {
+
 template<typename Sch, uint64_t Ghost, uint64_t... Sides, typename Encode>
 class data<kel::lbm::sch::Chunk<Sch,Ghost,Sides...>,Encode> final {
 public:
@@ -57,19 +58,21 @@ public:
 		return data<InnerSchema,Encode>::get_dims();
 	}
 
-	data<ValueSchema, Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
+	static constexpr auto to_ghost_index(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
 		std::decay_t<decltype(index)> ind;
 		for(uint64_t i = 0u; i < sizeof...(Sides); ++i){
 			ind.at({i}) = index.at({i}) + Ghost;
 		}
+		return ind;
+	}
+
+	data<ValueSchema, Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
+		std::decay_t<decltype(index)> ind = to_ghost_index(index);
 		return values_.at(ind);
 	}
 
 	const data<ValueSchema, Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
-		std::decay_t<decltype(index)> ind;
-		for(uint64_t i = 0u; i < sizeof...(Sides); ++i){
-			ind.at({i}) = index.at({i}) + Ghost;
-		}
+		std::decay_t<decltype(index)> ind = to_ghost_index(index);
 		return values_.at(ind);
 	}
 
diff --git a/lib/core/c++/collision.hpp b/lib/core/c++/collision.hpp
index 85c3af9..ed26a08 100644
--- a/lib/core/c++/collision.hpp
+++ b/lib/core/c++/collision.hpp
@@ -3,6 +3,7 @@
 #include "macroscopic.hpp"
 #include "component.hpp"
 #include "equilibrium.hpp"
+#include "chunk.hpp"
 
 namespace kel {
 namespace lbm {
diff --git a/lib/core/c++/lbm.hpp b/lib/core/c++/lbm.hpp
index 00f153a..5334e4e 100644
--- a/lib/core/c++/lbm.hpp
+++ b/lib/core/c++/lbm.hpp
@@ -13,6 +13,7 @@
 #include "equilibrium.hpp"
 #include "iterator.hpp"
 #include "macroscopic.hpp"
+#include "stream.hpp"
 #include "write_vtk.hpp"
 #include "util.hpp"
 
diff --git a/lib/core/c++/stream.hpp b/lib/core/c++/stream.hpp
new file mode 100644
index 0000000..d217373
--- /dev/null
+++ b/lib/core/c++/stream.hpp
@@ -0,0 +1,25 @@
+#pragma once
+
+#include "component.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct Stream {};
+}
+
+template<typename T, typename Descriptor, typename Encode>
+class component<T,Descriptor, cmpt::Stream, Encode> final {
+private:
+public:
+	static constexpr saw::string_literal name = "streaming";
+	static constexpr saw::string_literal after = "collide";
+	static constexpr saw::string_literal before = "";
+
+	template<typename CellFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step) const {
+
+	}
+};
+}
+}