End of day, building util funcs

3 files changed, 54 insertions, 2 deletions

diff --git a/c++/util.hpp b/c++/util.hpp
new file mode 100644
index 0000000..a215cc1
--- /dev/null
+++ b/c++/util.hpp
@@ -0,0 +1,13 @@
+#pragma once
+
+namespace kel {
+namespace lbm {
+template<typename T, typename Descriptor>
+struct is_neighbour {
+	template<typename CellFieldSchema>
+	static bool operator()(saw::data<CellFieldSchema>& latt,  saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>& index, saw::data<sch::UInt64> range) {
+		using dfi = df_info<T,Descriptor>;
+	}
+};
+}
+}
diff --git a/examples/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu.cpp
index e1e6333..b24ca38 100644
--- a/examples/cavity_2d_gpu.cpp
+++ b/examples/cavity_2d_gpu.cpp
@@ -197,7 +197,7 @@ void lbm_step(
   /**
 	* 1. Relaxation parameter \tau
 	*/
-	component<sch::T, sch::D2Q9, cmpt::BGK> coll{0.5384};
+	component<sch::T, sch::D2Q9, cmpt::BGK> coll{0.59};
 	component<sch::T, sch::D2Q9, cmpt::BounceBack> bb;
 
 	component<sch::T, sch::D2Q9, cmpt::MovingWall> bb_lid;
diff --git a/examples/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d.cpp
index 286b211..9498ba7 100644
--- a/examples/poiseulle_particles_channel_2d.cpp
+++ b/examples/poiseulle_particles_channel_2d.cpp
@@ -256,6 +256,28 @@ void set_geometry(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 		}
 
 	}, {{0u,0u}}, meta, {{4u,1u}});
+
+	/**
+	 * Set channel circular obstacle
+	 */
+	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+		auto& cell = latt(index);
+		auto& info = cell.template get<"info">();
+
+		saw::data<sch::FixedArray<sch::UInt64,2u>> area{{meta.at({0u})-4u, meta.at({1u})-4u}};
+		double pos_x = 0.8 * area.at({0u}).get();
+		double pos_y = 0.5 * area.at({1u}).get();
+
+		double dist_x = pos_x - index.at({0u}).get();
+		double dist_y = pos_y - index.at({1u}).get();
+
+		double r = 16.0;
+		double dist_sq = dist_x * dist_x + dist_y * dist_y;
+		if(dist_sq < r*r){
+			info({0u}).set(2u);
+		}
+
+	}, {{0u,0u}}, meta, {{2u,2u}});
 }
 
 void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
@@ -290,7 +312,7 @@ void add_particles(kel::lbm::particle_system<kel::lbm::sch::T,2u>& part_sys){
 	auto& p_mask = part.template get<"mask">();
 	{
 		particle_circle_geometry<sch::T> geo;
-		p_mask = geo.template generate_mask<sch::T>(16u);
+		p_mask = geo.template generate_mask<sch::T>(2u,1u);
 	}
 	auto& rigid_body = part.template get<"rigid_body">();
 	auto& p_size = part.template get<"size">();
@@ -417,6 +439,23 @@ void couple_particles_to_lattice(
 
 				// Add forces to put away from walls
 				/// 1. Check if neighbour is wall
+				for(saw::data<sch::UInt64> k{0u}; k.get() < sch::D2Q9::Q; ++k){
+					auto n_p_cell_pos = saw::data<sch::FixedArray<sch::UInt64,2u>> {{
+						p_cell_pos.at({{0u}}) + dfi::directions[k.get()][0u],
+						p_cell_pos.at({{1u}}) + dfi::directions[k.get()][1u]
+					}};
+
+					auto& n_cell = latt(n_p_cell_pos);
+					auto& n_info = n_cell.template get<"info">()({0u});
+
+					// If neighbour is wall, then add force pushing the particle away
+					if(n_info.get() <= 1u){
+						// add to p_acc
+						p_acc.at({{0u}}) = p_acc.at({{0u}}) + saw::data<sch::Int32>{10 * dfi::directions[dfi::opposite_index[k.get()]][0u]}.template cast_to<sch::T>();
+						p_acc.at({{1u}}) = p_acc.at({{1u}}) + saw::data<sch::Int32>{10 * dfi::directions[dfi::opposite_index[k.get()]][1u]}.template cast_to<sch::T>();
+						;
+					}
+				}
 				/// 2. Add force pushing away from wall
 
 				// Add forces to push away from other particles