Renaming velocity to momentum. for now

44 files changed, 1179 insertions, 200 deletions

diff --git a/examples/cavity_2d/cavity_2d.cpp b/examples/cavity_2d/cavity_2d.cpp
index 702661a..af22796 100644
--- a/examples/cavity_2d/cavity_2d.cpp
+++ b/examples/cavity_2d/cavity_2d.cpp
@@ -42,7 +42,7 @@ using CellStruct = Struct<
 
 template<typename T, uint64_t D>
 using MacroStruct = Struct<
-	Member<FixedArray<T,D>, "velocity">,
+	Member<FixedArray<T,D>, "momentum">,
 	Member<T, "pressure">
 >;
 
@@ -406,7 +406,7 @@ int main(){
 				auto dfs = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
 
 				auto& rho = macros.at({{i,j}}).template get<"pressure">();
-				auto& vel = macros.at({{i,j}}).template get<"velocity">();
+				auto& vel = macros.at({{i,j}}).template get<"momentum">();
 				compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel);
 			}, lattice);
 
diff --git a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
index 0fc490a..fe44baf 100644
--- a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
+++ b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
@@ -47,7 +47,7 @@ using CellStruct = CellFieldStruct<
 
 template<typename T, uint64_t D>
 using MacroStruct = Struct<
-	Member<FixedArray<T,D>, "velocity">,
+	Member<FixedArray<T,D>, "momentum">,
 	Member<T, "pressure">
 >;
 
@@ -362,7 +362,7 @@ int main(){
 				auto dfs_field = even_step ? dfs_old : dfs;
 
 				auto& rho = macros.at(index).template get<"pressure">();
-				auto& vel = macros.at(index).template get<"velocity">();
+				auto& vel = macros.at(index).template get<"momentum">();
 				compute_rho_u<sch::T,sch::D2Q9>(dfs_field[j],rho,vel);
 			}, {{0u,0u}}, meta);
 			std::string vtk_f_name{"tmp/cavity_2d_gpu_"};
diff --git a/examples/cavity_3d.cpp b/examples/cavity_3d.cpp
index 64356be..5ae4834 100644
--- a/examples/cavity_3d.cpp
+++ b/examples/cavity_3d.cpp
@@ -46,7 +46,7 @@ using CellStruct = Struct<
 
 template<typename T, uint64_t D>
 using MacroStruct = Struct<
-	Member<FixedArray<T,D>, "velocity">,
+	Member<FixedArray<T,D>, "momentum">,
 	Member<T, "pressure">
 >;
 
@@ -411,7 +411,7 @@ int main(){
 				auto dfs = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
 
 				auto& rho = macros.at({{i,j}}).template get<"pressure">();
-				auto& vel = macros.at({{i,j}}).template get<"velocity">();
+				auto& vel = macros.at({{i,j}}).template get<"momentum">();
 				compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel);
 			}, lattice);
 
diff --git a/examples/meta/meta.cpp b/examples/meta/meta.cpp
index a5e2ae5..7232c61 100644
--- a/examples/meta/meta.cpp
+++ b/examples/meta/meta.cpp
@@ -75,7 +75,7 @@ saw::error_or<void> meta_main(int argc, char** argv){
 				meta_lbm_info<sch::Prec,sch::Descriptor<1u,3u>>(args_str);
 				break;
 			default:
-				return saw::make_error<saw::err::critical>("Second positional needs to be a valid velocity set");
+				return saw::make_error<saw::err::critical>("Second positional needs to be a valid momentum set");
 			}
 		}
 		break;
@@ -89,7 +89,7 @@ saw::error_or<void> meta_main(int argc, char** argv){
 				meta_lbm_info<sch::Prec,sch::Descriptor<2u,9u>>(args_str);
 				break;
 			default:
-				return saw::make_error<saw::err::critical>("Second positional needs to be a valid velocity set");
+				return saw::make_error<saw::err::critical>("Second positional needs to be a valid momentum set");
 			}
 		}
 		break;
@@ -108,7 +108,7 @@ saw::error_or<void> meta_main(int argc, char** argv){
 				meta_lbm_info<sch::Prec,sch::Descriptor<3u,27u>>(args_str);
 				break;
 			default:
-				return saw::make_error<saw::err::critical>("Second positional needs to be a valid velocity set");
+				return saw::make_error<saw::err::critical>("Second positional needs to be a valid momentum set");
 			}
 		}
 		break;
diff --git a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp
index cc7649a..4c9c8bd 100644
--- a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp
+++ b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp
@@ -55,7 +55,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<ScalarChunk<T,Desc>, "porosity">,
 	Member<VelChunk<T,Desc>, "force">
