Adding 3D poiseulle and using http fetching

4 files changed, 161 insertions, 14 deletions

diff --git a/examples/poiseulle_3d/.nix/derivation.nix b/examples/poiseulle_3d/.nix/derivation.nix
index a06e14b..f00b6df 100644
--- a/examples/poiseulle_3d/.nix/derivation.nix
+++ b/examples/poiseulle_3d/.nix/derivation.nix
@@ -6,6 +6,7 @@
 , pname
 , version
 , kel-lbm
+, adaptive-cpp
 }:
 
 stdenv.mkDerivation {
@@ -26,6 +27,7 @@ stdenv.mkDerivation {
 		forstio.codec-json
 		forstio.remote
 		forstio.remote-sycl
+		adaptive-cpp
 		kel-lbm
 	];
 
diff --git a/examples/poiseulle_3d/SConscript b/examples/poiseulle_3d/SConscript
index 02338cd..078541a 100644
--- a/examples/poiseulle_3d/SConscript
+++ b/examples/poiseulle_3d/SConscript
@@ -11,6 +11,8 @@ dir_path = Dir('.').abspath
 
 # Environment for base library
 examples_env = env.Clone();
+examples_env['CXX'] = 'syclcc-clang';
+examples_env['CXXFLAGS'] += ['-O3'];
 
 examples_env.sources = sorted(glob.glob(dir_path + "/*.cpp"))
 examples_env.headers = sorted(glob.glob(dir_path + "/*.hpp"))
diff --git a/examples/poiseulle_3d/SConstruct b/examples/poiseulle_3d/SConstruct
index a7201cb..fe206e1 100644
--- a/examples/poiseulle_3d/SConstruct
+++ b/examples/poiseulle_3d/SConstruct
@@ -54,10 +54,13 @@ env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[],
         '-std=c++20',
         '-g',
         '-Wall',
-        '-Wextra'
+        '-Wextra',
+        '-isystem', 'AdaptiveCpp'
     ],
     LIBS=[
-        'forstio-core'
+        'forstio-core',
+        'acpp-rt',
+        'omp'
     ]
 );
 env.__class__.add_source_files = add_kel_source_files
diff --git a/examples/poiseulle_3d/poiseulle_3d.cpp b/examples/poiseulle_3d/poiseulle_3d.cpp
index 68e1bb7..9220a1f 100644
--- a/examples/poiseulle_3d/poiseulle_3d.cpp
+++ b/examples/poiseulle_3d/poiseulle_3d.cpp
@@ -2,6 +2,56 @@
 
 #include <forstio/remote/sycl/sycl.hpp>
 
