Fixing FpLBM, because I messed up guo forcing

2 files changed, 26 insertions, 1 deletions

diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp
index c9ed6e9..d8c075b 100644
--- a/modules/core/c++/fplbm.hpp
+++ b/modules/core/c++/fplbm.hpp
@@ -92,7 +92,7 @@ public:
 		auto& force_f = macros.template get<"force">();
 		auto& force = force_f.at(index);
 
-		vel = vel + force * half;
+		vel = vel + force * half / rho;
 
 		// compute_rho_u<T,Descriptor>(dfs,rho,vel);
 		auto eq = equilibrium<T,Descriptor>(rho,vel);
diff --git a/modules/core/c++/macroscopic.hpp b/modules/core/c++/macroscopic.hpp
index 768fef0..9db1538 100644
--- a/modules/core/c++/macroscopic.hpp
+++ b/modules/core/c++/macroscopic.hpp
@@ -33,5 +33,30 @@ void compute_rho_u (
 		vel().at({{i}}) = vel().at({{i}}) / rho().at({});
 	}
 }
+
+/**
+ * Calculate the macroscopic variables rho and u in Lattice Units.
+ */
+template<typename T, typename Desc>
+void compute_rho_momentum (
+		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>>> mom
+	)
+{
+	using dfi = df_info<T, Desc>;
+
+	rho().at({}).set(0);
+	for(uint64_t i = 0; i < Desc::D; ++i){
+		mom().at({{i}}).set(0);
+	}
+
+	for(size_t j = 0; j < Desc::Q; ++j){
+		rho().at({}) = rho().at({}) + dfs.at({j});
+		for(size_t i = 0; i < Desc::D; ++i){
+			mom().at({{i}}) = mom().at({{i}}) + saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[j][i])} * dfs.at({j});
+		}
+	}
+}
 }
 }