Merge branch 'dev'

3 files changed, 42 insertions, 20 deletions

diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
index 5e8f30a..97ecfde 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);
+			hlbm_one_part.apply(fields,macros,particles,index,t_i,{{128u}});
 		});
 	}).wait();
 
@@ -214,9 +214,9 @@ saw::error_or<void> step(
 			[&](){
 				uint64_t target_t_i = 8u;
 				if(t_i.get() < target_t_i){
-					return 1.0 + (0.0005 / target_t_i) * t_i.get();
+					return 1.0 + (0.005 / target_t_i) * t_i.get();
 				}
-				return 1.0005;
+				return 1.005;
 			}()
 		};
 		component<T,Desc,cmpt::ZouHeHorizontal<false>,encode::Sycl<saw::encode::Native>> flow_out{1.0};
@@ -380,13 +380,14 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 		}
 		sycl_q.wait();
 
-		if( i.get() % 32u == 0u){
+		if( i.get() % 32u == 33u){
 			{
 				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()){
@@ -394,7 +395,16 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 				}
 			}
 		}
-
+		{
+			auto eov = dev.copy_to_host(lbm_sycl_parts_data,*lbm_parts_data_ptr);
+			if(eov.is_error()){
+				return eov;
+			}
+			auto& p = lbm_parts_data_ptr->template get<"particles">().at({0u}).template get<"rigid_body">();
+			auto& ppos = p.template get<"position">();
+			auto& pacc = p.template get<"acceleration">();
+			std::cout<<"Particle: "<<ppos.at({{0u}}).get()<<", "<<ppos.at({{1u}}).get()<<"; "<<pacc.at({{0u}}).get()<<", "<<pacc.at({{1u}}).get()<<";"<<i.get()<<std::endl;
+		}
 		// Stream
 		sycl_q.submit([&](acpp::sycl::handler& h){
 			component<T,Desc,cmpt::Stream,encode::Sycl<saw::encode::Native>> stream;
@@ -417,7 +427,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 			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);
+		// print_progress_bar(i.get(), time_steps.get()-1u);
 	}
 
 	// After Loop
@@ -445,7 +455,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 	return saw::make_void();
 }
 
-using FloatT = kel::lbm::sch::Float32;
+using FloatT = kel::lbm::sch::Float64;
 
 int main(int argc, char** argv){
 	auto eov = lbm_main<FloatT,kel::lbm::sch::D2Q9>(argc, argv);
diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp
index e587b0a..4b1f2b4 100644
--- a/modules/core/c++/hlbm.hpp
+++ b/modules/core/c++/hlbm.hpp
@@ -33,7 +33,9 @@ public:
 		por.at({}) = 1.0;
 
 		auto& pnf = particle_N_f.at(index);
-		pnf = {};
+		for(uint64_t i{0u}; i < Descriptor::D; ++i){
+			pnf.at({{i}}) = 0.0;
+		}
 
 		auto& pnd = particle_D_f.at(index);
 		pnd.at({}) = 0.0;
@@ -105,7 +107,7 @@ private:
 */
 public:
 	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>& part_group, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step) const {
+	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& part_group, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step, saw::data<sch::UInt64> sub_steps) const {
 		/// Figure out how to access the particle list
 		// auto& p = particles.at(i);
 
@@ -134,12 +136,12 @@ public:
 
 			auto& pi = parts.at(index);
 			auto& pirb = pi.template get<"rigid_body">();
-			auto& pirb_pos = pirb.template get<"position">();
+			saw::data<sch::Vector<T,Desc::D>>& pirb_pos = pirb.template get<"position">();
 			auto& pirb_pos_old = pirb.template get<"position_old">();
 				
 			// TODO !!!! Divide by actual time step - for now it's ok
 			saw::data<sch::Scalar<T>> ts;
-			ts.at({}) = 1.0f;
+			ts.at({}) = one.at({}) / sub_steps.template cast_to<T>();
 
 			saw::data<sch::FixedArray<sch::UInt64,Desc::D>> start;
 			saw::data<sch::FixedArray<sch::UInt64,Desc::D>> stop;
@@ -155,25 +157,33 @@ public:
 			stop = aabb.template get<"b">();
 
 			saw::data<sch::Vector<T,Desc::D>> force_p;
-			
+			for(uint64_t i{0u}; i < Desc::D; ++i){
+				force_p.at({{i}}) = 0.0;
+			}
+
 			iterator<Desc::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);
-				auto& mpor = mpor_f.at(index_f);
-				auto rel_dist = saw::math::vectorize_data(index_f).template cast_to<T>() - pirb_pos;
+				saw::data<sch::Vector<T,Desc::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& mpor = mpor_f.at(index_f);
 				mpor = particle_porosity<T,Desc::D,1u,por::ParticleSpheroid<T>>::calculate(rel_dist,p_rad,eps);
-				if(mpor.at({}).get() < 1.0f){
+
+				if(mpor.at({}).get() >= 1.0f){
 					return;
 				}
 
 				saw::data<sch::Vector<T,Desc::D>> momentum;
+				for(uint64_t i{0u}; i < Desc::D; ++i){
+					momentum.at({{i}}) = 0.0;
+				}
 
-				for(uint64_t i = 0u; i < Desc::Q; ++i){
+				for(uint64_t i{0u}; i < Desc::Q; ++i){
 					saw::data<sch::Vector<T,Desc::D>> e_i;
 					saw::data<sch::FixedArray<sch::UInt64,Desc::D>> n_ind_i;
 					for(uint64_t k{0u}; k < Desc::D; ++k){
@@ -191,13 +201,15 @@ public:
 					momentum = momentum + ei_dfs;
 				}
 
-				force_p = force_p + momentum * mpor;
+				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({});
 
-			verlet_step_lambda<T,Desc::D>(pi,ts);
+			for(saw::data<sch::UInt64> i{0u}; i < sub_steps; ++i){
+				verlet_step_lambda<T,Desc::D>(pi,ts);
+			}
 			auto vel_p = (pirb_pos-pirb_pos_old) / ts;
 
 			iterator<Desc::D>::apply([&](const auto& index_f){
diff --git a/modules/core/c++/particle/porosity.hpp b/modules/core/c++/particle/porosity.hpp
index 4d95cce..07e5c9e 100644
--- a/modules/core/c++/particle/porosity.hpp
+++ b/modules/core/c++/particle/porosity.hpp
@@ -38,8 +38,8 @@ public:
 	static saw::data<sch::Scalar<T>> calculate(const saw::data<sch::Vector<T,D>>& lbm_rel_dist, saw::data<sch::Scalar<T>> rad, saw::data<sch::Scalar<T>> eps){
 		saw::data<sch::Scalar<T>> por;
 		
-		auto s_dist_2 = saw::math::dot(lbm_rel_dist,lbm_rel_dist);
-		auto s_dist = saw::math::sqrt(s_dist_2);
+		saw::data<sch::Scalar<T>> s_dist_2 = saw::math::dot(lbm_rel_dist,lbm_rel_dist);
+		saw::data<sch::Scalar<T>> s_dist = saw::math::sqrt(s_dist_2);
 
 		saw::data<sch::Scalar<T>> eps_h;
 		eps_h.at({}) = eps.at({}).get() / 2;