+namespace saw {
+template<typename Desc, typename CellT, typename Encode>
+class data<kel::lbm::sch::CellField<Desc, CellT>, encode::Sycl<Encode>> final {
+public:
+	using Schema = kel::lbm::sch::CellField<Desc,CellT>;
+	using MetaSchema = typename meta_schema<Schema>::MetaSchema;
+private:
+	data<schema::Array<CellT,Desc::D>, encode::Sycl<Encode>> inner_;
+public:
+	data() = delete;
+	data(const data<MetaSchema,Encode>& inner_meta__, acpp::sycl::queue& q):
+		inner_{inner_meta__,q}
+	{}
+
+	const data<MetaSchema, Encode> meta() const {
+		return inner_.dims();
+	}
+
+	template<uint64_t i>
+	data<schema::UInt64,Encode> get_dim_size() const {
+		static_assert(i < Desc::D, "Not enough dimensions");
+		return inner_.template get_dim_size<i>();
+	}
+
+	const data<CellT>& operator()(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index)const{
+		return inner_.at(index);
+	}
+
+	data<CellT>& operator()(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index){
+		return inner_.at(index);
+	}
+
+	const data<CellT>& at(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index)const{
+		return inner_.at(index);
+	}
+
+	data<CellT>& at(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index){
+		return inner_.at(index);
+	}
+
+	data<schema::UInt64,Encode> internal_size() const {
+		return inner_.internal_size();
+	}
+
+	data<CellT,Encode>* internal_data() {
+		return inner_.internal_data();
+	}
+};
+}
+
 namespace kel {
 namespace lbm {
 namespace sch {
@@ -21,8 +71,8 @@ using CellStruct = Struct<
 >;
 }
 
-template<typename StructFieldSchema>
-saw::error_or<void> set_geometry(saw::data<StructFieldSchema>& latt){
+template<typename StructFieldSchema, typename Encode>
+saw::error_or<void> set_geometry(saw::data<StructFieldSchema, Encode>& latt){
 	using namespace kel::lbm;
 
 	auto meta = latt.meta();
@@ -109,11 +159,96 @@ saw::error_or<void> set_geometry(saw::data<StructFieldSchema>& latt){
 	return saw::make_void();
 }
 
-saw::error_or<void> set_initial_conditions(saw::data<StructFieldSchema>& latt){
+template<typename StructFieldSchema, typename Encode>
+saw::error_or<void> set_initial_conditions(saw::data<StructFieldSchema, Encode>& latt){
+
+	saw::data<sch::T> rho{1.0};
+	saw::data<sch::FixedArray<sch::T,sch::Desc::D>> vel{};
+	vel.at({0u}) = {0.1f};
+	auto eq = equilibrium<sch::T,sch::Desc>(rho, vel);
+
+	auto meta = latt.meta();
+
+	/**
+	 * Set distribution
+	 */
+	iterator<sch::Desc::D>::apply([&](const saw::data<sch::FixedArray<sch::UInt64,sch::Desc::D>>& index){
+		auto& cell = latt.at(index);
+		auto& dfs = cell.template get<"dfs">();
+		auto& dfs_old = cell.template get<"dfs_old">();
+
+		for(saw::data<sch::UInt64> k = 0; k < saw::data<sch::UInt64>{sch::Desc::Q}; ++k){
+			dfs(k) = eq.at(k);
+			dfs_old(k) = eq.at(k);
+		}
+	}, {}, meta);
 	return saw::make_void();
 }
 
-saw::error_or<void> lbm_step(saw::data<StructFieldSchema>& latt){
+template<typename StructFieldSchema, typename Encode>
+saw::error_or<void> lbm_step(
+	saw::data<StructFieldSchema, Encode>& latt,
+	saw::data<sch::UInt64> time_step,
+	acpp::sycl::queue& sycl_q
+){
+	using dfi = df_info<sch::T,sch::Desc>;
+	auto meta = latt.meta();
+	bool even_step = ((time_step.get() % 2u) == 0u);
+
+	component<sch::T, sch::Desc, cmpt::BGK, saw::encode::Sycl<saw::encode::Native>> coll{0.5384};
+	component<sch::T, sch::Desc, cmpt::BounceBack, saw::encode::Sycl<saw::encode::Native>> bb;
+
+	/**
+	 * Collision
+	 */
+	iterator<sch::Desc::D>::apply([&](const saw::data<sch::FixedArray<sch::UInt64,sch::Desc::D>>& index){
+		auto& cell = latt.at(index);
+		auto& info = cell.template get<"info">();
+
+		switch(info({0u}).get()){
+			case 1u: {
+				bb.apply(latt, index, time_step);
+				break;
+			}
+			case 2u: {
+				coll.apply(latt, index, time_step);
+				break;
+			}
+			case 3u: {
+				//inlet.apply(latt, index, time_step);
+				break;
+			}
+			case 4u: {
+				//outlet.apply(latt, index, time_step);
+				break;
+			}
+			default:
+				break;
+		}
+
+	}, {{0u,0u}}, meta);
+
+	// Stream
+	iterator<sch::Desc::D>::apply([&](const saw::data<sch::FixedArray<sch::UInt64,sch::Desc::D>>& index){
+		auto& cell = latt(index);
+		auto& df_new = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">();
+		auto& info_new = cell.template get<"info">();
+
+		if(info_new({0u}).get() > 0u){
+			for(uint64_t k = 0u; k < sch::Desc::Q; ++k){
+				auto dir = dfi::directions[dfi::opposite_index[k]];
+				auto index_dir = index;
+				for(uint64_t i = 0; i < sch::Desc::D; ++i){
+					index_dir.at({i}) = index_dir.at({i}) + dir[i];
+				}
+				auto& cell_dir_old = latt.at(index_dir);
+
+				auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">();
+
+				df_new({k}) = df_old({k});
+			}
+		}
+	}, {{0u,0u}}, meta);
 	return saw::make_void();
 }
 
@@ -126,31 +261,36 @@ saw::error_or<void> real_main(int argc, char** argv){
 	};
 
 	print_lbm_meta<sch::T,sch::Desc>(conv, {1e-3f});
+	acpp::sycl::queue sycl_q;
 
 	/**
 	 * Set the meta and initialize the lattice
 	 */
 	saw::data<sch::FixedArray<sch::UInt64,sch::Desc::D>> meta{{2048u,128u,128u}};
-	saw::data<sch::CellField<sch::Desc,sch::CellStruct>, encode::Sycl<encode::Native>> lattice{meta};
+	saw::data<sch::CellField<sch::Desc,sch::CellStruct>, saw::encode::Sycl<saw::encode::Native>> lattice{meta,sycl_q};
 
 	/**
 	 * Set the geometry
 	 */
-	auto eov = set_geometry(lattice);
-	if(eov.is_error()){
-		return eov;
+	{
+		auto eov = set_geometry(lattice);
+		if(eov.is_error()){
+			return eov;
+		}
 	}
 
 	/**
 	 * Set initial conditions in the simulations
 	 */
+	{
 	auto eov = set_initial_conditions(lattice);
-	if(eov.is_error()){
-		return eov;
+		if(eov.is_error()){
+			return eov;
+		}
 	}
 
-	for(uint64_t i = 0u; i < 1024u*32; ++i){
-		auto eov = lbm_step(lattice);
+	for(saw::data<sch::UInt64> ts{0u}; ts < saw::data<sch::UInt64>{1024u*32u}; ++ts){
+		auto eov = lbm_step(lattice,ts,sycl_q);
 		if(eov.is_error()){
 			return eov;
 		}