From e67393d7fcc6e9aa1dfe10e3a6a0d087b9621678 Mon Sep 17 00:00:00 2001
From: "Claudius \"keldu\" Holeksa" <mail@keldu.de>
Date: Wed, 8 Jul 2026 17:32:12 +0200
Subject: Force calculated

---
 examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp |  1 -
 examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp  |  1 -
 examples/stokes_drag_particle_2d_psm_gpu/sim.cpp   |  3 ---
 modules/core/c++/psm.hpp                           | 16 ++++++++++++++++
 4 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp
index 9e96219..0d6ebec 100644
--- a/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp
+++ b/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp
@@ -135,7 +135,6 @@ saw::error_or<void> setup_initial_conditions(
 			por_f.at(index).at({}) = {0};
 			rho.at({}) = {1};
 			auto& vel = vel_f.at(index);
-			vel.at({{0u}}) = 0.0002f;
 			auto eq = equilibrium<T,Desc>(rho,vel);
 
 			df = eq;
diff --git a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp
index 9c28d88..6166aa1 100644
--- a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp
+++ b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp
@@ -132,7 +132,6 @@ saw::error_or<void> setup_initial_conditions(
 			por_f.at(index).at({}) = {1};
 			rho.at({}) = {1};
 			auto& vel = vel_f.at(index);
-			vel.at({{0u}}) = 0.0002f;
 			auto eq = equilibrium<T,Desc>(rho,vel);
 
 			df = eq;
diff --git a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp
index 135bde4..3220a63 100644
--- a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp
+++ b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp
@@ -134,7 +134,6 @@ saw::error_or<void> setup_initial_conditions(
 			por_f.at(index).at({}) = {1};
 			rho.at({}) = {1};
 			auto& vel = vel_f.at(index);
-			vel.at({{0u}}) = 0.0002f;
 
 			auto eq = equilibrium<T,Desc>(rho,vel);
 
@@ -178,7 +177,6 @@ saw::error_or<void> step(
 		component<T,Desc,cmpt::PSM,encode::Sycl<saw::encode::Native>> collision{0.8};
 		component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb;
 		component<T,Desc,cmpt::AntiBounceBack<0u>,encode::Sycl<saw::encode::Native>> abb;
-		component<T,Desc,cmpt::HlbmOneParticleMomentumExchange,encode::Sycl<saw::encode::Native>> one_part_moment;
 
 		saw::data<sch::Scalar<T>> eps;
 		eps.at({}).set(1.5f);
@@ -223,7 +221,6 @@ saw::error_or<void> step(
 				break;
 				case 2u:
 					opa.apply(macros,index,t_i);
-					one_part_moment.apply(fields,macros,index,t_i);
 					collision.apply(fields,macros,index,t_i);
 					fg.apply(fields,macros,index,t_i);
 					break;
diff --git a/modules/core/c++/psm.hpp b/modules/core/c++/psm.hpp
index 10e8f6e..30a00da 100644
--- a/modules/core/c++/psm.hpp
+++ b/modules/core/c++/psm.hpp
@@ -60,6 +60,22 @@ public:
 			uint64_t i_opp = dfi::opposite_index[i];
 			dfs.at({i}) = dfs_cpy.at({i}) + frequency_ * (eq.at(i) - dfs_cpy.at({i})) * porous.at({}) + (dfs_cpy.at({i_opp}) - dfs_cpy.at({i}) ) * flip_porous;
 		}
+
+		auto& force_f = macros.template get<"force">();
+		auto& force = force_f.at(index);
+		for(uint64_t k{0u}; k < Descriptor::D; ++k){
+			force.at({{k}}).set(0);
+		}
+
+		for(uint64_t i{0u}; i < Descriptor::Q; ++i){
+			uint64_t i_opp = dfi::opposite_index[i];
+			auto dfs_diff = dfs.at({i}) - dfs.at({i_opp});
+			for(uint64_t k{0u}; k < Descriptor::D; ++k){
+				force.at({{k}}) = force.at({{k}}) + dfs_diff * dfi::directions[i][k];
+			}
+		}
+
+		force = force * porous;
 	}
 };
 
-- 
cgit v1.2.3