Dangling changes

Fixed ForcedGuo collision

5 files changed, 31 insertions, 19 deletions

diff --git a/c++/boundary.hpp b/c++/boundary.hpp
index 3cf6f3f..bc4859c 100644
--- a/c++/boundary.hpp
+++ b/c++/boundary.hpp
@@ -51,7 +51,7 @@ public:
 		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 = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
 
 		using dfi = df_info<T,Descriptor>;
 
diff --git a/c++/collision.hpp b/c++/collision.hpp
index df9a734..cd6cf14 100644
--- a/c++/collision.hpp
+++ b/c++/collision.hpp
@@ -17,11 +17,11 @@ template<typename T, typename Descriptor>
 class component<T, Descriptor, cmpt::BGK> {
 private:
 	typename saw::native_data_type<T>::type relaxation_;
-	typename saw::native_data_type<T>::type frequency_;
+	saw::data<T> frequency_;
 public:
 	component(typename saw::native_data_type<T>::type relaxation__):
 		relaxation_{relaxation__},
-		frequency_{1 / relaxation_}
+		frequency_{typename saw::native_data_type<T>::type(1) / relaxation_}
 	{}
 
 	using Component = cmpt::BGK;
@@ -57,7 +57,8 @@ public:
 		auto eq = equilibrium<T,Descriptor>(rho,vel);
 
 		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
-			dfs_old({i}).set(dfs_old({i}).get() + (1.0 / relaxation_) * (eq.at(i).get() - dfs_old({i}).get()));
+			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()));
 		}
 	}
 };
@@ -66,11 +67,11 @@ template<typename T, typename Descriptor>
 class component<T, Descriptor, cmpt::BGKGuo> {
 private:
 	typename saw::native_data_type<T>::type relaxation_;
-	typename saw::native_data_type<T>::type frequency_;
+	saw::data<T> frequency_;
 public:
 	component(typename saw::native_data_type<T>::type relaxation__):
 		relaxation_{relaxation__},
-		frequency_{1 / relaxation_}
+		frequency_{typename saw::native_data_type<T>::type(1) / relaxation_}
 	{}
 
 	using Component = cmpt::BGKGuo;
@@ -89,7 +90,7 @@ public:
 
 	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);
+		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">();
@@ -105,22 +106,21 @@ public:
 		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};
-			saw::data<T> ci_dot_F{0};
-			saw::data<T> u_dot_F{0};
 			for(uint64_t d = 0u; d < Descriptor::D; ++d){
-				ci_dot_u = ci_dot_u + vel.at({d}) * dfi::directions[i][d];
-				ci_dot_F = ci_dot_F + force({d}) * dfi::directions[i][d];
-				u_dot_F = u_dot_F + vel.at({d}) * force({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));
 
-			// auto df_old = dfs_old({i});
-			saw::data<T> F_i = saw::data<T>{dfi::weights[i]} * (ci_dot_F / dfi::cs2 + ci_dot_u * ci_dot_F / (dfi::cs2 * dfi::cs2));
-			if(F_i.get() != 0){
-				std::cerr<<"AAAH";
+			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}) + saw::data<T>{frequency_} * (eq.at(i) - dfs_old({i}) ) - F_i * saw::data<T>{1.0f - 0.5f * frequency_};
+
+			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++/particle/particle.hpp b/c++/particle/particle.hpp
index a55b06e..ae7769a 100644
--- a/c++/particle/particle.hpp
+++ b/c++/particle/particle.hpp
@@ -27,7 +27,7 @@ using ParticleMask = Struct<
 template<typename T, uint64_t D>
 using Particle = Struct<
 	Member<ParticleRigidBody<T,D>, "rigid_body">,
-	Member<ParticleMask<Float32,D>, "mask">,
+	Member<ParticleMask<T,D>, "mask">,
 	Member<T, "size">
 >;
 }
@@ -84,6 +84,7 @@ public:
 
 	template<typename LbmLattice>
 	void update_particle_border(saw::data<LbmLattice>& latt){
+		(void) latt;
 		for(auto& iter : particles_){
 			auto& par = iter;
 
diff --git a/c++/statistics.hpp b/c++/statistics.hpp
new file mode 100644
index 0000000..c07ccb7
--- /dev/null
+++ b/c++/statistics.hpp
@@ -0,0 +1,11 @@
+#pragma once
+
+namespace kel {
+namespace lbm {
+template<typename T>
+class statistics {
+private:
+public:
+};
+}
+}
diff --git a/c++/write_vtk.hpp b/c++/write_vtk.hpp
index eb2545a..c9ecc99 100644
--- a/c++/write_vtk.hpp
+++ b/c++/write_vtk.hpp
@@ -81,7 +81,7 @@ struct lbm_vtk_writer<sch::Array<sch::Struct<sch::Member<StructT,StructN>...> ,
 	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();
+		// 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);