Merge branch 'dev'

2 files changed, 10 insertions, 5 deletions

diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
index 97ecfde..f89f5da 100644
--- a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
+++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
@@ -201,7 +201,7 @@ saw::error_or<void> step(
 				index.at({{i}}).set(idx[i]);
 			}
 
-			hlbm_one_part.apply(fields,macros,particles,index,t_i,{{128u}});
+			hlbm_one_part.apply(fields,macros,particles,index,t_i,{16u});
 		});
 	}).wait();
 
@@ -380,14 +380,13 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 		}
 		sycl_q.wait();
 
-		if( i.get() % 32u == 33u){
+		if( i.get() % 32u == 0u){
 			{
 				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",i.get(), *lbm_macro_data_ptr);
 				if(eov.is_error()){
diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp
index 4b1f2b4..bb21ae5 100644
--- a/modules/core/c++/hlbm.hpp
+++ b/modules/core/c++/hlbm.hpp
@@ -123,6 +123,7 @@ public:
 
 		auto& part_spheroid_group = part_group;
 		auto& mvel_f = macros.template get<"velocity">();
+		auto& mrho_f = macros.template get<"density">();
 		auto& mpor_f = macros.template get<"porosity">();
 		
 		auto& particle_N_f = field.template get<"particle_N">();
@@ -160,6 +161,7 @@ public:
 			for(uint64_t i{0u}; i < Desc::D; ++i){
 				force_p.at({{i}}) = 0.0;
 			}
+			auto vel_p_old = (pirb_pos-pirb_pos_old) / ts;
 
 			iterator<Desc::D>::apply([&](const auto& index_f) -> void{
 				// ask for the d_k value here.
@@ -193,19 +195,23 @@ public:
 
 					uint64_t i_opp = dfi::opposite_index[i];
 
-					saw::data<T> dfs_added = dfs.at({i}) - dfs_old_f.at(n_ind_i).at({i_opp});
+					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& rho = mrho_f.at(index_f);
 				force_p = force_p + momentum;
 			},start,stop);
 			
 			auto& pirb_acc = pirb.template get<"acceleration">();
-			pirb_acc = force_p / part_spheroid_group.template get<"total_mass">().at({});
+			pirb_acc = - force_p / part_spheroid_group.template get<"total_mass">().at({});
 
 			for(saw::data<sch::UInt64> i{0u}; i < sub_steps; ++i){
 				verlet_step_lambda<T,Desc::D>(pi,ts);