Merge branch 'dev'HEADmaster

8 files changed, 174 insertions, 49 deletions

diff --git a/default.nix b/default.nix
index dc75b57..02f9213 100644
--- a/default.nix
+++ b/default.nix
@@ -151,6 +151,17 @@ in rec {
 			inherit pname version stdenv forstio adaptive-cpp;
 			inherit kel;
 		};
+
+		poiseulle_particles_2d_all_gpu = pkgs.symlinkJoin {
+			name = "poiseulle_particles_2d_all_gpu";
+			paths = [
+				examples.poiseulle_particles_2d_bgk_gpu
+				examples.poiseulle_particles_2d_psm_gpu
+				examples.poiseulle_particles_2d_hlbm_gpu
+				examples.poiseulle_particles_2d_fplbm_gpu
+				examples.poiseulle_particles_2d_ibm_gpu
+			];
+		};
 		
 		poiseulle_moving_particle_2d_psm_gpu = pkgs.callPackage ./examples/poiseulle_moving_particle_2d_psm_gpu/.nix/derivation.nix {
 			inherit pname version stdenv forstio adaptive-cpp;
diff --git a/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp b/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp
index 0c10d38..cae8a4d 100644
--- a/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp
+++ b/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp
@@ -1,6 +1,7 @@
 #include <kel/lbm/sycl/lbm.hpp>
 #include <kel/lbm/lbm.hpp>
 #include <kel/lbm/particle.hpp>
+#include <kel/lbm/particle_opa.hpp>
 
 #include <forstio/io/io.hpp>
 #include <forstio/remote/filesystem/easy.hpp>
@@ -157,25 +158,6 @@ saw::error_or<void> setup_initial_conditions(
 		df_f.get_dims(),
 		{{1u,1u}}
 	);
-	
-	iterator<Desc::D>::apply(
-		[&](auto& index){
-			saw::data<sch::Vector<T,Desc::D>> middle, ind_vec;
-			middle.at({{0u}}) = dim_x * 0.5;
-			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();
 }
@@ -190,6 +172,8 @@ saw::error_or<void> step(
 	auto& q = dev.get_handle();
 	auto& info_f = fields.template get<"info">();
 	auto& porous_f = macros.template get<"porosity">();
+	sch::data<sch::Scalar<T>> eps;
+	eps.at({}).set(1.5f);
 
 	q.submit([&](acpp::sycl::handler& h){
 		component<T,Desc,cmpt::PSM,encode::Sycl<saw::encode::Native>> collision{0.8};
@@ -199,7 +183,7 @@ saw::error_or<void> step(
 		component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> flow_in{1.0};
 		component<T,Desc,cmpt::ZouHeHorizontal<false>,encode::Sycl<saw::encode::Native>> flow_out{1.0};
 
-		component<T,Desc,cmpt::OneParticleAt, encode::Sycl<saw::encode::Native>> opa{{},{}};
+		component<T,Desc,cmpt::OneParticleAt, encode::Sycl<saw::encode::Native>> opa{{},{},eps};
 
 		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;
@@ -371,6 +355,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 			}
 		}
 		*/
+		/*
 		if(i.get() % 32u == 0u){
 			{
 				auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr);
@@ -385,6 +370,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 				}
 			}
 		}
+		*/
 		// Stream
 		sycl_q.submit([&](acpp::sycl::handler& h){
 			component<T,Desc,cmpt::Stream,encode::Sycl<saw::encode::Native>> stream;
@@ -424,6 +410,18 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 			return eov;
 		}
 	}
+	{
+		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",time_steps.get(), *lbm_macro_data_ptr);
+		if(eov.is_error()){
+			return eov;
+		}
+	}
 
 	sycl_q.wait();
 	return saw::make_void();
diff --git a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp
index 9cc77ae..73964b7 100644
--- a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp
@@ -366,6 +366,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 			}
 		}
 		*/
+		/*
 		if(i.get() % 32u == 0u){
 			{
 				auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr);
@@ -380,6 +381,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 				}
 			}
 		}
+		*/
 		// Stream
 		sycl_q.submit([&](acpp::sycl::handler& h){
 			component<T,Desc,cmpt::Stream,encode::Sycl<saw::encode::Native>> stream;
@@ -419,6 +421,18 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 			return eov;
 		}
 	}
+	{
+		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",time_steps.get(), *lbm_macro_data_ptr);
+		if(eov.is_error()){
+			return eov;
+		}
+	}
 
 	sycl_q.wait();
 	return saw::make_void();
diff --git a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
index 0337158..b3ea7ae 100644
--- a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
@@ -1,6 +1,7 @@
 #include <kel/lbm/sycl/lbm.hpp>
 #include <kel/lbm/lbm.hpp>
 #include <kel/lbm/particle.hpp>
+#include <kel/lbm/particle/particle_opa.hpp>
 
 #include <forstio/io/io.hpp>
 #include <forstio/remote/filesystem/easy.hpp>
@@ -157,7 +158,29 @@ saw::error_or<void> setup_initial_conditions(
 		df_f.get_dims(),
 		{{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.25;
+		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;
@@ -176,6 +199,7 @@ saw::error_or<void> setup_initial_conditions(
 		{},// 0-index
 		df_f.get_dims()
 	);
+	*/
 
 	return saw::make_void();
 }
@@ -203,6 +227,18 @@ saw::error_or<void> step(
 
 		component<T,Desc,cmpt::Equilibrium,encode::Sycl<saw::encode::Native>> equi{rho_b,vel_b};
 
+
+		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.25;
+			p_pos.at({{1u}}) = dim_y * 0.5;
+		}
+		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 = 64u;
@@ -230,6 +266,7 @@ saw::error_or<void> step(
 					bb.apply(fields,index,t_i);
 				break;
 				case 2u:
+					opa.apply(macros,index,t_i);
 					collision.apply(fields,macros,index,t_i);
 					break;
 				case 3u:
@@ -385,6 +422,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 			}
 		}
 		*/
+		/*
 		if(i.get() % 32u == 0u){
 			{
 				auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr);
@@ -399,6 +437,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
 				}
 			}
 		}
+		*/
 		// Stream
 		sycl_q.submit([&](acpp::sycl::handler& h){
 			component<T,Desc,cmpt::Stream,encode::Sycl<saw::encode::Native>> stream;
diff --git a/lib/core/c++/particle/particle_opa.hpp b/lib/core/c++/particle/particle_opa.hpp
index 4588a55..a5e6c28 100644
--- a/lib/core/c++/particle/particle_opa.hpp
+++ b/lib/core/c++/particle/particle_opa.hpp
@@ -1,7 +1,9 @@
 #pragma once
 
-#include "common.hpp"
 #include "../component.hpp"
+#include "common.hpp"
+#include "schema.hpp"
+#include "porosity.hpp"
 
 namespace kel {
 namespace lbm {
@@ -35,11 +37,12 @@ public:
 
 		auto& porous = porous_f.at(index);
 
-
 		auto pos_ind = saw::math::vectorize_data(index);
 
-		auto diff = pos_ind - pos_;
+		// 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_);
 	}
 };
 }
diff --git a/lib/core/c++/particle/porosity.hpp b/lib/core/c++/particle/porosity.hpp
index f555cae..11185c5 100644
--- a/lib/core/c++/particle/porosity.hpp
+++ b/lib/core/c++/particle/porosity.hpp
@@ -1,16 +1,21 @@
 #pragma once
 
-#include "particle.hpp"
+#include "common.hpp"
+#include "schema.hpp"
 #include "../math/n_closest.hpp"
 
 namespace kel {
 namespace lbm {
+namespace por {
+template<typename T>
+struct ParticleSpheroid {};
+}
+
 template<typename T, uint64_t D, typename Coll>
 class particle_porosity {
 public:
 
-	
-	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,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">();
@@ -27,20 +32,47 @@ public:
 	}
 };
 
-
-
-template<typename T, uint64_t D, typename saw::native_data_type<T>::type radius, typename saw::native_data_type<T>::type eps>
-class particle_porosity<T, D, coll::ParticleCollisionSpheroid<T,radius, eps>> final {
+template<typename T, uint64_t D>
+class particle_porosity<T, D, por::ParticleSpheroid<T>> final {
 public:
-	saw::data<sch::Scalar<T>> calculate(const saw::data<sch::Vector<T,D>>& lbm_pos, saw::data<sch::Scalar<T>> rad) const {
-		saw::data<sch::Scalar<T>> pos;
-		
+	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>> eps_h;
+		eps_h.at({}) = eps.at({}).get() / 2;
+
+		auto rad_low = rad - eps_h;
+		auto rad_high = rad + eps_h;
+
+		if(s_dist.at({}) <= rad_low.at({})){
+			por.at({}).set(1);
+			return por;
+		}
+		if(s_dist.at({}) >= rad_high.at({})){
+			por.at({}).set(0);
+			return por;
+		}
+
+		{
+			typename saw::native_data_type<T>::type inner = (std::numbers::pi / ( eps.at({}).get() )) * (s_dist - rad_low).at({}).get();
+			auto cos_inner = std::cos(inner);
+			auto cos_i_2 = cos_inner * cos_inner;
+			por.at({}).set(cos_i_2);
+		}
+
+		return por;
 	}
+};
 
-	saw::data<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T,radius,eps> > >& part_group, uint64_t i, const saw::data<sch::Vector<T,D>>& lbm_pos) const {
+
+template<typename T, uint64_t D>
+class particle_porosity<T, D, 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) {
 		saw::data<sch::Scalar<T>> por;
-		por.at({});
 
 		auto& parts = part_group.template get<"particles">();
 		auto& pi = parts.at({i});
@@ -51,18 +83,19 @@ public:
 		// Basically the queried position minus the center of the particle
 		auto dist = lbm_pos - pi_rb_pos;
 
-		saw::data<sch::Scalar<T>> dist_len;
-		for(uint64_t i{0u}; i < D; ++i){
-			auto& d_i = dist.at({{i}});
-			dist_len.at({}) = d_i * d_i;
-		}
-		dist_len.at({}).set(std::sqrt(dist_len.at({}).get());
+		saw::data<sch::Scalar<T>> dist_len = saw::math::dot(dist,dist);
+		dist_len.at({}).set(std::sqrt(dist_len.at({}).get()));
 
-		saw::data<sch::Scalar<T>> rad_d{radius};
-		saw::data<sch::Scalar<T>> eps_half{eps *0.5};
+		auto& coll = part_group.template get<"collision">();
+		const saw::data<sch::Scalar<T>>& rad_d = coll.template get<"radius">();
+		const auto& eps = part_group.template get<"epsilon">();
+		// Move this somewhere
 		
-		saw::data<sch::Scalar<T>> rad_d_eps_p = rad_d + eps_half;
-		saw::data<sch::Scalar<T>> rad_d_eps_n = rad_d - eps_half;
+		saw::data<sch::Scalar<T>> eps_h;
+		eps_h.at({}) = eps.at({}).get() / 2;
+
+		saw::data<sch::Scalar<T>> rad_d_eps_p = rad_d + eps_h;
+		saw::data<sch::Scalar<T>> rad_d_eps_n = rad_d - eps_h;
 
 		// saw::data<sch::Scalar<T>> rad_d_eps_p_2 = rad_d_eps_p * rad_d_eps_p;
 		// saw::data<sch::Scalar<T>> rad_d_eps_n_2 = rad_d_eps_n * rad_d_eps_n;
@@ -71,19 +104,19 @@ public:
 
 		saw::data<sch::Scalar<T>> inside;
 		if(dist_len.at({}) <= rad_d_eps_n.at({})){
-			inside.at({}).set(1);
+			inside.at({}).set(0);
 			return inside;
 		}
 		if(dist_len.at({}) > rad_d_eps_p.at({})){
-			inside.at({}).set(0);
-			return inside
+			inside.at({}).set(1);
+			return inside;
 		}
 
 		/*
 		 * cos^2 ( ||x-X(t)||_2 - (R-eps/2) )
 		 */
 		{
-			typename saw::native_data_type<T>::type inner = (std::pi / ( 2.0 * eps )) * (dist_len - rad_d_eps_n).at({}).get();
+			typename saw::native_data_type<T>::type inner = (std::numbers::pi / (eps)) * (dist_len - rad_d_eps_n).at({}).get();
 			auto cos_inner = std::cos(inner);
 			auto cos_i_2 = cos_inner * cos_inner;
 			inside.at({}).set(cos_i_2);
diff --git a/lib/core/c++/particle/schema.hpp b/lib/core/c++/particle/schema.hpp
index 18a697a..714a16f 100644
--- a/lib/core/c++/particle/schema.hpp
+++ b/lib/core/c++/particle/schema.hpp
@@ -56,6 +56,7 @@ template<typename T, uint64_t D, typename CollisionType = coll::Spheroid<T>>
 using ParticleGroup = Struct<
 	Member<Array<T,D>, "mask">,
 	Member<FixedArray<typename CollisionType::Schema,1u>, "collision">,
+	Member<FixedArray<Scalar<T>,1u>, "epsilon">,
 	Member<FixedArray<Scalar<T>,1u>, "mask_step">,
 	Member<FixedArray<Scalar<T>,1u>, "density">,
 	Member<FixedArray<Vector<T,D>,1u>, "center_of_mass">,
diff --git a/scripts/python/csv_l2_norm.py b/scripts/python/csv_l2_norm.py
new file mode 100755
index 0000000..88f4817
--- /dev/null
+++ b/scripts/python/csv_l2_norm.py
@@ -0,0 +1,26 @@
+#!/usr/bin/env python3
+
+import argparse
+import numpy as np
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Compute the L2 Norm");
+    parser.add_argument("one");
+    parser.add_argument("two");
+    args = parser.parse_args();
+
+    a = np.loadtxt(args.one, delimiter=",");
+    b = np.loadtxt(args.two, delimiter=",");
+
+    if a.shape != b.shape:
+        raise ValueError(f"Shape mismatch: {a.shape} vs {b.shape}");
+    #endif
+
+    print(np.linalg.norm(a-b));
+
+#enddef
+
+if __name__ == "__main__":
+    main();
+#enddef