@@ -119,7 +119,7 @@ saw::error_or<void> init(
 	auto& dfs_old_f = fields.template get<"dfs_old">();
 	auto& dfs_f = fields.template get<"dfs">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 	auto& force_f = macros.template get<"force">();
 	
diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
index 9e812e2..bd754bf 100644
--- a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
+++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
@@ -46,7 +46,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<ScalarChunk<T,Desc>, "porosity">
 >;
@@ -109,7 +109,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 	
 	iterator<Desc::D>::apply(
diff --git a/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp
index 2c4c9d1..2daa452 100644
--- a/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp
+++ b/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp
@@ -46,7 +46,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<ScalarChunk<T,Desc>, "porosity">,
 	Member<VelChunk<T,Desc>, "force">
@@ -110,7 +110,7 @@ saw::error_or<void> init(
 	auto& dfs_old_f = fields.template get<"dfs_old">();
 	auto& dfs_f = fields.template get<"dfs">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 	auto& force_f = macros.template get<"force">();
 	
diff --git a/examples/poiseulle_2d/poiseulle_2d.cpp b/examples/poiseulle_2d/poiseulle_2d.cpp
index 3bc7a6c..f9db3df 100644
--- a/examples/poiseulle_2d/poiseulle_2d.cpp
+++ b/examples/poiseulle_2d/poiseulle_2d.cpp
@@ -41,7 +41,7 @@ using CellStruct = Struct<
 
 template<typename T, uint64_t D>
 using MacroStruct = Struct<
-	Member<FixedArray<T,D>, "velocity">,
+	Member<FixedArray<T,D>, "momentum">,
 	Member<T, "pressure">
 >;
 
@@ -375,7 +375,7 @@ int main(int argc, char** argv){
 				auto& dfs = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">();
 
 				auto& rho = macros.at(index).template get<"pressure">();
-				auto& vel = macros.at(index).template get<"velocity">();
+				auto& vel = macros.at(index).template get<"momentum">();
 				compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel);
 				rho = rho * saw::data<sch::T>{dfi::cs2};
 
diff --git a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp
index 1aaa5c0..9963c59 100644
--- a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp
+++ b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp
@@ -36,7 +36,7 @@ using CellStruct = Struct<
 >;
 
 using MacroStruct = Struct<
-	Member<FixedArray<Float64,D2Q9::D>, "velocity">,
+	Member<FixedArray<Float64,D2Q9::D>, "momentum">,
 	Member<Float64, "pressure">
 >;
 
@@ -245,7 +245,7 @@ void step(
 				auto& macro_c = macro_cells[acc_id];
 
 				saw::data<sch::T>& rho = macro_c.template get<"pressure">();
-				saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"velocity">();
+				saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"momentum">();
 
 				compute_rho_u<sch::T,Desc>(dfs_old,rho,vel);
 				auto eq = equilibrium<sch::T,Desc>(rho,vel);
@@ -262,7 +262,7 @@ void step(
 				auto& macro_c = macro_cells[acc_id];
 
 				saw::data<sch::T>& rho = macro_c.template get<"pressure">();
-				saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"velocity">();
+				saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"momentum">();
 
 				compute_rho_u<sch::T,Desc>(dfs_old,rho,vel);
 				auto eq = equilibrium<sch::T,Desc>(rho,vel);
@@ -279,7 +279,7 @@ void step(
 				auto& macro_c = macro_cells[acc_id];
 
 				saw::data<sch::T>& rho = macro_c.template get<"pressure">();
-				saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"velocity">();
+				saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"momentum">();
 
 				compute_rho_u<sch::T,Desc>(dfs_old,rho,vel);
 				auto eq = equilibrium<sch::T,Desc>(rho,vel);
diff --git a/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp b/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp
index 1693c5e..1af0023 100644
--- a/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp
+++ b/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp
@@ -47,7 +47,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<ScalarChunk<T,Desc>, "porosity">
 >;
@@ -123,7 +123,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 	
 	iterator<Desc::D>::apply(
diff --git a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp
index 84ce22b..9ed1646 100644
--- a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp
@@ -45,7 +45,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">
 >;
 
@@ -102,7 +102,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	
 	iterator<Desc::D>::apply(
 		[&](auto& index){
diff --git a/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp
index c6aac92..f45e10e 100644
--- a/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp
@@ -49,7 +49,7 @@ using ForceChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<VelChunk<T,Desc>, "force">,
 	Member<ScalarChunk<T,Desc>, "porosity">
@@ -109,7 +109,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 
 	iterator<Desc::D>::apply(
diff --git a/examples/poiseulle_particles_2d_gpu/common.hpp b/examples/poiseulle_particles_2d_gpu/common.hpp
index 90b7687..b321344 100644
--- a/examples/poiseulle_particles_2d_gpu/common.hpp
+++ b/examples/poiseulle_particles_2d_gpu/common.hpp
@@ -49,7 +49,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<ScalarChunk<T,Desc>, "porosity">
 >;
diff --git a/examples/poiseulle_particles_2d_gpu/init.hpp b/examples/poiseulle_particles_2d_gpu/init.hpp
index cdce25e..b7fdeb0 100644
--- a/examples/poiseulle_particles_2d_gpu/init.hpp
+++ b/examples/poiseulle_particles_2d_gpu/init.hpp
@@ -57,7 +57,7 @@ saw::error_or<void> init(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 	
 	iterator<Desc::D>::apply(
diff --git a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp
index 6cadbd0..0061fa2 100644
--- a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp
@@ -46,7 +46,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<VelChunk<T,Desc>, "force">,
 	Member<ScalarChunk<T,Desc>, "porosity">
@@ -106,7 +106,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 	
 	iterator<Desc::D>::apply(
diff --git a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp
index 640c9ab..d63515d 100644
--- a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp
@@ -46,7 +46,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<VelChunk<T,Desc>, "force">
 >;
@@ -113,7 +113,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& force_f = macros.template get<"force">();
 	
 	iterator<Desc::D>::apply(
@@ -218,7 +218,7 @@ saw::error_or<void> step(
 			
 	q.submit([&](acpp::sycl::handler& h){
 		h.parallel_for(acpp::sycl::range<1u>{1u}, [=](acpp::sycl::id<1u> idx){
-			auto& vel_f = macros.template get<"velocity">();
+			auto& vel_f = macros.template get<"momentum">();
 			auto& dense_f = macros.template get<"density">();
 
 			auto& ps = particles;
diff --git a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
index c2ab675..760f70a 100644
--- a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
@@ -47,7 +47,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<VelChunk<T,Desc>, "force">,
 	Member<ScalarChunk<T,Desc>, "porosity">
@@ -124,7 +124,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 	
 	iterator<Desc::D>::apply(
diff --git a/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp
index 1e6e75f..24832f2 100644
--- a/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp
+++ b/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp
@@ -68,7 +68,7 @@ using CellStruct = Struct<
 
 template<typename T, uint64_t D>
 using MacroStruct = Struct<
-	Member<Vector<T,D>, "velocity">,
+	Member<Vector<T,D>, "momentum">,
 	Member<T, "pressure">,
 	Member<UInt16, "particle">,
 	Member<FixedArray<T,D>, "force">
@@ -369,7 +369,7 @@ void couple_particles_to_lattice(
 
 			// Interpolate this from close U cells.
 			// For now pick the closest U
-			auto& closest_u = macros.at(p_cell_pos).template get<"velocity">();
+			auto& closest_u = macros.at(p_cell_pos).template get<"momentum">();
 			// auto p_shift = closest_u;
 
 			// p_pos - p_pos_old
@@ -380,7 +380,7 @@ void couple_particles_to_lattice(
 				p_vel_rel.at({{1u}}) = p_vel_rel.at({{1u}}) / p_mass;
 			}
 
-			// Add relative velocity to particle acceleration (Since our timestep is 1 we don't need to do anything else)
+			// Add relative momentum to particle acceleration (Since our timestep is 1 we don't need to do anything else)
 
 			p_acc = p_acc + p_vel_rel;
 
@@ -463,7 +463,7 @@ void couple_particles_to_lattice(
 			}
 
 
-			/* So I want to include the relative velocity from both particles, but this is a bit hard considering I just assumed 0 velocity from the other party
+			/* So I want to include the relative momentum from both particles, but this is a bit hard considering I just assumed 0 momentum from the other party
 			else if(n_macro_cell_particle.get() != (i.get() + 1u) and n_macro_cell_particle.get() > 0u){
 				// Generally compare
 			}
@@ -721,7 +721,7 @@ saw::error_or<void> real_main(int argc, char** argv){
 				auto& macro_cell = macros.at(index);
 
 				auto& rho = macro_cell.template get<"pressure">();
-				auto& vel = macro_cell.template get<"velocity">();
+				auto& vel = macro_cell.template get<"momentum">();
 				auto& force = macro_cell.template get<"force">();
 
 				auto& part_mask = macro_cell.template get<"particle">();
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix b/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix
new file mode 100644
index 0000000..d5361c2
--- /dev/null
+++ b/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix
@@ -0,0 +1,44 @@
+{ lib
+, stdenv
+, scons
+, clang-tools
+, forstio
+, python3
+, pname
+, version
+, adaptive-cpp
+, kel
+, particle_coupling ? "fplbm"
+}:
+
+let
+	slip_txt = if slip then "slip" else "noslip";
+in stdenv.mkDerivation {
+	pname = "${pname}-examples-stdkr_2d_gpu_${particle_coupling}";
+	inherit version;
+	src = ./..;
+
+	nativeBuildInputs = [
+		scons
+		clang-tools
+		python3
+	];
+
+	buildInputs = [
+		forstio.core
+		forstio.async
+		forstio.codec
+		forstio.codec-unit
+		forstio.io
+    forstio.remote
+    forstio.remote-filesystem
+		forstio.codec-json
+		adaptive-cpp
+		kel.lbm.core
+		kel.lbm.sycl
+	];
+
+	preferLocalBuild = true;
+
+	outputs = [ "out" "dev" ];
+}
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/SConscript b/examples/schaefer_turek_durst_krause_rannbacher/SConscript
new file mode 100644
index 0000000..4483d58
--- /dev/null
+++ b/examples/schaefer_turek_durst_krause_rannbacher/SConscript
@@ -0,0 +1,34 @@
+#!/bin/false
+
+import os
+import os.path
+import glob
+
+
+Import('env')
+
+dir_path = Dir('.').abspath
+
+# Environment for base library
+examples_env = env.Clone();
+examples_env['CXX'] = 'syclcc-clang';
+examples_env['CXXFLAGS'] += ['-O3'];
+
+examples_env.sources = sorted(glob.glob(dir_path + "/*.cpp"))
+examples_env.headers = sorted(glob.glob(dir_path + "/*.hpp"))
+
+env.sources += examples_env.sources;
+env.headers += examples_env.headers;
+
+# Cavity2D
+examples_objects = [];
+examples_env.add_source_files(examples_objects, ['sim.cpp'], shared=False);
+examples_env.poiseulle_2d_gpu = examples_env.Program('#bin/poiseulle_particles_2d_hlbm_gpu', [examples_objects]);
+
+# Set Alias
+env.examples = [
+	examples_env.poiseulle_2d_gpu
+];
+env.Alias('examples', env.examples);
+env.targets += ['examples'];
+env.Install('$prefix/bin/', env.examples);
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/SConstruct b/examples/schaefer_turek_durst_krause_rannbacher/SConstruct
new file mode 100644
index 0000000..0611b67
--- /dev/null
+++ b/examples/schaefer_turek_durst_krause_rannbacher/SConstruct
@@ -0,0 +1,81 @@
+#!/usr/bin/env python3
+
+import sys
+import os
+import os.path
+import glob
+import re
+
+
+if sys.version_info < (3,):
+    def isbasestring(s):
+        return isinstance(s,basestring)
+else:
+    def isbasestring(s):
+        return isinstance(s, (str,bytes))
+
+def add_kel_source_files(self, sources, filetype, lib_env=None, shared=False, target_post=""):
+
+    if isbasestring(filetype):
+        dir_path = self.Dir('.').abspath
+        filetype = sorted(glob.glob(dir_path+"/"+filetype))
+
+    for path in filetype:
+        target_name = re.sub( r'(.*?)(\.cpp|\.c\+\+)', r'\1' + target_post, path )
+        if shared:
+            target_name+='.os'
+            sources.append( self.SharedObject( target=target_name, source=path ) )
+        else:
+            target_name+='.o'
+            sources.append( self.StaticObject( target=target_name, source=path ) )
+    pass
+
+def isAbsolutePath(key, dirname, env):
+	assert os.path.isabs(dirname), "%r must have absolute path syntax" % (key,)
+
+env_vars = Variables(
+	args=ARGUMENTS
+)
+
+env_vars.Add('prefix',
+	help='Installation target location of build results and headers',
+	default='/usr/local/',
+	validator=isAbsolutePath
+)
+
+env_vars.Add('build_examples',
+	help='If examples should be built',
+	default="true"
+)
+
+env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[],
+    CPPDEFINES=['SAW_UNIX'],
+    CXXFLAGS=[
+        '-std=c++20',
+        '-g',
+        '-Wall',
+        '-Wextra'
+    ],
+    LIBS=[
+        'forstio-core',
+        'forstio-async',
+        'forstio-io'
+    ]
+);
+env.__class__.add_source_files = add_kel_source_files
+env.Tool('compilation_db');
+env.cdb = env.CompilationDatabase('compile_commands.json');
+
+env.objects = [];
+env.sources = [];
+env.headers = [];
+env.targets = [];
+
+Export('env')
+SConscript('SConscript')
+
+env.Alias('cdb', env.cdb);
+env.Alias('all', [env.targets]);
+env.Default('all');
+
+env.Alias('install', '$prefix')
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/common.hpp b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp
new file mode 100644
index 0000000..2d51614
--- /dev/null
+++ b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp
@@ -0,0 +1,67 @@
+#pragma once
+
+#include <kel/lbm/sycl/lbm.hpp>
+#include <kel/lbm/lbm.hpp>
+#include <kel/lbm/particle.hpp>
+#include <kel/lbm/particle/particle.hpp>
+
+#include <forstio/io/io.hpp>
+#include <forstio/remote/filesystem/easy.hpp>
+#include <forstio/codec/json/json.hpp>
+#include <forstio/codec/simple.hpp>
+
+namespace kel {
+namespace lbm {
+
+/**
+ * length 2.2m
+ * viscosity 10^-3
+ */
+
+constexpr uint64_t dim_y = 256ul;
+constexpr uint64_t dim_x = dim_y * 20ul;
+
+constexpr uint64_t particle_amount = 1ul;
+
+namespace sch {
+using namespace saw::schema;
+
+using InfoChunk = Chunk<UInt8, 0u, dim_x, dim_y>;
+
+template<typename T, typename Desc>
+using DfChunk = Chunk<FixedArray<T,Desc::Q>, 1u, dim_x, dim_y>;
+
+template<typename T, typename Desc>
+using ScalarChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
+
+template<typename T, typename Desc>
+using VectorChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>;
+
+template<typename T, typename Desc>
+using ChunkStruct = Struct<
+	Member<InfoChunk, "info">,
+	Member<DfChunk<T,Desc>, "dfs">,
+	Member<DfChunk<T,Desc>, "dfs_old">,
+	Member<VectorChunk<T,Desc>, "particle_N">,
+	Member<ScalarChunk<T,Desc>, "particle_D">
+>;
+
+template<typename T, typename Desc>
+using VelChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>;
+
+template<typename T>
+using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
+
+template<typename T, typename Desc>
+using MacroStruct = Struct<
+	Member<VelChunk<T,Desc>, "momentum">,
+	Member<RhoChunk<T>, "density">,
+	Member<ScalarChunk<T,Desc>, "porosity">
+>;
+
+template<typename T, typename Desc>
+using ParticleSpheroidGroup = ParticleGroup<T,Desc::D,1u,coll::Spheroid<T>>;
+}
+
+}
+}
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/init.hpp b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp
new file mode 100644
index 0000000..b7fdeb0
--- /dev/null
+++ b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp
@@ -0,0 +1,127 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+
+template<typename T, typename Desc, typename Coll>
+saw::error_or<void> init(
+		const converter<T>& conv,
+		saw::data<sch::ChunkStruct<T,Desc>>& fields,
+		saw::data<sch::MacroStruct<T,Desc>>& macros,
+		saw::data<sch::ParticleSpheroidGroup<T,Desc>>& particles
+){
+	(void) conv;
+
+	auto& info_f = fields.template get<"info">();
+	auto& porous_f = macros.template get<"porosity">();
+	// Set everything as walls
+	iterator<Desc::D>::apply(
+		[&](auto& index){
+			info_f.at(index).set(1u);
+		},
+		{},
+		info_f.get_dims(),
+		{}
+	);
+	// Fluid
+	iterator<Desc::D>::apply(
+		[&](auto& index){
+			info_f.at(index).set(2u);
+		},
+		{},
+		info_f.get_dims(),
+		{{1u,1u}}
+	);
+	
+	// Inflow
+	iterator<Desc::D>::apply(
+		[&](auto& index){
+			info_f.at(index).set(3u);
+		},
+		{{0u,0u}},
+		{{1u,dim_y}},
+		{{0u,1u}}
+	);
+	
+	// Outflow
+	iterator<Desc::D>::apply(
+		[&](auto& index){
+			info_f.at(index).set(4u);
+		},
+		{{dim_x-1u,0u}},
+		{{dim_x, dim_y}},
+		{{0u,1u}}
+	);
+	//
+	auto& df_f = fields.template get<"dfs_old">();
+	auto& rho_f = macros.template get<"density">();
+	auto& vel_f = macros.template get<"momentum">();
+	auto& por_f = macros.template get<"porosity">();
+	
+	iterator<Desc::D>::apply(
+		[&](auto& index){
+			auto& df = df_f.at(index);
+			auto& rho = rho_f.at(index);
+			por_f.at(index).at({}) = {1};
+			rho.at({}) = {1};
+			auto& vel = vel_f.at(index);
+			auto eq = equilibrium<T,Desc>(rho,vel);
+
+			df = eq;
+		},
+		{},// 0-index
+		df_f.get_dims()
+	);
+
+	iterator<Desc::D>::apply(
+		[&](auto& index){
+			auto& df = df_f.at(index);
+			auto& rho = rho_f.at(index);
+			rho.at({}) = {1};
+			auto& vel = vel_f.at(index);
+			if(info_f.at(index).get() == 2u){
+				vel.at({{0u}}) = 0.0;
+			}
+			auto eq = equilibrium<T,Desc>(rho,vel);
+
+			df = eq;
+		},
+		{},// 0-index
+		df_f.get_dims(),
+		{{1u,1u}}
+	);
+	
+	iterator<Desc::D>::apply(
+		[&](auto& index){
+			saw::data<sch::Vector<T,Desc::D>> middle, ind_vec;
+			middle.at({{0u}}) = dim_x * 0.25;
+			middle.at({{1u}}) = dim_y * 0.5;
+
+			ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to<T>();
+			ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to<T>();
+
+			auto dist = middle - ind_vec;
+			auto dist_2 = saw::math::dot(dist,dist);
+			if(dist_2.at({}).get() < dim_y*dim_y*0.01){
+				porous_f.at(index).at({}) = 0.0;
+			}
+		},
+		{},// 0-index
+		df_f.get_dims()
+	);
+
+	{
+		saw::data<sch::Scalar<T>> radius_p;
+		radius_p.at({}).set(2);
+		saw::data<sch::Scalar<T>> dense_p;
+		dense_p.at({}).set(1);
+		particles = create_spheroid_particle_group<T,Desc::D,1u>(radius_p, dense_p, {{16u}});
+	}
+
+	return saw::make_void();
+}
+
+}
+}
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp
new file mode 100644
index 0000000..c7f1014
--- /dev/null
+++ b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp
@@ -0,0 +1,226 @@
+#include "common.hpp"
+#include "init.hpp"
+#include "step.hpp"
+
+#include <forstio/codec/args.hpp>
+
+/**
+ * For deciding what parameters to use maybe?
+ */
+namespace args {
+using namespace saw::schema;
+
+using LbmArgsStruct = Struct<
+	Member<String, "coupling">
+>;
+
+using LbmArgs = Args<
+	LbmArgsStruct,
+	Tuple<>
+>;
+}
+
+template<typename T, typename Desc, typename Coll>
+saw::error_or<void> lbm_main(const saw::data<args::LbmArgs>& args){
+	using namespace kel::lbm;
+
+	auto& an = args.template get<"args">();
+
+	using dfi = df_info<T,Desc>;
+
+	auto eo_lbm_dir = output_directory();
+	if(eo_lbm_dir.is_error()){
+		return std::move(eo_lbm_dir.get_error());
+	}
+	auto& lbm_dir = eo_lbm_dir.get_value();
+
+	auto out_dir = lbm_dir / "stdkr_2d_gpu" / an.template get<"coupling">().stl_view();
+
+	{
+		std::error_code ec;
+		std::filesystem::create_directories(out_dir,ec);
+		if(ec != std::errc{}){
+			return saw::make_error<saw::err::critical>("Could not create output directory");
+		}
+	}
+
+	converter<T> conv {
+		// delta_x
+		{{1.0}},
+		// delta_t
+		{{1.0}}
+	};
+	
+	print_lbm_meta<T,Desc>(conv,{0.1},{1e-4},{0.4 * dim_y});
+
+	// saw::data<sch::FixedArray<sch::UInt64,Desc::D>> meta{{dim_x,dim_y}};
+	auto lbm_data_ptr = saw::heap<saw::data<sch::ChunkStruct<T,Desc>>>();
+	auto lbm_macro_data_ptr = saw::heap<saw::data<sch::MacroStruct<T,Desc>>>();
+	auto lbm_particle_data_ptr = saw::heap<saw::data<sch::ParticleSpheroidGroup<T,Desc>>>();
+
+	std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl;
+
+	auto eo_aio = saw::setup_async_io();
+	if(eo_aio.is_error()){
+		return std::move(eo_aio.get_error());
+	}
+	auto& aio = eo_aio.get_value();
+	saw::wait_scope wait{aio.event_loop};
+
+	bool krun = true;
+	bool print_status = false;
+	aio.event_port.on_signal(saw::Signal::Terminate).then([&](){
+		krun = false;
+	}).detach();
+	aio.event_port.on_signal(saw::Signal::User1).then([&](){
+		print_status = true;
+	}).detach();
+
+	device dev;
+
+	auto& sycl_q = dev.get_handle();
+
+	sycl_q.wait();
+	{
+		auto eov = init<T,Desc,Coll>(conv,*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_data_ptr);
+		if(eov.is_error()){
+			return eov;
+		}
+	}
+	{
+		auto eov = write_vtk_file(out_dir,"initial_state",0u,*lbm_data_ptr);
+		if(eov.is_error()){
+			return eov;
+		}
+	}
+
+	saw::data<sch::ChunkStruct<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_data{sycl_q};
+	saw::data<sch::MacroStruct<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_macro_data{sycl_q};
+	saw::data<sch::ParticleSpheroidGroup<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_particle_data{sycl_q};
+	sycl_q.wait();
+
+	{
+		auto eov = dev.copy_to_device(*lbm_data_ptr,lbm_sycl_data);
+		if(eov.is_error()){
+			return eov;
+		}
+	}
+	{
+		auto eov = dev.copy_to_device(*lbm_macro_data_ptr,lbm_sycl_macro_data);
+		if(eov.is_error()){
+			return eov;
+		}
+	}
+	sycl_q.wait();
+	auto lsd_view = make_view(lbm_sycl_data);
+	auto lsdm_view = make_view(lbm_sycl_macro_data);
+	auto lsdp_view = make_view(lbm_sycl_particle_data);
+	
+	saw::data<sch::UInt64> time_steps{64u*1024ul};
+
+	auto& info_f = lsd_view.template get<"info">();
+
+	for(saw::data<sch::UInt64> i{0u}; i < time_steps and krun; ++i){
+		// BC + Collision
+		{
+			auto eov = step<T,Desc,Coll>(conv,lsd_view,lsdm_view,lsdp_view,i,dev);
+			if(eov.is_error()){
+				return eov;
+			}
+		}
+
+		sycl_q.wait();
+		if(i.get() % 64u == 0u){
+			{
+				auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr);
+				if(eov.is_error()){
+					return eov;
+				}
+			}
+			{
+				auto eov = write_csv_file(out_dir,"m",i.get(), *lbm_macro_data_ptr);
+				if(eov.is_error()){
+					return eov;
+				}
+			}
+		}
+		// Stream
+		sycl_q.submit([&](acpp::sycl::handler& h){
+			component<T,Desc,cmpt::Stream,encode::Sycl<saw::encode::Native>> stream;
+
+			h.parallel_for(acpp::sycl::range<Desc::D>{dim_x,dim_y}, [=](acpp::sycl::id<Desc::D> idx){
+				saw::data<sch::FixedArray<sch::UInt64,Desc::D>> index;
+				for(uint64_t i = 0u; i < Desc::D; ++i){
+					index.at({{i}}).set(idx[i]);
+				}
+
+				auto info = info_f.at(index).get();
+
+				if(info > 0u){
+					stream.apply(lsd_view,index,i);
+				}
+			});
+		}).wait();
+		wait.poll();
+
+		if(print_status){
+			std::cout<<"Status: "<<i.get()<<" of "<<time_steps.get()<<" - "<<(i.template cast_to<sch::Float64>().get() * 100 / time_steps.get())<<"%"<<std::endl;
+			print_status = false;
+		}
+		print_progress_bar(i.get(), time_steps.get()-1u);
+	}
+	// After Loop
+	sycl_q.wait();
+	{
+		auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr);
+		if(eov.is_error()){
+			return eov;
+		}
+	}
+	{
+		auto eov = write_vtk_file(out_dir,"m",time_steps.get(), *lbm_macro_data_ptr);
+		if(eov.is_error()){
+			return eov;
+		}
+	}
+
+	sycl_q.wait();
+	return saw::make_void();
+}
+
+using FloatT = kel::lbm::sch::Float32;
+using DescT = kel::lbm::sch::D2Q9;
+
+saw::error_or<void> k_main(int argc, char** argv){
+	using namespace kel::lbm;
+	
+	saw::data<args::LbmArgs> args;
+	{
+		auto eov = saw::parse_args(args,argc,argv);
+		return eov;
+	}
+	auto& an = args.template get<"args">();
+	auto& coupling = an.template get<"coupling">();
+
+	if(coupling.stl_view() == "hlbm"){
+		return lbm_main<FloatT,DescT,method::Hlbm>(args);
+	}
+
+	return saw::make_error<saw::err::critical>("Invalid coupling");
+}
+
+
+int main(int argc, char** argv){
+	auto eov = k_main(argc, argv);
+	if(eov.is_error()){
+		auto& err = eov.get_error();
+		std::cerr<<"[Error] "<<err.get_category();
+		auto err_msg = err.get_message();
+		if(err_msg.size() > 0u){
+			std::cerr<<" - "<<err_msg;
+		}
+		std::cerr<<std::endl;
+		return err.get_id();
+	}
+	return 0;
+}
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/step.hpp b/examples/schaefer_turek_durst_krause_rannbacher/step.hpp
new file mode 100644
index 0000000..37eb4e7
--- /dev/null
+++ b/examples/schaefer_turek_durst_krause_rannbacher/step.hpp
@@ -0,0 +1,97 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+
+template<typename T, typename Desc, typename Coll>
+saw::error_or<void> step(
+		const converter<T>& conv,
+		saw::data<sch::Ptr<sch::ChunkStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& fields,
+		saw::data<sch::Ptr<sch::MacroStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& macros,
+		saw::data<sch::Ptr<sch::ParticleSpheroidGroup<T,Desc>>,encode::Sycl<saw::encode::Native>>& particles,
+		saw::data<sch::UInt64> t_i,
+		device& dev
+){
+	auto& q = dev.get_handle();
+	auto& info_f = fields.template get<"info">();
+	auto& porous_f = macros.template get<"porosity">();
+	
+	q.submit([&](acpp::sycl::handler& h){
+		component<T,Desc,cmpt::HlbmReset,encode::Sycl<saw::encode::Native>> hlbm_reset;
+		h.parallel_for(acpp::sycl::range<Desc::D>{dim_x,dim_y}, [=](acpp::sycl::id<Desc::D> idx){
+			saw::data<sch::FixedArray<sch::UInt64,Desc::D>> index;
+			for(uint64_t i = 0u; i < Desc::D; ++i){
+				index.at({{i}}).set(idx[i]);
+			}
+
+			hlbm_reset.apply(fields,macros,index,t_i);
+		});
+	}).wait();
+
+	q.submit([&](acpp::sycl::handler& h){
+		component<T,Desc,cmpt::HlbmOneParticle,encode::Sycl<saw::encode::Native>> hlbm_one_part;
+
+		h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){
+			saw::data<sch::FixedArray<sch::UInt64,1u>> index;
+			for(uint64_t i = 0u; i < 1u; ++i){
+				index.at({{i}}).set(idx[i]);
+			}
+
+			hlbm_one_part.apply(fields,macros,particles,index,t_i,{16u});
+		});
+	}).wait();
+
+	// auto coll_ev = 
+	q.submit([&](acpp::sycl::handler& h){
+		component<T,Desc,cmpt::Hlbm,encode::Sycl<saw::encode::Native>> collision{0.8};
+		component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb;
+
+		component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> flow_in{
+			[&](){
+				uint64_t target_t_i = 16u;
+				if(t_i.get() < target_t_i){
+					return 1.0 + (0.01 / target_t_i) * t_i.get();
+				}
+				return 1.01;
+			}()
+		};
+		component<T,Desc,cmpt::ZouHeHorizontal<false>,encode::Sycl<saw::encode::Native>> flow_out{1.0};
+
+		h.parallel_for(acpp::sycl::range<Desc::D>{dim_x,dim_y}, [=](acpp::sycl::id<Desc::D> idx){
+			saw::data<sch::FixedArray<sch::UInt64,Desc::D>> index;
+			for(uint64_t i = 0u; i < Desc::D; ++i){
+				index.at({{i}}).set(idx[i]);
+			}
+
+			auto info = info_f.at(index);
+			
+			switch(info.get()){
+				case 0u:
+				break;
+				case 1u:
+					bb.apply(fields,index,t_i);
+				break;
+				case 2u:
+					collision.apply(fields,macros,index,t_i);
+					break;
+				case 3u:
+					flow_in.apply(fields,index,t_i);
+					collision.apply(fields,macros,index,t_i);
+					break;
+				case 4u:
+					flow_out.apply(fields,index,t_i);
+					collision.apply(fields,macros,index,t_i);
+				break;
+				default:
+				break;
+			}
+		});
+	}).wait();
+
+	return saw::make_void();
+}
+
+}
+}
diff --git a/examples/settling_cubes_2d_ibm_gpu/sim.cpp b/examples/settling_cubes_2d_ibm_gpu/sim.cpp
index 33712b5..d7c51bc 100644
--- a/examples/settling_cubes_2d_ibm_gpu/sim.cpp
+++ b/examples/settling_cubes_2d_ibm_gpu/sim.cpp
@@ -45,7 +45,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<VectorChunk<T,Desc>, "force">
 >;
@@ -93,7 +93,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template gStarted hearing about similar cases not long after, many of which were successful on the part of the crooks. Scary stuff. et<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	
 	iterator<Desc::D>::apply(
 		[&](auto& index){
@@ -163,7 +163,7 @@ saw::error_or<void> step(
 
 	auto& parts = p_group.template get<"particles">();
 	auto& p_mask = p_group.template get<"mask">();
-	auto& vels = macros.template get<"velocity">();
+	auto& vels = macros.template get<"momentum">();
 	auto& forces = macros.template get<"force">();
 
 	auto p_meta = parts.meta();
@@ -214,7 +214,7 @@ saw::error_or<void> step(
 					// Lagrange indicator position
 					auto p_pos_lag = p_pos + transformed_pos;
 
-					// Pick the closest velocity and clamp it
+					// Pick the closest momentum and clamp it
 					saw::data<sch::FixedArray<sch::UInt64,Desc::D>> p_cell_pos;
 					saw::data<sch::Vector<sch::UInt64,Desc::D>> p_cell_pos_vec;
 					for(uint64_t i{0u}; i < Desc::D; ++i){
diff --git a/examples/settling_spheres_2d_hlbm_gpu/sim.cpp b/examples/settling_spheres_2d_hlbm_gpu/sim.cpp
index f7e47b6..5913a1f 100644
--- a/examples/settling_spheres_2d_hlbm_gpu/sim.cpp
+++ b/examples/settling_spheres_2d_hlbm_gpu/sim.cpp
@@ -49,7 +49,7 @@ using PorChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<PorChunk<T>, "porosity">
 >;
@@ -91,7 +91,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template gStarted hearing about similar cases not long after, many of which were successful on the part of the crooks. Scary stuff. et<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	
 	iterator<Desc::D>::apply(
 		[&](auto& index){
@@ -168,7 +168,7 @@ saw::error_or<void> step(
 
 	auto& parts = p_group.template get<"particles">();
 	auto& p_mask = p_group.template get<"mask">();
-	auto& vels = macros.template get<"velocity">();
+	auto& vels = macros.template get<"momentum">();
 	auto& forces = macros.template get<"force">();
 
 	auto p_meta = parts.meta();
diff --git a/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp b/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp
index 9d2ec12..bbe8cfd 100644
--- a/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp
+++ b/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp
@@ -49,7 +49,7 @@ using ForceChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">
 >;
 
@@ -122,7 +122,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 
 	iterator<Desc::D>::apply(
 		[&](auto& index){
diff --git a/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp
index 9e96219..9cf974a 100644
--- a/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp
+++ b/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp
@@ -49,7 +49,7 @@ using ForceChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<VelChunk<T,Desc>, "force">,
 	Member<ScalarChunk<T,Desc>, "porosity">
@@ -125,7 +125,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 
 	iterator<Desc::D>::apply(
diff --git a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp
index 9c28d88..b63a003 100644
--- a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp
+++ b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp
@@ -46,7 +46,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<VelChunk<T,Desc>, "force">,
 	Member<ScalarChunk<T,Desc>, "porosity">
@@ -122,7 +122,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 	
 	iterator<Desc::D>::apply(
diff --git a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp
index 75cbe54..6b86c68 100644
--- a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp
+++ b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp
@@ -47,7 +47,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<VelChunk<T,Desc>, "force">,
 	Member<ScalarChunk<T,Desc>, "porosity">
@@ -124,7 +124,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& vel_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 	
 	iterator<Desc::D>::apply(
diff --git a/modules/core/c++/boundary.hpp b/modules/core/c++/boundary.hpp
index 7f185ad..e7567c9 100644
--- a/modules/core/c++/boundary.hpp
+++ b/modules/core/c++/boundary.hpp
@@ -133,14 +133,14 @@ template<typename FP, typename Descriptor, typename Encode>
 class component<FP, Descriptor, cmpt::Equilibrium, Encode> final {
 private:
 	saw::data<sch::Scalar<FP>> density_;
-	saw::data<sch::Vector<FP,Descriptor::D>> velocity_;
+	saw::data<sch::Vector<FP,Descriptor::D>> momentum_;
 public:
 	component(
 		saw::data<sch::Scalar<FP>> density__,
-		saw::data<sch::Vector<FP,Descriptor::D>> velocity__
+		saw::data<sch::Vector<FP,Descriptor::D>> momentum__
 	):
 		density_{density__},
-		velocity_{velocity__}
+		momentum_{momentum__}
 	{}
 
 	template<typename CellFieldSchema>
@@ -151,7 +151,7 @@ public:
 		
 		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
 
-		auto eq = equilibrium<FP,Descriptor>(density_,velocity_);
+		auto eq = equilibrium<FP,Descriptor>(density_,momentum_);
 
 		dfs_old_f.at(index) = eq;
 	}
@@ -726,12 +726,12 @@ public:
 template<typename T, typename Descriptor, bool East, typename Encode>
 class component<T, Descriptor, cmpt::ZouHeVelocityX<East>, Encode> final {
 private:
-	saw::data<sch::Vector<T,Descriptor::D>> velocity_;
+	saw::data<sch::Vector<T,Descriptor::D>> momentum_;
 public:
 	component(
-		saw::data<sch::Vector<T,Descriptor::D>> velocity__
+		saw::data<sch::Vector<T,Descriptor::D>> momentum__
 	):
-		velocity_{velocity__}
+		momentum_{momentum__}
 	{}
 
 	template<typename CellFieldSchema>
@@ -753,20 +753,20 @@ public:
 		}
 
 		saw::data<sch::Scalar<T>> rho;
-		rho.at({}) = (dfs.at({0u}) + dfs.at({4u}) + dfs.at({3u}) + dir_sum * 2) / (velocity_.at({{0u}}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0)});
+		rho.at({}) = (dfs.at({0u}) + dfs.at({4u}) + dfs.at({3u}) + dir_sum * 2) / (momentum_.at({{0u}}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0)});
 
 		if constexpr (East) {
-			dfs.at({2u}) = dfs.at({1u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(2.0 / 3.0)} * rho.at({}) * velocity_.at({{0u}});
-			dfs.at({6u}) = dfs.at({5u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}})
-				+ (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * velocity_.at({{1u}})) * 0.5f;
-			dfs.at({8u}) = dfs.at({7u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}})
-				+ (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * velocity_.at({{1u}})) * 0.5f;
+			dfs.at({2u}) = dfs.at({1u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(2.0 / 3.0)} * rho.at({}) * momentum_.at({{0u}});
+			dfs.at({6u}) = dfs.at({5u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}})
+				+ (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * momentum_.at({{1u}})) * 0.5f;
+			dfs.at({8u}) = dfs.at({7u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}})
+				+ (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * momentum_.at({{1u}})) * 0.5f;
 		}else if constexpr ( not East ){
-			dfs.at({1u}) = dfs.at({2u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(2.0 / 3.0)} * rho.at({}) * velocity_.at({{0u}});
-			dfs.at({5u}) = dfs.at({6u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}})
-				+ (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * velocity_.at({{1u}}))*0.5f;
-			dfs.at({7u}) = dfs.at({8u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}})
-				+ (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * velocity_.at({{1u}}))*0.5f;
+			dfs.at({1u}) = dfs.at({2u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(2.0 / 3.0)} * rho.at({}) * momentum_.at({{0u}});
+			dfs.at({5u}) = dfs.at({6u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}})
+				+ (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * momentum_.at({{1u}}))*0.5f;
+			dfs.at({7u}) = dfs.at({8u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}})
+				+ (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * momentum_.at({{1u}}))*0.5f;
 		}
 	}
 };
diff --git a/modules/core/c++/collision.hpp b/modules/core/c++/collision.hpp
index 023f61f..39514b8 100644
--- a/modules/core/c++/collision.hpp
+++ b/modules/core/c++/collision.hpp
@@ -75,7 +75,7 @@ public:
 		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
 
 		auto& rho_f = macros.template get<"density">();
-		auto& vel_f = macros.template get<"velocity">();
+		auto& vel_f = macros.template get<"momentum">();
 
 		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
 		saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);
@@ -135,7 +135,7 @@ public:
 		auto& dfs = dfs_old_f.at(index);
 
 		auto& rho_f = macros.template get<"density">();
-		auto& vel_f = macros.template get<"velocity">();
+		auto& vel_f = macros.template get<"momentum">();
 	
 		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
 		
diff --git a/modules/core/c++/component/psm.hpp b/modules/core/c++/component/psm.hpp
index cada53c..5832760 100644
--- a/modules/core/c++/component/psm.hpp
+++ b/modules/core/c++/component/psm.hpp
@@ -39,7 +39,7 @@ public:
 		auto& porous_f = macros.template get<"porosity">();
 
 		auto& rho_f = macros.template get<"density">();
-		auto& vel_f = macros.template get<"velocity">();
+		auto& vel_f = macros.template get<"momentum">();
 		
 		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
 		saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);
diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp
index e2527ce..382d067 100644
--- a/modules/core/c++/fplbm.hpp
+++ b/modules/core/c++/fplbm.hpp
@@ -46,7 +46,7 @@ public:
 
 			auto& rho_f = macros.template get<"density">();
 			auto& rho = rho_f.at(index);
-			auto& vel_f = macros.template get<"velocity">();
+			auto& vel_f = macros.template get<"momentum">();
 			auto& vel = vel_f.at(index);
 
 			compute_rho_u<T,Descriptor>(dfs,rho,vel);
@@ -86,7 +86,7 @@ public:
 		auto& rho_f = macros.template get<"density">();
 		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
 
-		auto& vel_f = macros.template get<"velocity">();
+		auto& vel_f = macros.template get<"momentum">();
 		saw::data<sch::Vector<T,Descriptor::D>> vel = vel_f.at(index);
 
 		auto& force_f = macros.template get<"force">();
@@ -148,7 +148,7 @@ public:
 		// auto& dfs = dfs_old_f.at(index);
 
 		auto& rho_f = macros.template get<"density">();
-		auto& vel_f = macros.template get<"velocity">();
+		auto& vel_f = macros.template get<"momentum">();
 		auto& por_f = macros.template get<"porosity">();
 		auto& force_f = macros.template get<"force">();
 
@@ -245,7 +245,7 @@ public:
 		auto& dfs = dfs_old_f.at(index);
 
 		auto& rho_f = macros.template get<"density">();
-		auto& vel_f = macros.template get<"velocity">();
+		auto& vel_f = macros.template get<"momentum">();
 		auto& por_f = macros.template get<"porosity">();
 
 		// Temporary find a better way
diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp
index cc21dcf..e8ea66a 100644
--- a/modules/core/c++/hlbm.hpp
+++ b/modules/core/c++/hlbm.hpp
@@ -70,7 +70,7 @@ public:
 
 		auto& porosity_f = macros.template get<"porosity">();
 		auto& rho_f = macros.template get<"density">();
-		auto& vel_f = macros.template get<"velocity">();
+		auto& vel_f = macros.template get<"momentum">();
 
 		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
 		saw::data<sch::Vector<T,Desc::D>>& vel = vel_f.at(index);
@@ -125,7 +125,7 @@ public:
 		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
 
 		auto& part_spheroid_group = part_group;
-		auto& mvel_f = macros.template get<"velocity">();
+		auto& mvel_f = macros.template get<"momentum">();
 		auto& mrho_f = macros.template get<"density">();
 		auto& mpor_f = macros.template get<"porosity">();
 		
diff --git a/modules/core/c++/macroscopic.hpp b/modules/core/c++/macroscopic.hpp
index 19ab3e9..768fef0 100644
--- a/modules/core/c++/macroscopic.hpp
+++ b/modules/core/c++/macroscopic.hpp
@@ -28,6 +28,7 @@ void compute_rho_u (
 		}
 	}
 
+	// TODO - Change to momentum, but that needs checks in every impacted function
 	for(size_t i = 0; i < Desc::D; ++i){
 		vel().at({{i}}) = vel().at({{i}}) / rho().at({});
 	}
diff --git a/modules/core/c++/particle.hpp b/modules/core/c++/particle.hpp
index 2c60884..efcc3fb 100644
--- a/modules/core/c++/particle.hpp
+++ b/modules/core/c++/particle.hpp
@@ -25,7 +25,7 @@ public:
 
 		// Compute forces
 
-		// Update particle velocity
+		// Update particle momentum
 		verlet_step_lambda_old<T,Descriptor::D>(p,{1.0});  
 
 		// Update porosity over lattice nodes
diff --git a/modules/core/c++/psm.hpp b/modules/core/c++/psm.hpp
index 15ad0f8..2daff5e 100644
--- a/modules/core/c++/psm.hpp
+++ b/modules/core/c++/psm.hpp
@@ -55,7 +55,7 @@ public:
 		auto& porous_f = macros.template get<"porosity">();
 
 		auto& rho_f = macros.template get<"density">();
-		auto& vel_f = macros.template get<"velocity">();
+		auto& vel_f = macros.template get<"momentum">();
 		
 		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
 		saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);
@@ -105,116 +105,418 @@ public:
 	component() = default;
 
 	template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema>
-	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& pg, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step, saw::data<sch::UInt64> sub_steps) const {
-	
-		using dfi = df_info<T,Descriptor>;
-		bool is_even = ((time_step.get() % 2) == 0);
-		
-		saw::data<sch::Scalar<T>> one;
-		one.at({}) = 1.0;
+void apply(
+    const saw::data<CellFieldSchema, Encode>& field,
+    const saw::data<MacroFieldSchema, Encode>& macros,
+    const saw::data<ParticleSchema, Encode>& pg,
+    saw::data<sch::FixedArray<sch::UInt64,1u>> index,
+    saw::data<sch::UInt64> time_step,
+    saw::data<sch::UInt64> sub_steps
+) const {
 
-		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+    using dfi = df_info<T,Descriptor>;
 
-		auto& por_f = macros.template get<"porosity">();
-		auto& rho_f = macros.template get<"density">();
-		auto& vel_f = macros.template get<"velocity">();
-		auto& force_f = macros.template get<"force">();
-	
-		{
-			auto parts = pg.template get<"particles">();
-			auto parts_size = parts.meta().at({0u});
-
-			auto& p_coll = pg.template get<"collision">().at({});
-			auto& p_rad = p_coll.template get<"radius">();
-
-			auto& pi = parts.at(index);
-			auto& pirb = pi.template get<"rigid_body">();
-			saw::data<sch::Vector<T,Descriptor::D>>& pirb_pos = pirb.template get<"position">();
-			auto& pirb_pos_old = pirb.template get<"position_old">();
-				
-			saw::data<sch::Scalar<T>> ts;
-			ts.at({}) = one.at({}) / sub_steps.template cast_to<T>();
-
-			saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> start;
-			saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> stop;
-
-			auto eo_aabb = particle_aabb<typename ParticleSchema::ValueType>::calculate(pg,{{0u}},vel_f.meta());
-			if(eo_aabb.is_error()){
-				return;
-			}
-			auto& aabb = eo_aabb.get_value();
+    bool is_even = ((time_step.get() % 2u) == 0u);
 
-			/// Ok, I iterate over the space which covers our particle? So lower bounds to upper bounds
-			start = aabb.template get<"a">();
-			stop = aabb.template get<"b">();
+    /*
+     * ------------------------------------------------------------
+     * Fields
+     * ------------------------------------------------------------
+     */
 
-			saw::data<sch::Vector<T,Descriptor::D>> force_p;
-			for(uint64_t i{0u}; i < Descriptor::D; ++i){
-				force_p.at({{i}}) = 0.0;
-			}
-			auto vel_p_old = (pirb_pos-pirb_pos_old) / ts;
-
-			iterator<Descriptor::D>::apply([&](const auto& index_f) -> void{
-				// ask for the d_k value here.
-				// For every value im iterating over I need sth
-				// std::cout<<"Pos: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl;
-			
-				auto& dfs = dfs_old_f.at(index_f);
-				saw::data<sch::Vector<T,Descriptor::D>> rel_dist = saw::math::vectorize_data(index_f).template cast_to<T>() - pirb_pos;
-				saw::data<sch::Scalar<T>> eps;
-				eps.at({}) = 1.5f;
-				
-				auto& por = por_f.at(index_f);
-				por = particle_porosity<T,Descriptor::D,1u,por::ParticleSpheroid<T>>::calculate(rel_dist,p_rad,eps);
-
-				if(por.at({}).get() >= 1.0f){
-					return;
-				}
-
-				saw::data<sch::Vector<T,Descriptor::D>> momentum;
-				for(uint64_t i{0u}; i < Descriptor::D; ++i){
-					momentum.at({{i}}) = 0.0;
-				}
-
-				for(uint64_t i{0u}; i < Descriptor::Q; ++i){
-					saw::data<sch::Vector<T,Descriptor::D>> e_i;
-					saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> n_ind_i;
-					for(uint64_t k{0u}; k < Descriptor::D; ++k){
-						e_i.at({{k}}) = (dfi::directions[i])[k];
-						n_ind_i.at({k}) = index_f.at({k}) + (dfi::directions[i])[k];
-					}
-
-					uint64_t i_opp = dfi::opposite_index[i];
-
-					auto u_p_e = saw::math::dot(e_i,vel_p_old);
-
-					saw::data<T> dfs_added = dfs.at({i})*(saw::data<T>{1}-u_p_e.at({})) + dfs_old_f.at(n_ind_i).at({i_opp})*(saw::data<T>{1}+u_p_e.at({}));
-					saw::data<sch::Scalar<T>> dfs_added_v;
-					dfs_added_v.at({}) = dfs_added;
-					auto ei_dfs = e_i * dfs_added_v;
-
-					momentum = momentum + ei_dfs;
-				}
-				// technically needs to adjust for rotation as well
-				
-				auto& force = force_f.at(index_f);
-				auto& rho = rho_f.at(index_f);
-				// To Fluid
-				
-				auto flip_por = one - por;
-				force = momentum * flip_por;
-				// To Particle
-				force_p = force_p - force;
-			},start,stop);
-			
-			auto& pirb_acc = pirb.template get<"acceleration">();
-			pirb_acc = force_p / pg.template get<"total_mass">().at({});
-
-			for(saw::data<sch::UInt64> i{0u}; i < sub_steps; ++i){
-				verlet_step_lambda<T,Descriptor::D>(pi,ts);
-			}
-		}
-	}
+    auto& dfs_old_f =
+        (is_even)
+            ? field.template get<"dfs_old">()
+            : field.template get<"dfs">();
+
+    auto& por_f   = macros.template get<"porosity">();
+    auto& rho_f   = macros.template get<"density">();
+    auto& vel_f   = macros.template get<"momentum">();
+    auto& force_f = macros.template get<"force">();
+
+    saw::data<sch::Scalar<T>> one;
+    one.at({}) = 1.0;
+
+		saw::data<sch::Scalar<T>> half;
+		half.at({}) = 0.5;
+
+    /*
+     * ------------------------------------------------------------
+     * Particle
+     * ------------------------------------------------------------
+     */
+
+    {
+        auto parts = pg.template get<"particles">();
+
+        auto& p_coll =
+            pg.template get<"collision">().at({});
+
+        auto& p_rad =
+            p_coll.template get<"radius">();
+
+        auto& pi = parts.at(index);
+
+        auto& pirb =
+            pi.template get<"rigid_body">();
+
+        saw::data<sch::Vector<T,Descriptor::D>>& pirb_pos =
+            pirb.template get<"position">();
+
+        auto& pirb_pos_old =
+            pirb.template get<"position_old">();
+
+
+        /*
+         * --------------------------------------------------------
+         * Particle timestep
+         * --------------------------------------------------------
+         */
+
+        saw::data<sch::Scalar<T>> ts;
+
+        ts.at({}) =
+            one.at({})
+            / sub_steps.template cast_to<T>();
+
+
+        /*
+         * --------------------------------------------------------
+         * Particle momentum
+         * --------------------------------------------------------
+         */
+
+        auto vel_p_old =
+            (pirb_pos - pirb_pos_old) / ts;
+
+
+        /*
+         * --------------------------------------------------------
+         * Find particle AABB
+         * --------------------------------------------------------
+         */
+
+        auto eo_aabb =
+            particle_aabb<
+                typename ParticleSchema::ValueType
+            >::calculate(
+                pg,
+                {{0u}},
+                vel_f.meta()
+            );
+
+        if(eo_aabb.is_error()){
+            return;
+        }
+
+        auto& aabb = eo_aabb.get_value();
+
+        saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> start;
+        saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> stop;
+
+        start = aabb.template get<"a">();
+        stop  = aabb.template get<"b">();
+
+
+        /*
+         * --------------------------------------------------------
+         * Total force acting on particle
+         * --------------------------------------------------------
+         */
+
+        saw::data<sch::Vector<T,Descriptor::D>> force_p;
+
+        for(uint64_t k = 0u; k < Descriptor::D; ++k){
+            force_p.at({{k}}) = 0.0;
+        }
+
+
+        /*
+         * --------------------------------------------------------
+         * Iterate over particle AABB
+         * --------------------------------------------------------
+         */
+
+        iterator<Descriptor::D>::apply(
+            [&](const auto& index_f) -> void {
+
+                /*
+                 * ------------------------------------------------
+                 * Distribution function at current cell
+                 * ------------------------------------------------
+                 */
+
+                auto& dfs =
+                    dfs_old_f.at(index_f);
+
+
+                /*
+                 * ------------------------------------------------
+                 * Particle relative position
+                 * ------------------------------------------------
+                 */
+
+                saw::data<sch::Vector<T,Descriptor::D>> rel_dist =
+                    saw::math::vectorize_data(index_f)
+                    .template cast_to<T>()
+                    - pirb_pos;
+
+
+                /*
+                 * ------------------------------------------------
+                 * PSM particle porosity
+                 * ------------------------------------------------
+                 */
+
+                saw::data<sch::Scalar<T>> eps;
+                eps.at({}) = 1.5;
+
+                auto& por =
+                    por_f.at(index_f);
+
+                por =
+                    particle_porosity<
+                        T,
+                        Descriptor::D,
+                        1u,
+                        por::ParticleSpheroid<T>
+                    >::calculate(
+                        rel_dist,
+                        p_rad,
+                        eps
+                    );
+
+
+                /*
+                 * Outside particle.
+                 */
+
+                if(por.at({}).get() >= 1.0){
+                    return;
+                }
+
+
+                /*
+                 * ------------------------------------------------
+                 * Momentum exchange for this cell
+                 * ------------------------------------------------
+                 */
+
+                saw::data<sch::Vector<T,Descriptor::D>> momentum;
+
+                for(uint64_t k = 0u; k < Descriptor::D; ++k){
+                    momentum.at({{k}}) = 0.0;
+                }
+
+
+                /*
+                 * ------------------------------------------------
+                 * D3Q27 / general Q momentum exchange
+                 *
+                 * We exclude i=0 because the rest population has
+                 * no momentum.
+                 *
+                 * The original code summed both directions of
+                 * every opposite pair. Therefore we apply 1/2
+                 * after the sum.
+                 * ------------------------------------------------
+                 */
+
+                for(uint64_t i = 1u; i < Descriptor::Q; ++i){
+
+                    uint64_t i_opp =
+                        dfi::opposite_index[i];
+
+
+                    /*
+                     * Direction vector e_i.
+                     */
+
+                    saw::data<sch::Vector<T,Descriptor::D>> e_i;
+
+                    saw::data<
+                        sch::FixedArray<
+                            sch::UInt64,
+                            Descriptor::D
+                        >
+                    > n_ind_i;
+
+
+                    for(uint64_t k = 0u;
+                        k < Descriptor::D;
+                        ++k){
+
+                        e_i.at({{k}}) =
+                            dfi::directions[i][k];
+
+                        /*
+                         * Neighbor in direction e_i.
+                         */
+
+                        n_ind_i.at({k}) =
+                            index_f.at({k})
+                            + dfi::directions[i][k];
+                    }
+
+
+                    /*
+                     * ------------------------------------------------
+                     * Particle momentum projected onto e_i
+                     * ------------------------------------------------
+                     */
+
+                    auto u_p_e =
+                        saw::math::dot(
+                            e_i,
+                            vel_p_old
+                        );
+
+
+                    /*
+                     * ------------------------------------------------
+                     * Momentum-exchange population
+                     *
+                     * f_i(x)
+                     *
+                     * + f_-i(x + e_i)
+                     *
+                     * with moving-wall correction.
+                     * ------------------------------------------------
+                     */
+
+                    saw::data<T> dfs_added =
+                          dfs.at({i})
+                            * (
+                                saw::data<T>{1.0}
+                                - u_p_e.at({})
+                            )
+
+                        + dfs_old_f
+                            .at(n_ind_i)
+                            .at({i_opp})
+                            * (
+                                saw::data<T>{1.0}
+                                + u_p_e.at({})
+                            );
+
+
+                    /*
+                     * Convert scalar into vector multiplication.
+                     */
+
+                    saw::data<sch::Scalar<T>> dfs_added_v;
+
+                    dfs_added_v.at({}) =
+                        dfs_added;
+
+
+                    /*
+                     * ------------------------------------------------
+                     * Delta p_i
+                     *
+                     * ------------------------------------------------
+                     */
+
+                    momentum =
+                        momentum
+                        + e_i * dfs_added_v;
+                }
+
+
+                /*
+                 * ------------------------------------------------
+                 * Since both i and i_opp were included above,
+                 * compensate for double counting.
+                 *
+                 * This is the important change from your original
+                 * implementation.
+                 * ------------------------------------------------
+                 */
+
+
+
+                momentum =
+                    momentum * half;
+
+
+                /*
+                 * ------------------------------------------------
+                 * PSM solid fraction
+                 *
+                 * por = 1   -> fluid
+                 * por = 0   -> particle
+                 *
+                 * Therefore:
+                 *
+                 * flip_por = 1 - por
+                 * ------------------------------------------------
+                 */
+
+                auto flip_por =
+                    one - por;
+
+
+                /*
+                 * ------------------------------------------------
+                 * Force transferred TO FLUID
+                 * ------------------------------------------------
+                 */
+
+                auto& force =
+                    force_f.at(index_f);
+
+                force =
+                    momentum * flip_por;
+
+
+                /*
+                 * ------------------------------------------------
+                 * Force transferred TO PARTICLE
+                 *
+                 * Newton's third law.
+                 * ------------------------------------------------
+                 */
+
+                force_p =
+                    force_p - force;
+
+            },
+            start,
+            stop
+        );
+
+
+        /*
+         * ------------------------------------------------------------
+         * Particle acceleration
+         * ------------------------------------------------------------
+         */
+
+        auto& pirb_acc =
+            pirb.template get<"acceleration">();
+
+        pirb_acc =
+            force_p
+            / pg.template get<"total_mass">().at({});
+
+
+        /*
+         * ------------------------------------------------------------
+         * Integrate particle
+         * ------------------------------------------------------------
+         */
+
+        for(saw::data<sch::UInt64> i{0u};
+            i < sub_steps;
+            ++i){
+
+            verlet_step_lambda<T,Descriptor::D>(
+                pi,
+                ts
+            );
+        }
+    }
+}
 };
 
 }
diff --git a/modules/core/c++/write_csv.hpp b/modules/core/c++/write_csv.hpp
index a60e208..aa28613 100644
--- a/modules/core/c++/write_csv.hpp
+++ b/modules/core/c++/write_csv.hpp
@@ -91,7 +91,7 @@ struct lbm_csv_writer<sch::Chunk<T,Ghost,D...>> {
 template<typename T, uint64_t D>
 struct lbm_csv_writer<sch::Vector<T,D>> {
 	static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::Vector<T,D>>& field){
-		static_assert(D > 0, "Non-dimensionality is bad for velocity.");
+		static_assert(D > 0, "Non-dimensionality is bad for momentum.");
 
 		// csv_file<<"VECTORS "<<name<<" float\n";
 		for(uint64_t i = 0u; i < D; ++i){
diff --git a/modules/core/c++/write_vtk.hpp b/modules/core/c++/write_vtk.hpp
index e852172..67e8a2f 100644
--- a/modules/core/c++/write_vtk.hpp
+++ b/modules/core/c++/write_vtk.hpp
@@ -115,7 +115,7 @@ struct lbm_vtk_writer<sch::Vector<T,D>> {
 		return saw::make_void();
 	}
 	static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::Vector<T,D>>& field){
-		static_assert(D > 0, "Non-dimensionality is bad for velocity.");
+		static_assert(D > 0, "Non-dimensionality is bad for momentum.");
 		static_assert(D <= 3, "4th dimension as well. Mostly due to vtk.");
 
 		// vtk_file<<"VECTORS "<<name<<" float\n";
diff --git a/modules/core/tests/equilibrium.cpp b/modules/core/tests/equilibrium.cpp
index 20a1f08..7d36600 100644
--- a/modules/core/tests/equilibrium.cpp
+++ b/modules/core/tests/equilibrium.cpp
@@ -20,7 +20,7 @@ void check_equilibrium(){
 	auto eq = lbm::equilibrium<lbm::sch::Float64,Descriptor>(rho,vel);
 
 	for(saw::data<lbm::sch::UInt64> i{0u}; i.get() < Descriptor::Q; ++i){
-		SAW_EXPECT(eq.at(i).get() == dfi::weights[i.get()], std::string{"No velocity and normalized rho should be exactly the weights: "} + std::to_string(eq.at(i).get()) + std::string{" "} + std::to_string(dfi::weights[i.get()]));
+		SAW_EXPECT(eq.at(i).get() == dfi::weights[i.get()], std::string{"No momentum and normalized rho should be exactly the weights: "} + std::to_string(eq.at(i).get()) + std::string{" "} + std::to_string(dfi::weights[i.get()]));
 	}
 }
 
diff --git a/modules/core/tests/vtk_write.cpp b/modules/core/tests/vtk_write.cpp
index c99d752..1de25bf 100644
--- a/modules/core/tests/vtk_write.cpp
+++ b/modules/core/tests/vtk_write.cpp
@@ -31,7 +31,7 @@ using CellStruct = Struct<
 
 template<typename T, uint64_t D>
 using MacroStruct = Struct<
-	Member<FixedArray<T,D>, "velocity">,
+	Member<FixedArray<T,D>, "momentum">,
 	Member<T, "pressure">
 >;
 
@@ -47,7 +47,7 @@ SAW_TEST("VTK Write test example"){
 	saw::data<sch::Array<sch::MacroStruct<sch::T,2>, 2>> cells{{{2u,2u}}};
 
 	auto& cell_0 = cells.at({{{0,0}}});
-	cell_0.template get<"velocity">()= {{0.5,-0.1}};
+	cell_0.template get<"momentum">()= {{0.5,-0.1}};
 	cell_0.template get<"pressure">().set(1.1);
 
 	auto eov = lbm::impl::lbm_vtk_writer<sch::Array<sch::MacroStruct<sch::T,2>, 2>>::apply(sstream, cells);
@@ -56,7 +56,7 @@ SAW_TEST("VTK Write test example"){
 	// I want to print it to see it for myself. For now I have no tooling to more easily view associated and potentially generated files
 	std::cout<<sstream.str()<<std::endl;
 
-	static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS velocity float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n";
+	static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS momentum float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n";
 	SAW_EXPECT(sstream.str() == comparison_str, "Expected different encoding");
 
 	// using Type = typename parameter_pack_type<i,T...>::type;
@@ -74,7 +74,7 @@ SAW_TEST("VTK Write raw test example"){
 	auto meta = cells.dims();
 
 	auto& cell_0 = cells.at({{{0,0}}});
-	cell_0.template get<"velocity">()= {{0.5,-0.1}};
+	cell_0.template get<"momentum">()= {{0.5,-0.1}};
 	cell_0.template get<"pressure">().set(1.1);
 
 	auto eov = lbm::impl::lbm_vtk_writer_raw<sch::MacroStruct<sch::T,2u>,2u>::apply(sstream, &cell_0, meta);
@@ -83,7 +83,7 @@ SAW_TEST("VTK Write raw test example"){
 	// I want to print it to see it for myself. For now I have no tooling to more easily view associated and potentially generated files
 	std::cout<<sstream.str()<<std::endl;
 
-	static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS velocity float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n";
+	static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS momentum float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n";
 	SAW_EXPECT(sstream.str() == comparison_str, "Expected different encoding");
 
 	// using Type = typename parameter_pack_type<i,T...>::type;
diff --git a/poc/triangle_stream/common.hpp b/poc/triangle_stream/common.hpp
index 927f8f8..848d337 100644
--- a/poc/triangle_stream/common.hpp
+++ b/poc/triangle_stream/common.hpp
@@ -37,7 +37,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">
 >;
 }