Merge branch 'dev'

9 files changed, 84 insertions, 68 deletions

diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
index daa397a..75dd666 100644
--- a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
+++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
@@ -57,14 +57,14 @@ using MacroStruct = Struct<
 //>;
 
 template<typename T, typename Desc>
-using ParticleSpheroidGroup = ParticleGroup<T,Desc::D,ParticleCollisionSpheroid<T,2.0>>;
+using ParticleSpheroidGroup = ParticleGroup<T,Desc::D,1u,coll::Spheroid<T>>;
 }
 
 template<typename T, typename Desc>
 saw::error_or<void> setup_initial_conditions(
-		saw::data<sch::ChunkStruct<T,Desc>>& fields,
-		saw::data<sch::MacroStruct<T,Desc>>& macros,
-		saw::data<sch::ParticleSpheroidGroup<T,Desc>>& particles
+	saw::data<sch::ChunkStruct<T,Desc>>& fields,
+	saw::data<sch::MacroStruct<T,Desc>>& macros,
+	saw::data<sch::ParticleSpheroidGroup<T,Desc>>& particles
 ){
 	auto& info_f = fields.template get<"info">();
 	auto& porous_f = macros.template get<"porosity">();
@@ -145,6 +145,12 @@ saw::error_or<void> setup_initial_conditions(
 		{{1u,1u}}
 	);
 
+	saw::data<sch::Scalar<T>> rad;
+	rad.at({}) = 2.0f;
+	saw::data<sch::Scalar<T>> dense;
+	dense.at({}) = 1.0f;
+	particles = create_spheroid_particle_group<T,Desc::D,1u>(rad,dense,{16u});
+
 	return saw::make_void();
 }
 
@@ -152,7 +158,7 @@ template<typename T, typename Desc>
 saw::error_or<void> step(
 		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::Ptr<sch::ParticleSpheroidGroup<T,Desc>>,encode::Sycl<saw::encode::Native>>& particles,
 		saw::data<sch::UInt64> t_i,
 		device& dev
 ){
@@ -160,26 +166,35 @@ saw::error_or<void> step(
 	auto& info_f = fields.template get<"info">();
 	auto& porous_f = macros.template get<"porosity">();
 
+	q.submit([&](acpp::sycl::handler& h){
+		component<T,Desc,cmpt::HlbmParticle,encode::Sycl<saw::encode::Native>> hlbm_one_part;
+
+		h.parallel_for(acpp::sycl::range<1u>{1u}, [=](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);
+
+			// auto aabb = particle_aabb<>::cacl();
+
+		});
+	}).wait();
+
 	// auto coll_ev = 
 	q.submit([&](acpp::sycl::handler& h){
 		component<T,Desc,cmpt::Hlbm,encode::Sycl<saw::encode::Native>> collision{0.65};
 		component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb;
 		component<T,Desc,cmpt::AntiBounceBack<0u>,encode::Sycl<saw::encode::Native>> abb;
 
-		saw::data<sch::Scalar<T>> rho_b;
-		rho_b.at({}) = 1.0;
-		saw::data<sch::Vector<T,Desc::D>> vel_b;
-		vel_b.at({{0u}}) = 0.015;
-
-		component<T,Desc,cmpt::Equilibrium,encode::Sycl<saw::encode::Native>> equi{rho_b,vel_b};
-
 		component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> flow_in{
 			[&](){
-				uint64_t target_t_i = 64u;
+				uint64_t target_t_i = 8u;
 				if(t_i.get() < target_t_i){
-					return 1.0 + (0.0002 / target_t_i) * t_i.get();
+					return 1.0 + (0.0005 / target_t_i) * t_i.get();
 				}
-				return 1.0002;
+				return 1.0005;
 			}()
 		};
 		component<T,Desc,cmpt::ZouHeHorizontal<false>,encode::Sycl<saw::encode::Native>> flow_out{1.0};
@@ -265,6 +280,8 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 	// 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_parts_data_ptr = saw::heap<saw::data<sch::ParticleSpheroidGroup<T,Desc>>>();
+	//saw::data<sch::Ptr<sch::ParticleSpheroidGroup<T,Desc>>,encode::Sycl<saw::encode::Native>>& particles
 	
 	std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl;
 
@@ -290,7 +307,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 
 	sycl_q.wait();
 	{
-		auto eov = setup_initial_conditions<T,Desc>(*lbm_data_ptr,*lbm_macro_data_ptr);
+		auto eov = setup_initial_conditions<T,Desc>(*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_parts_data_ptr);
 		if(eov.is_error()){
 			return eov;
 		}
@@ -304,6 +321,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 
 	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_parts_data{sycl_q};
 	sycl_q.wait();
 
 	{
@@ -318,35 +336,30 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 			return eov;
 		}
 	}
+	{
+		auto eov = dev.copy_to_device(*lbm_parts_data_ptr,lbm_sycl_parts_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_parts_data);
+
 	saw::data<sch::UInt64> time_steps{16u*4096ul};
 	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>(lsd_view,lsdm_view,i,dev);
+			auto eov = step<T,Desc>(lsd_view,lsdm_view,lsdp_view,i,dev);
 			if(eov.is_error()){
 				return eov;
 			}
 		}
 		sycl_q.wait();
-		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_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;
@@ -386,6 +399,12 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 			return eov;
 		}
 	}
+	{
+		auto eov = write_csv_file(out_dir,"m",time_steps.get(), *lbm_macro_data_ptr);
+		if(eov.is_error()){
+			return eov;
+		}
+	}
 
 	sycl_q.wait();
 	return saw::make_void();
diff --git a/examples/settling_spheres_2d_hlbm_gpu/sim.cpp b/examples/settling_spheres_2d_hlbm_gpu/sim.cpp
index a004a92..f7e47b6 100644
--- a/examples/settling_spheres_2d_hlbm_gpu/sim.cpp
+++ b/examples/settling_spheres_2d_hlbm_gpu/sim.cpp
@@ -13,7 +13,6 @@ namespace lbm {
 
 constexpr uint64_t dim_x = 512ul;
 constexpr uint64_t dim_y = dim_x * 2ul;
-
 constexpr uint64_t particle_amount = 16ul;
 
 namespace sch {
@@ -59,15 +58,9 @@ template<typename T, typename Desc>
 using ParticleSpheroidGroup = ParticleGroup<
 	T,
 	Desc::D,
-	sch::ParticleCollisionSpheroid<T,2.0f>
->;
-
-template<typename T, typename Desc>
-using ParticleGroups = Tuple<
-	ParticleSpheroidGroup<T,Desc>
+	particle_amount,
+	coll::Spheroid<T>
 >;
-
-
 }
 
 template<typename T, typename Desc>
@@ -186,9 +179,6 @@ saw::error_or<void> step(
 			for(uint64_t i = 0u; i < Desc::D; ++i){
 				index.at({{i}}).set(idx[i]);
 			}
-
-			// Reset the force to zero
-			forces.at(index) = {};
 		});
 	}).wait();
 
diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp
index 9a89697..eccc2cd 100644
--- a/modules/core/c++/hlbm.hpp
+++ b/modules/core/c++/hlbm.hpp
@@ -105,7 +105,7 @@ class component<T, Desc, cmpt::HlbmParticle, Encode> final {
 private:
 	template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema, uint64_t i>
 	void apply_i(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& part_groups, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step) const {
-		
+		// if constexpr ( i < )
 	}	
 public:
 	/*
@@ -125,7 +125,7 @@ public:
 		auto& mvel = macros.template get<"velocity">();
 		{
 			auto& parts = part_spheroid_group.template get<"particles">();
-			auto parts_size = parts.meta().at({0u});
+			auto parts_size = parts.dims().at({0u});
 
 			auto& pi = parts.at(index);
 			auto& pirb = pi.template get<"rigid_body">();
@@ -136,6 +136,8 @@ public:
 
 			auto aabb = particle_aabb<ParticleSchema>::calculate(part_spheroid_group,index,mvel.meta());
 			/// 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">();
 
 			iterator<Desc::D>::apply([&](const auto& index){
 				// ask for the d_k value here.
diff --git a/modules/core/c++/particle/aabb.hpp b/modules/core/c++/particle/aabb.hpp
index 18426d5..983fd97 100644
--- a/modules/core/c++/particle/aabb.hpp
+++ b/modules/core/c++/particle/aabb.hpp
@@ -10,20 +10,24 @@ class particle_aabb final {
 	static_assert(saw::always_false<PGroup>, "Not supported");
 };
 
-template<typename T, uint64_t D>
+template<typename T, uint64_t D, uint64_t PC>
 class particle_aabb<
-	sch::ParticleGroup<T,D,coll::Spheroid<T>>
+	sch::ParticleGroup<T,D,PC,coll::Spheroid<T>>
 > final {
 public:
-	using Schema = sch::ParticleGroup<T,D,coll::Spheroid<T>>;
+	using Schema = sch::ParticleGroup<T,D,PC,coll::Spheroid<T>>;
 
 	using AABB = sch::Struct<
 		sch::Member<sch::FixedArray<sch::UInt64,D>,"a">,
 		sch::Member<sch::FixedArray<sch::UInt64,D>,"b">
 	>;
 public:
-	template<typename Encode>
-	static constexpr saw::data<AABB> calculate(const saw::data<Schema,Encode>& p_grp, const saw::data<sch::FixedArray<sch::UInt64,1u>>& i, const saw::data<sch::FixedArray<sch::UInt64,D>>& meta){
+	template<typename Sch, typename Encode>
+	static constexpr saw::error_or<saw::data<AABB>> calculate(const saw::data<Sch,Encode>& p_grp, const saw::data<sch::FixedArray<sch::UInt64,1u>>& i, const saw::data<sch::FixedArray<sch::UInt64,D>>& meta){
+		static_assert(PC > 0u, "Can't calculate from no particles");
+		if(not (i.at({0u}).get() < PC) ){
+			return saw::make_error<saw::err::critical>("Too large i in particle_aabb");
+		}
 
 		saw::data<AABB> aabb;
 		auto& parts = p_grp.template get<"particles">();
@@ -35,7 +39,7 @@ public:
 		auto& a = aabb.template get<"a">();
 		auto& b = aabb.template get<"b">();
 
-		const saw::data<sch::Scalar<T>>& rad_d = p_grp.template get<"collision">().at({0u}).template get<"radius">();
+		const saw::data<sch::Scalar<T>>& rad_d = p_grp.template get<"collision">().at(i).template get<"radius">();
 
 		saw::data<sch::Vector<T,D>> lower;
 		saw::data<sch::Vector<T,D>> upper;
@@ -51,5 +55,6 @@ public:
 
 	}
 };
+
 }
 }
diff --git a/modules/core/c++/particle/particle.hpp b/modules/core/c++/particle/particle.hpp
index 8e75e5a..3c3630a 100644
--- a/modules/core/c++/particle/particle.hpp
+++ b/modules/core/c++/particle/particle.hpp
@@ -13,13 +13,13 @@
 namespace kel {
 namespace lbm {
 
-template<typename T, uint64_t D>
-saw::data<sch::ParticleGroup<T,D, coll::Spheroid<T>>> create_spheroid_particle_group(
+template<typename T, uint64_t D,uint64_t PartAmount>
+saw::data<sch::ParticleGroup<T, D, PartAmount, coll::Spheroid<T>>> create_spheroid_particle_group(
 		saw::data<sch::Scalar<T>> radius_p,
 		saw::data<sch::Scalar<T>> density_p,
 		const saw::data<sch::UInt64>& mask_resolution
 ){
-	saw::data<sch::ParticleGroup<T,D,coll::Spheroid<T>>> part;
+	saw::data<sch::ParticleGroup<T,D,PartAmount,coll::Spheroid<T>>> part;
 
 	auto& rad_s = part.template get<"collision">().at({0u}).template get<"radius">();
 	rad_s = radius_p;
diff --git a/modules/core/c++/particle/particle_opa.hpp b/modules/core/c++/particle/particle_opa.hpp
index 30a8ade..ce1ecb7 100644
--- a/modules/core/c++/particle/particle_opa.hpp
+++ b/modules/core/c++/particle/particle_opa.hpp
@@ -41,7 +41,7 @@ public:
 		// Write out value
 		auto diff = pos_ind.template cast_to<T>() - pos_;
 		auto diff_dot = saw::math::dot(diff,diff);
-		porous = particle_porosity<T,Descriptor::D,por::ParticleSpheroid<T>>::calculate(diff, rad_, eps_);
+		porous = particle_porosity<T,Descriptor::D,1u,por::ParticleSpheroid<T>>::calculate(diff, rad_, eps_);
 	}
 };
 }
diff --git a/modules/core/c++/particle/porosity.hpp b/modules/core/c++/particle/porosity.hpp
index 2f39d2e..353cd58 100644
--- a/modules/core/c++/particle/porosity.hpp
+++ b/modules/core/c++/particle/porosity.hpp
@@ -11,11 +11,11 @@ template<typename T>
 struct ParticleSpheroid {};
 }
 
-template<typename T, uint64_t D, typename Coll>
+template<typename T, uint64_t D, uint64_t PC, typename Coll>
 class particle_porosity {
 public:
 
-	static saw::data<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,Coll>>& part_group, uint64_t p_i, const saw::data<sch::Vector<T,D>>& lbm_pos){
+	static saw::data<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,PC,Coll>>& part_group, uint64_t p_i, const saw::data<sch::Vector<T,D>>& lbm_pos){
 		auto& mask = part_group.template get<"mask">();
 
 		auto& particles = part_group.template get<"particles">();
@@ -33,7 +33,7 @@ public:
 };
 
 template<typename T, uint64_t D>
-class particle_porosity<T, D, por::ParticleSpheroid<T>> final {
+class particle_porosity<T, D, 1u, por::ParticleSpheroid<T>> final {
 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;
@@ -68,10 +68,10 @@ public:
 };
 
 
-template<typename T, uint64_t D>
-class particle_porosity<T, D, coll::Spheroid<T>> final {
+template<typename T, uint64_t D, uint64_t PC>
+class particle_porosity<T, D, PC, coll::Spheroid<T>> final {
 public:
-	static saw::data<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,coll::Spheroid<T> > >& part_group, uint64_t i, const saw::data<sch::Vector<T,D>>& lbm_pos) {
+	static saw::data<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,PC,coll::Spheroid<T> > >& part_group, uint64_t i, const saw::data<sch::Vector<T,D>>& lbm_pos) {
 		saw::data<sch::Scalar<T>> por;
 
 		auto& parts = part_group.template get<"particles">();
diff --git a/modules/core/c++/particle/schema.hpp b/modules/core/c++/particle/schema.hpp
index 714a16f..221137a 100644
--- a/modules/core/c++/particle/schema.hpp
+++ b/modules/core/c++/particle/schema.hpp
@@ -52,7 +52,7 @@ using Particle = Struct<
 	// Member<Array<Float64,D>, "mask">,
 >;
 
-template<typename T, uint64_t D, typename CollisionType = coll::Spheroid<T>>
+template<typename T, uint64_t D, uint64_t PartAmount, typename CollisionType = coll::Spheroid<T>>
 using ParticleGroup = Struct<
 	Member<Array<T,D>, "mask">,
 	Member<FixedArray<typename CollisionType::Schema,1u>, "collision">,
@@ -61,7 +61,7 @@ using ParticleGroup = Struct<
 	Member<FixedArray<Scalar<T>,1u>, "density">,
 	Member<FixedArray<Vector<T,D>,1u>, "center_of_mass">,
 	Member<FixedArray<Scalar<T>,1u>, "total_mass">,
-	Member<Array<Particle<T,D>,1u>, "particles">
+	Member<FixedArray<Particle<T,D>,PartAmount>, "particles">
 >;
 }
 }
diff --git a/modules/core/tests/particles.cpp b/modules/core/tests/particles.cpp
index d11fa4f..e060097 100644
--- a/modules/core/tests/particles.cpp
+++ b/modules/core/tests/particles.cpp
@@ -278,15 +278,14 @@ SAW_TEST("Spheroid Particle / AABB"){
 
 	using T = sch::Float64;
 	constexpr uint64_t D = 2u;
+	constexpr uint64_t PC = 1u;
 
 	saw::data<sch::Scalar<T>> radi,dense;
 	radi.at({}) = 2.2f;
 	dense.at({}) = 1.0f;
-	auto spheroid_pgrp = lbm::create_spheroid_particle_group<T,D>(radi,dense,{8u});
+	auto spheroid_pgrp = lbm::create_spheroid_particle_group<T,D,PC>(radi,dense,{8u});
 
 	auto& rb_array = spheroid_pgrp.template get<"particles">();
-	rb_array = {1u};
-
 	auto& rb = rb_array.at({0u}).template get<"rigid_body">();
 
 	auto& rb_pos = rb.template get<"position">();
@@ -295,8 +294,9 @@ SAW_TEST("Spheroid Particle / AABB"){
 		rb_pos.at({{1u}}).set(2.9f);
 	}
 
-	auto res = lbm::particle_aabb<lbm::sch::ParticleGroup<T,D,lbm::coll::Spheroid<T>>>::calculate(spheroid_pgrp,{{0u}},{{64u,32u}});
-
+	auto eo_res = lbm::particle_aabb<lbm::sch::ParticleGroup<T,D,PC,lbm::coll::Spheroid<T>>>::calculate(spheroid_pgrp,{{0u}},{{64u,32u}});
+	SAW_EXPECT(eo_res.is_value(), "particle_aabb has error value");
+	auto& res = eo_res.get_value();
 	{
 		auto res_check = [&](const auto& tbp){
 			lbm::iterator<1u>::apply([&](const auto& index){