Merge branch 'dev'

4 files changed, 103 insertions, 81 deletions

diff --git a/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp
index fb31f47..98c72a2 100644
--- a/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp
+++ b/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp
@@ -135,6 +135,7 @@ saw::error_or<void> setup_initial_conditions(
 			por_f.at(index).at({}) = {0};
 			rho.at({}) = {1};
 			auto& vel = vel_f.at(index);
+			vel.at({{0u}}) = 0.0002f;
 			auto eq = equilibrium<T,Desc>(rho,vel);
 
 			df = eq;
@@ -149,9 +150,6 @@ saw::error_or<void> setup_initial_conditions(
 			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;
@@ -214,17 +212,30 @@ saw::error_or<void> step(
 		}
 		component<T,Desc,cmpt::OneParticleAt,encode::Sycl<saw::encode::Native>> opa{p_pos,rad,eps};
 
-		component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> flow_in{
-			[&](){
-				uint64_t target_t_i = 16u;
+		component<T,Desc,cmpt::ZouHeVelocityX<true>,encode::Sycl<saw::encode::Native>> flow_in{
+			[&]() -> saw::data<sch::Vector<T,Desc::D>>{
+				saw::data<sch::Vector<T,Desc::D>> vel;
+				uint64_t target_t_i = 8u;
 				if(t_i.get() < target_t_i){
-					return 1.0 + (0.00002 / target_t_i) * t_i.get();
+					vel.at({{0u}}) = (0.0002 / target_t_i) * t_i.get();
+				}else{
+					vel.at({{0u}}) = 0.0002;
 				}
-				return 1.00002;
+				return vel;
+			}()
+		};
+		component<T,Desc,cmpt::ZouHeVelocityX<false>,encode::Sycl<saw::encode::Native>> flow_out{
+			[&]() -> saw::data<sch::Vector<T,Desc::D>>{
+				saw::data<sch::Vector<T,Desc::D>> vel;
+				uint64_t target_t_i = 8u;
+				if(t_i.get() < target_t_i){
+					vel.at({{0u}}) = (0.0002 / target_t_i) * t_i.get();
+				}else{
+					vel.at({{0u}}) = 0.0002;
+				}
+				return vel;
 			}()
 		};
-		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;
@@ -365,7 +376,9 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 	sycl_q.wait();
 	auto lsd_view = make_view(lbm_sycl_data);
 	auto lsdm_view = make_view(lbm_sycl_macro_data);
-	saw::data<sch::UInt64> time_steps{16u*4096ul};
+	
+	saw::data<sch::UInt64> time_steps{32u*4096ul};
+	
 	auto& info_f = lsd_view.template get<"info">();
 
 	for(saw::data<sch::UInt64> i{0u}; i < time_steps and krun; ++i){
diff --git a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp
index 9a04211..e0cbda2 100644
--- a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp
+++ b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp
@@ -132,6 +132,7 @@ saw::error_or<void> setup_initial_conditions(
 			por_f.at(index).at({}) = {1};
 			rho.at({}) = {1};
 			auto& vel = vel_f.at(index);
+			vel.at({{0u}}) = 0.0002f;
 			auto eq = equilibrium<T,Desc>(rho,vel);
 
 			df = eq;
@@ -146,9 +147,6 @@ saw::error_or<void> setup_initial_conditions(
 			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;
@@ -214,17 +212,30 @@ saw::error_or<void> step(
 		
 		component<T,Desc,cmpt::ForceGather, encode::Sycl<saw::encode::Native>> fg;
 
-		component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> flow_in{
-			[&](){
-				uint64_t target_t_i = 16u;
+		component<T,Desc,cmpt::ZouHeVelocityX<true>,encode::Sycl<saw::encode::Native>> flow_in{
+			[&]() -> saw::data<sch::Vector<T,Desc::D>>{
+				saw::data<sch::Vector<T,Desc::D>> vel;
+				uint64_t target_t_i = 8u;
 				if(t_i.get() < target_t_i){
-					return 1.0 + (0.00002 / target_t_i) * t_i.get();
+					vel.at({{0u}}) = (0.0002 / target_t_i) * t_i.get();
+				}else{
+					vel.at({{0u}}) = 0.0002;
 				}
-				return 1.00002;
+				return vel;
+			}()
+		};
+		component<T,Desc,cmpt::ZouHeVelocityX<false>,encode::Sycl<saw::encode::Native>> flow_out{
+			[&]() -> saw::data<sch::Vector<T,Desc::D>>{
+				saw::data<sch::Vector<T,Desc::D>> vel;
+				uint64_t target_t_i = 8u;
+				if(t_i.get() < target_t_i){
+					vel.at({{0u}}) = (0.0002 / target_t_i) * t_i.get();
+				}else{
+					vel.at({{0u}}) = 0.0002;
+				}
+				return vel;
 			}()
 		};
-		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;
@@ -367,7 +378,9 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 	sycl_q.wait();
 	auto lsd_view = make_view(lbm_sycl_data);
 	auto lsdm_view = make_view(lbm_sycl_macro_data);
-	saw::data<sch::UInt64> time_steps{16u*4096ul};
+	
+	saw::data<sch::UInt64> time_steps{32u*4096ul};
+	
 	auto& info_f = lsd_view.template get<"info">();
 
 	for(saw::data<sch::UInt64> i{0u}; i < time_steps and krun; ++i){
diff --git a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp
index 9781691..9d678a8 100644
--- a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp
+++ b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp
@@ -134,6 +134,8 @@ saw::error_or<void> setup_initial_conditions(
 			por_f.at(index).at({}) = {1};
 			rho.at({}) = {1};
 			auto& vel = vel_f.at(index);
+			vel.at({{0u}}) = 0.0002f;
+
 			auto eq = equilibrium<T,Desc>(rho,vel);
 
 			df = eq;
@@ -148,9 +150,6 @@ saw::error_or<void> setup_initial_conditions(
 			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;
@@ -160,48 +159,6 @@ saw::error_or<void> setup_initial_conditions(
 		{{1u,1u}}
 	);
 
-	saw::data<sch::Scalar<T>> eps;
-	eps.at({}).set(1.5f);
-	saw::data<sch::Scalar<T>> rad;
-	rad.at({}).set(dim_y*0.1);
-	saw::data<sch::Vector<T,Desc::D>> p_pos;
-	{
-		p_pos.at({{0u}}) = dim_x * 0.5;
-		p_pos.at({{1u}}) = dim_y * 0.5;
-	}
-
-	/*
-	component<T,Desc,cmpt::OneParticleAt> opa{p_pos,rad,eps};
-	iterator<Desc::D>::apply(
-		[&](auto& index){
-			opa.apply(macros,index,{});
-		},
-		{},// 0-index
-		df_f.get_dims()
-	);
-	*/
-	
-	/*
-	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()
-	);
-	*/
-
 	return saw::make_void();
 }
 
@@ -233,16 +190,30 @@ saw::error_or<void> step(
 		}
 		component<T,Desc,cmpt::OneParticleAt,encode::Sycl<saw::encode::Native>> opa{p_pos,rad,eps};
 
-		component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> flow_in{
-			[&](){
-				uint64_t target_t_i = 16u;
+		component<T,Desc,cmpt::ZouHeVelocityX<true>,encode::Sycl<saw::encode::Native>> flow_in{
+			[&]() -> saw::data<sch::Vector<T,Desc::D>>{
+				saw::data<sch::Vector<T,Desc::D>> vel;
+				uint64_t target_t_i = 8u;
+				if(t_i.get() < target_t_i){
+					vel.at({{0u}}) = (0.0002 / target_t_i) * t_i.get();
+				}else{
+					vel.at({{0u}}) = 0.0002;
+				}
+				return vel;
+			}()
+		};
+		component<T,Desc,cmpt::ZouHeVelocityX<false>,encode::Sycl<saw::encode::Native>> flow_out{
+			[&]() -> saw::data<sch::Vector<T,Desc::D>>{
+				saw::data<sch::Vector<T,Desc::D>> vel;
+				uint64_t target_t_i = 8u;
 				if(t_i.get() < target_t_i){
-					return 1.0 + (0.00002 / target_t_i) * t_i.get();
+					vel.at({{0u}}) = (0.0002 / target_t_i) * t_i.get();
+				}else{
+					vel.at({{0u}}) = 0.0002;
 				}
-				return 1.00002;
+				return vel;
 			}()
 		};
-		component<T,Desc,cmpt::ZouHeHorizontal<false>,encode::Sycl<saw::encode::Native>> flow_out{1.0};
 
 		component<T,Desc,cmpt::ForceGather, encode::Sycl<saw::encode::Native>> fg;
 
@@ -389,7 +360,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 	auto lsd_view = make_view(lbm_sycl_data);
 	auto lsdm_view = make_view(lbm_sycl_macro_data);
 
-	saw::data<sch::UInt64> time_steps{16u*4096ul};
+	saw::data<sch::UInt64> time_steps{32u*4096ul};
 
 	auto& info_f = lsd_view.template get<"info">();
 
diff --git a/modules/core/c++/boundary.hpp b/modules/core/c++/boundary.hpp
index 01ae7b5..22a6f34 100644
--- a/modules/core/c++/boundary.hpp
+++ b/modules/core/c++/boundary.hpp
@@ -221,14 +221,14 @@ public:
 	}
 };
 
-/*
-template<typename FP, typename Descriptor, bool East, typename Encode>
-class component<FP, Descriptor, cmpt::ZouHeVelocityX<East>, Encode> final {
+
+template<typename T, typename Descriptor, bool East, typename Encode>
+class component<T, Descriptor, cmpt::ZouHeVelocityX<East>, Encode> final {
 private:
-	saw::data<sch::Vector<FP,Descriptor::D>> velocity_;
+	saw::data<sch::Vector<T,Descriptor::D>> velocity_;
 public:
 	component(
-		saw::data<sch::Vector<FP,Descriptor::D>> velocity__
+		saw::data<sch::Vector<T,Descriptor::D>> velocity__
 	):
 		velocity_{velocity__}
 	{}
@@ -237,13 +237,38 @@ public:
 	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step) const {
 		
 		bool is_even = ((time_step.get() % 2u) == 0u);
-		using dfi = df_info<FP,Descriptor>;
+		using dfi = df_info<T,Descriptor>;
 		
 		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
 
+
+		auto& dfs = dfs_old_f.at(index);
+
+		saw::data<T> dir_sum;
+		if constexpr (East) {
+			dir_sum = dfs.at({1u}) + dfs.at({5u}) + dfs.at({7u});
+		} else if constexpr(not East){
+			dir_sum = dfs.at({2u}) + dfs.at({6u}) + dfs.at({8u});
+		}
+
+		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)});
+
+		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;
+		}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;
+		}
 	}
 };
-*/
 
 }
 }