Adding 3D poiseulle and using http fetching

10 files changed, 196 insertions, 34 deletions

diff --git a/default.nix b/default.nix
index f2fb618..5f7dc1b 100644
--- a/default.nix
+++ b/default.nix
@@ -4,13 +4,6 @@
 }:
 
 let
-  forstio = (import ((builtins.fetchGit {
-    url = "git@git.keldu.de:forstio/forstio";
-    ref = "dev";
-  }).outPath + "/default.nix"){
-    inherit stdenv;
-    inherit clang-tools;
-	}).forstio;
 
 	adaptive-cpp = pkgs.callPackage .nix/adaptive-cpp.nix {
 		inherit stdenv;
@@ -18,6 +11,15 @@ let
 		lld = pkgs.lld_19;
 	};
 
+	forstio = (import ((builtins.fetchTarball {
+		url = "https://git.keldu.de/forstio-forstio/snapshot/dev.tar.gz";
+		sha256 = "sha256:1mbak3ksppkij0325p6gz8mibbss0gfj4w85fs4f1jkz8wpdbkgq";
+	}) + "/default.nix"){
+		inherit stdenv;
+		inherit clang-tools;
+		inherit adaptive-cpp;
+	}).forstio;
+
 	pname = "kel-lbm";
 	version = "0.0.2";
 in rec {
@@ -52,6 +54,16 @@ in rec {
 			inherit pname version stdenv forstio adaptive-cpp;
 			kel-lbm = lbm;
 		};
+
+		poiseulle_2d = pkgs.callPackage ./examples/poiseulle_3d/.nix/derivation.nix {
+			inherit pname version stdenv forstio;
+			kel-lbm = lbm;
+		};
+
+		poiseulle_3d = pkgs.callPackage ./examples/poiseulle_3d/.nix/derivation.nix {
+			inherit pname version stdenv forstio adaptive-cpp;
+			kel-lbm = lbm;
+		};
 	};
 
 	debug = {
diff --git a/examples/cavity_2d_gpu/.nix/derivation.nix b/examples/cavity_2d_gpu/.nix/derivation.nix
index 01c2a10..0b604f5 100644
--- a/examples/cavity_2d_gpu/.nix/derivation.nix
+++ b/examples/cavity_2d_gpu/.nix/derivation.nix
@@ -7,7 +7,7 @@
 , python3
 , pname
 , version
-, adaptive-cpp-dev
+, adaptive-cpp
 , kel-lbm
 }:
 
@@ -29,7 +29,7 @@ stdenv.mkDerivation {
 		forstio.codec-unit
 		forstio.remote
 		forstio.codec-json
-		adaptive-cpp-dev
+		adaptive-cpp
 		kel-lbm
 	];
 
diff --git a/examples/poiseulle_2d/poiseulle_2d.cpp b/examples/poiseulle_2d/poiseulle_2d.cpp
index 8c31cb8..3bc7a6c 100644
--- a/examples/poiseulle_2d/poiseulle_2d.cpp
+++ b/examples/poiseulle_2d/poiseulle_2d.cpp
@@ -251,7 +251,7 @@ void lbm_step(
 	/**
 	 * Collision
 	 */
-	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+	iterator<sch::D2Q9::D>::apply([&](const saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>>& index){
 		auto& cell = latt(index);
 		auto& info = cell.template get<"info">();
 
@@ -279,7 +279,7 @@ void lbm_step(
 	}, {{0u,0u}}, meta);
 
 	// Stream
-	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
+	iterator<sch::D2Q9::D>::apply([&](const saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::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">();
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;
 		}
diff --git a/lib/c++/boundary.hpp b/lib/c++/boundary.hpp
index 5cc7705..c10784c 100644
--- a/lib/c++/boundary.hpp
+++ b/lib/c++/boundary.hpp
@@ -44,8 +44,8 @@ public:
 	 * Raw setup
 	 */
 	template<typename CellFieldSchema>
-	void apply(saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, uint64_t time_step){
-		bool is_even = ((time_step % 2u) == 0u);
+	void apply(saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step){
+		bool is_even = ((time_step.get() % 2u) == 0u);
 
 		// This is a ref
 		auto& cell = field(index);
@@ -88,10 +88,10 @@ public:
 	{}
 
 	template<typename CellFieldSchema>
-	void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, uint64_t time_step){
+	void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step){
 		using dfi = df_info<FP,Descriptor>;
 
-		bool is_even = ((time_step % 2) == 0);
+		bool is_even = ((time_step.get() % 2) == 0);
 		auto& cell = field(index);
 
 		auto& info = cell.template get<"info">();
diff --git a/lib/c++/collision.hpp b/lib/c++/collision.hpp
index 9ab542b..73e86ae 100644
--- a/lib/c++/collision.hpp
+++ b/lib/c++/collision.hpp
@@ -45,9 +45,9 @@ public:
 	 * Raw setup
 	 */
 	template<typename CellFieldSchema>
-	void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, uint64_t time_step){
-		bool is_even = ((time_step % 2) == 0);
-		auto& cell = field(index);
+	void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step){
+		bool is_even = ((time_step.get() % 2) == 0);
+		auto& cell = field.at(index);
 
 		auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
 		// auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
@@ -90,8 +90,8 @@ public:
 	static constexpr saw::string_literal before = "streaming";
 
 	template<typename CellFieldSchema>
-	void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, uint64_t time_step){
-		bool is_even = ((time_step % 2) == 0);
+	void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, saw::data<sch::UInt64> time_step){
+		bool is_even = ((time_step.get() % 2) == 0);
 		auto& cell = field(index);
 
 		auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
diff --git a/lib/c++/descriptor.hpp b/lib/c++/descriptor.hpp
index 1e3402f..c6938e3 100644
--- a/lib/c++/descriptor.hpp
+++ b/lib/c++/descriptor.hpp
@@ -311,6 +311,9 @@ public:
 	}
 };
 
+/*
+ * Create a cellfield
+ */
 template<typename Desc, typename CellT, typename Encode>
 class data<kel::lbm::sch::CellField<Desc, CellT>, Encode> final {
 public: