HLBM kinda broken again?

2 files changed, 13 insertions, 24 deletions

diff --git a/examples/poiseulle_particles_2d_gpu/sim.cpp b/examples/poiseulle_particles_2d_gpu/sim.cpp
index 691ec93..6d5452e 100644
--- a/examples/poiseulle_particles_2d_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_gpu/sim.cpp
@@ -26,11 +26,14 @@ 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<ScalarChunk<T,Desc>, "particle_N">,
+	Member<VectorChunk<T,Desc>, "particle_N">,
 	Member<ScalarChunk<T,Desc>, "particle_D">
 >;
 
@@ -102,11 +105,13 @@ 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& 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);
@@ -324,7 +329,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 		}
 	}
 	sycl_q.wait();
-	saw::data<sch::UInt64> time_steps{particle_size};
+	saw::data<sch::UInt64> time_steps{256ul};
 
 	for(saw::data<sch::UInt64> i{0u}; i < time_steps and krun; ++i){
 		{
@@ -348,15 +353,6 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 			}*/
 		}
 		{
-			std::string file_name = "df_";
-			file_name += std::to_string(i.get());
-			file_name += ".vtk";
-			auto eov = write_vtk_file(out_dir/file_name, *lbm_data_ptr);
-			if(eov.is_error()){
-				return eov;
-			}
-		}
-		{
 			auto eov = step<T,Desc>(lsd_view,lsdm_view,i,dev);
 			if(eov.is_error()){
 				return eov;
@@ -391,15 +387,6 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 			return eov;
 		}
 	}
-	{
-		std::string file_name = "df_";
-		file_name += std::to_string(time_steps.get());
-		file_name += ".vtk";
-		auto eov = write_vtk_file(out_dir/file_name, *lbm_data_ptr);
-		if(eov.is_error()){
-			return eov;
-		}
-	}
 
 	/*
 	iterator<Desc::D>::apply(
diff --git a/lib/core/c++/hlbm.hpp b/lib/core/c++/hlbm.hpp
index 8f267d5..0373125 100644
--- a/lib/core/c++/hlbm.hpp
+++ b/lib/core/c++/hlbm.hpp
@@ -40,15 +40,17 @@ public:
 		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
 		saw::data<sch::Vector<T,Descriptor::D>> vel = vel_f.at(index);
 
-		compute_rho_u<T,Descriptor::D>(dfs_old_f.at(index), rho, vel);
+		compute_rho_u<T,Descriptor>(dfs_old_f.at(index), rho, vel);
 
 		auto& porosity = porosity_f.at(index);
-		auto flip_porosity = saw::data<sch::Scalar<T>>{1.0} - porosity;
+		saw::data<sch::Scalar<T>> one;
+		one.at({}) = 1.0;
+		auto flip_porosity = one - porosity;
 
 		auto& N = particle_N_f.at(index);
 		auto& D = particle_D_f.at(index);
 		// Convex combination of velocities
-		vel = vel * flip_porosity + N * porosity / D;
+		vel = vel * porosity + N * flip_porosity / D;
 		// Equilibrium
 		auto eq = equilibrium<T,Descriptor>(rho,vel);
 
@@ -57,7 +59,7 @@ public:
 		}
 
 		porosity.at({}) = 1.0;
-
+		D.at({}) = 0.0;
 	}
 };