Dangling changes

Fixed ForcedGuo collision

1 files changed, 16 insertions, 16 deletions

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_);
 		}
 	}
 };