Weird missing vtk writes

7 files changed, 87 insertions, 104 deletions

diff --git a/lib/core/c++/collision.hpp b/lib/core/c++/collision.hpp
index ed26a08..a05e263 100644
--- a/lib/core/c++/collision.hpp
+++ b/lib/core/c++/collision.hpp
@@ -61,6 +61,28 @@ public:
 			dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i}));
 		}
 	}
+	
+	template<typename CellFieldSchema, typename MacroFieldSchema>
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+		bool is_even = ((time_step.get() % 2) == 0);
+
+		// 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& rho_f = macros.template get<"density">();
+		auto& vel_f = macros.template get<"velocity">();
+	
+		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
+		saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);
+
+		compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel);
+		
+		auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i}));
+		}
+	}
 };
 
 template<typename T, typename Descriptor, typename Encode>
diff --git a/lib/core/c++/equilibrium.hpp b/lib/core/c++/equilibrium.hpp
index bb55d00..7d1324e 100644
--- a/lib/core/c++/equilibrium.hpp
+++ b/lib/core/c++/equilibrium.hpp
@@ -5,7 +5,7 @@
 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){
+saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<sch::Scalar<T>> rho, saw::data<sch::Vector<T,Descriptor::D>> vel){
 	using dfi = df_info<T, Descriptor>;
 
 	saw::data<sch::FixedArray<T,Descriptor::Q>> eq;
@@ -17,7 +17,7 @@ saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<T> rho, saw::
 	//  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);
+		vel_vel = vel_vel + vel.at({{j}}) * vel.at({{j}});
 	}
 
 	/**
@@ -27,13 +27,13 @@ saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<T> rho, saw::
 		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])});
+			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() * 
+			dfi::weights[i] * rho.at({}).get() * 
 			(
 			 1.0
 			 + vel_c_cs2.get()
diff --git a/lib/core/c++/macroscopic.hpp b/lib/core/c++/macroscopic.hpp
index 1532873..19ab3e9 100644
--- a/lib/core/c++/macroscopic.hpp
+++ b/lib/core/c++/macroscopic.hpp
@@ -4,113 +4,32 @@
 
 namespace kel {
 namespace lbm {
-template<typename T, typename Desc>
-void compute_rho_u (
-		const saw::data<sch::FixedArray<T,Desc::Q>>& dfs,
-		saw::data<T>& rho,
-		saw::data<sch::FixedArray<T,Desc::D>>& vel
-	)
-{
-	using dfi = df_info<T, Desc>;
-
-	rho = 0;
-	for(size_t j = 0; j < Desc::D; ++j){
-		vel.at({{j}}).set(0);
-	}
-
-	for(size_t j = 0; j < Desc::Q; ++j){
-		rho = rho + dfs.at({j});
-		for(size_t i = 0; i < Desc::D; ++i){
-			vel.at({{i}}) = vel.at({{i}}) + dfi::directions[j][i] * dfs.at({j}).get();
-		}
-	}
-
-	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,
-		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,
+		const saw::data<sch::FixedArray<T,Desc::Q>>& dfs,
+		saw::ref<saw::data<sch::Scalar<T>>> rho,
 		saw::ref<saw::data<sch::Vector<T,Desc::D>>> vel
 	)
 {
 	using dfi = df_info<T, Desc>;
 
-	rho().set(0);
+	rho().at({}).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);
+		rho().at({}) = rho().at({}) + dfs.at({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);
+			vel().at({{i}}) = vel().at({{i}}) + saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[j][i])} * dfs.at({j});
 		}
 	}
 
 	for(size_t i = 0; i < Desc::D; ++i){
-		vel().at({{i}}) = vel().at({{i}}) / rho();
+		vel().at({{i}}) = vel().at({{i}}) / rho().at({});
 	}
 }
 }
diff --git a/lib/core/c++/stream.hpp b/lib/core/c++/stream.hpp
index 8c0342b..dc7cfb3 100644
--- a/lib/core/c++/stream.hpp
+++ b/lib/core/c++/stream.hpp
@@ -17,19 +17,44 @@ public:
 	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 {
+	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<T,Descriptor>;
+		auto g_index = index;
+
+		for(uint64_t i = 0u; i < Descriptor::D; ++i){
+			++g_index.at({{i}});
+		}
+
 		bool is_even = ((time_step.get() % 2) == 0);
 
 		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+		auto& dfs_new_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">();
 		auto info_f = field.template get<"info">();
 
 		auto info_meta = info_f.get_dims();
 
+		/*
 		bool border = false;
 		for(uint64_t i = 0u; i < Descriptor::D; ++i){
 			auto ind_i = index.at({i});
 			border |= (ind_i.get()) == 0u or (ind_i == info_meta.at({i}));
 		}
+
+		if (not border){
+		}
+		*/
+
+		auto& dfs_new = dfs_new_f.ghost_at(g_index);
+
+		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			auto dir = dfi::directions[dfi::opposite_index[i]];
+			auto g_index_nb = g_index;
+			for(uint64_t z = 0u; z < Descriptor::D; ++z){
+				g_index_nb.at({z}).set(g_index.at({z}).get() + dir[z]);
+			}
+
+			dfs_new.at({i}) = dfs_old_f.ghost_at(g_index_nb).at({i});
+		}
 	}
 };
 }
