Merge branch 'dev'

1 files changed, 22 insertions, 32 deletions

diff --git a/lib/core/c++/boundary.hpp b/lib/core/c++/boundary.hpp
index adb473d..4dbbdf8 100644
--- a/lib/core/c++/boundary.hpp
+++ b/lib/core/c++/boundary.hpp
@@ -174,14 +174,13 @@ class component<FP, Descriptor, cmpt::ZouHeHorizontal<Dir>, Encode> final {
 private:
 	saw::data<FP> rho_setting_;
 public:
-	component(const saw::data<FP>& density_setting__):
-		rho_setting_{density_setting__}
+	component(const saw::data<FP>& rho_setting__):
+		rho_setting_{rho_setting__}
 	{}
 
 	template<typename CellFieldSchema>
 	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);
 
@@ -191,42 +190,33 @@ public:
 		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()];
-
-			if(c_k[0u]*known_dir <= 0){
-				sum_df += dfs_old.at({k});
-			}
-		}
-
-		/**
-		 * Get the sum of the unknown dfs and precalculate the direction
-		 */
-		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()];
-			if(c_k[0u]*known_dir < 0){
-				sum_unknown_dfs += dfs_old.at({k}) * c_k[0u];
+		auto rho_vel_x = [&]() -> saw::data<FP> {
+			if constexpr (Dir){
+				auto S = dfs_old.at({0u})
+					+ dfs_old.at({3u}) + dfs_old.at({4u})
+					+ (dfs_old.at({1u}) + dfs_old.at({5u}) + dfs_old.at({7u})) * 2u;
+				return rho_setting_ - S;
 			}
-		}
-
-		auto vel_x = sum_unknown_dfs / rho_setting_;
+			else if constexpr (not Dir) {
+				auto S = dfs_old.at({0u})
+					+ dfs_old.at({3u}) + dfs_old.at({4u})
+					+ (dfs_old.at({2u}) + dfs_old.at({6u}) + dfs_old.at({8u})) * 2u;
+				return S - rho_setting_;
+			}
+			return {};
+		}();
 
 		static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly");
 
 		if constexpr (Dir) {
-			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}));
+			dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data<FP>{2.0 / 3.0} * rho_vel_x;
+			dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data<FP>{1.0 / 6.0} * rho_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_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u}));
 		}else if constexpr (not Dir){
-			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}));
+			dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data<FP>{2.0 / 3.0} * rho_vel_x;
+			dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data<FP>{1.0 / 6.0} * rho_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_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u}));
 		}
 	}
 };