diff --git a/lib/core/c++/write_vtk.hpp b/lib/core/c++/write_vtk.hpp
index e84cf9d..3c60a66 100644
--- a/lib/core/c++/write_vtk.hpp
+++ b/lib/core/c++/write_vtk.hpp
@@ -100,6 +100,8 @@ struct lbm_vtk_writer<sch::Chunk<T,Ghost,D...>> {
 			}
 		}
 
+		vtk_file<<"\n";
+
 		return saw::make_void();
 	}
 };
@@ -129,6 +131,21 @@ struct lbm_vtk_writer<sch::Vector<T,D>> {
 	}
 };
 
+template<typename T>
+struct lbm_vtk_writer<sch::Scalar<T>> {
+	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::Scalar<T>>& field){
+		// vtk_file<<"VECTORS "<<name<<" float\n";
+		vtk_file<<field.at({}).get();
+		vtk_file<<"\n";
+		return saw::make_void();
+	}
+};
+
 template<typename... MemberT, saw::string_literal... Keys, uint64_t... Ghost, uint64_t... Dims>
 struct lbm_vtk_writer<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>> final  {
 	template<uint64_t i>
@@ -176,7 +193,6 @@ struct lbm_vtk_writer<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,
 			static_assert(saw::ct_multiply<uint64_t,Dims...>::value > 0u, "Invalid Dim size resulting in length 0u");
 
 			for(saw::data<sch::UInt64> i{0u}; i.get() < sizeof...(Dims); ++i){
-				pd_size = pd_size * meta.at(i);
 				vtk_file << " " << meta.at(i).get(); 
 			}
 			for(saw::data<sch::UInt64> i{sizeof...(Dims)}; i.get() < 3u; ++i){
diff --git a/lib/core/tests/equilibrium.cpp b/lib/core/tests/equilibrium.cpp
index 9201e55..20a1f08 100644
--- a/lib/core/tests/equilibrium.cpp
+++ b/lib/core/tests/equilibrium.cpp
@@ -11,10 +11,11 @@ void check_equilibrium(){
 
 	using dfi = lbm::df_info<lbm::sch::Float64,Descriptor>;
 
-	saw::data<lbm::sch::Float64> rho{1.0};
-	saw::data<lbm::sch::FixedArray<lbm::sch::Float64,Descriptor::D>> vel;
+	saw::data<lbm::sch::Scalar<lbm::sch::Float64>> rho;
+	rho.at({}).set(1.0);
+	saw::data<lbm::sch::Vector<lbm::sch::Float64,Descriptor::D>> vel;
 	for(saw::data<lbm::sch::UInt64> i{0u}; i.get() < Descriptor::D; ++i){
-		vel.at(i) = {0.0};
+		vel.at({{i}}) = {0.0};
 	}
 	auto eq = lbm::equilibrium<lbm::sch::Float64,Descriptor>(rho,vel);
 
diff --git a/lib/sycl/c++/data.hpp b/lib/sycl/c++/data.hpp
index 6a4f068..10c500d 100644
--- a/lib/sycl/c++/data.hpp
+++ b/lib/sycl/c++/data.hpp
@@ -52,11 +52,11 @@ public:
 		return saw::data<schema::FixedArray<schema::UInt64, sizeof...(Dims)>>{{Dims...}};
 	}
 
-	constexpr data<Sch,Encode>& at(data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index){
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index){
 		return values_[kel::lbm::flatten_index<schema::UInt64,sizeof...(Dims)>::apply(index,get_dims()).get()];
 	}
 	
-	constexpr data<Sch,Encode>& at(data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index) const{
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index) const{
 		return values_[kel::lbm::flatten_index<schema::UInt64,sizeof...(Dims)>::apply(index,get_dims()).get()];
 	}
 
@@ -88,11 +88,11 @@ public:
 		return saw::data<schema::FixedArray<schema::UInt64, sizeof...(Dims)>>{{Dims...}};
 	}
 
-	constexpr data<Sch,Encode>& at(data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index){
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index){
 		return values_[kel::lbm::flatten_index<schema::UInt64,sizeof...(Dims)>::apply(index,get_dims()).get()];
 	}
 	
-	constexpr data<Sch,Encode>& at(data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index) const{
+	constexpr data<Sch,Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Dims)>>& index) const{
 		return values_[kel::lbm::flatten_index<schema::UInt64,sizeof...(Dims)>::apply(index,get_dims()).get()];
 	}
 
@@ -115,11 +115,11 @@ public:
 		values_{q__}
 	{}
 
-	data<ValueSchema, kel::lbm::encode::Sycl<Encode>>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
+	constexpr data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
 		return values_.at(index);
 	}
 
-	data<ValueSchema, kel::lbm::encode::Sycl<Encode>>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
+	constexpr data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
 		return values_.at(index);
 	}
 
@@ -178,11 +178,11 @@ public:
 		values_{values__}
 	{}
 
-	data<ValueSchema, kel::lbm::encode::Sycl<Encode>>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
+	constexpr data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){
 		return values_.at(index);
 	}
 
-	data<ValueSchema, kel::lbm::encode::Sycl<Encode>>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
+	constexpr data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const {
 		return values_.at(index);
 	}