Dangling changes

3 files changed, 107 insertions, 341 deletions

diff --git a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
index 964bfde..f69b908 100644
--- a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
+++ b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
@@ -9,49 +9,6 @@
 #include <chrono>
 #include <cmath>
 
-namespace saw {
-namespace encode {
-template<typename InnerEnc>
-struct Sycl {
-};
-}
-
-template<typename Sch, uint64_t Dim>
-class data<schema::Array<Sch, Dim>, encode::Sycl<encode::Native>> {
-public:
-	using Schema = schema::Array<Sch,Dim>;
-private:
-	using SyclKelAllocator = acpp::sycl::usm_allocator<data<Sch, encode::Native>, acpp::sycl::usm::alloc::shared>;
-	std::vector<data<Sch, encode::Native>, SyclKelAllocator> data_;
-	data<schema::UInt64, encode::Native> size_;
-public:
-	data(const data<Schema, encode::Native>& host_data__):
-		data_{&host_data__.at({0u}),host_data__.size().get()},
-		size_{host_data__.size()}
-	{}
-
-	auto& get_handle() {
-		return data_;
-	}
-
-	const auto& get_handle() const {
-		return data_;
-	}
-
-	data<schema::UInt64, encode::Native> internal_size() const {
-		return size_;
-	}
-
-	data<Sch, encode::Native>* internal_data() {
-		return &data_[0u];
-	}
-
-	const data<Sch, encode::Native>* internal_data() const {
-		return data_.data();
-	}
-};
-}
-
 namespace kel {
 namespace lbm {
 namespace sch {
@@ -417,6 +374,8 @@ int main(){
 	}
 	auto stop = std::chrono::steady_clock::now();
 	std::cout<<std::format("{:%H:%M:%S}",(stop-start))<<std::endl;
+
+	sycl_q.wait();
 	sycl::free(info,sycl_q);
 	sycl::free(dfs,sycl_q);
 	sycl::free(dfs_old,sycl_q);
diff --git a/examples/poiseulle_2d_gpu/.nix/derivation.nix b/examples/poiseulle_2d_gpu/.nix/derivation.nix
new file mode 100644
index 0000000..7f5c2b0
--- /dev/null
+++ b/examples/poiseulle_2d_gpu/.nix/derivation.nix
@@ -0,0 +1,39 @@
+{ lib
+, stdenv
+, scons
+, clang-tools
+, forstio
+, python3
+, pname
+, version
+, adaptive-cpp
+, kel-lbm
+}:
+
+stdenv.mkDerivation {
+	pname = pname + "-examples-" + "poiseulle_2d_gpu";
+	inherit version;
+	src = ./..;
+
+	nativeBuildInputs = [
+		scons
+		clang-tools
+		python3
+	];
+
+	buildInputs = [
+		forstio.core
+		forstio.async
+		forstio.codec
+		forstio.codec-unit
+		forstio.remote
+		forstio.codec-json
+		adaptive-cpp
+		kel-lbm.core
+#		kel-lbm.sycl
+	];
+
+	preferLocalBuild = true;
+
+	outputs = [ "out" "dev" ];
+}
diff --git a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp
index a319122..d52af6a 100644
--- a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp
+++ b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp
@@ -1,10 +1,14 @@
-#include "../c++/lbm.hpp"
-#include "../c++/collision.hpp"
-#include "../c++/boundary.hpp"
-#include "../c++/iterator.hpp"
+#include <kel/lbm/lbm.hpp>
+#include <AdaptiveCpp/sycl/sycl.hpp>
 
 #include <forstio/codec/data.hpp>
 
+template<typename T>
+using SyclHostAlloc = acpp::sycl::usm_allocator<saw::data<T>, acpp::sycl::usm::alloc::host>;
+
+template<typename T>
+using SyclDeviceAlloc = acpp::sycl::usm_allocator<saw::data<T>, acpp::sycl::usm::alloc::device>;
+
 namespace kel {
 namespace lbm {
 namespace sch {
@@ -19,24 +23,19 @@ using namespace saw::schema;
  * D factor
  * Q factor
  */
-using T = Float32;
-using D2Q5 = Descriptor<2u,5u>;
+using T = Float64;
 using D2Q9 = Descriptor<2u,9u>;
 
-template<typename Desc>
-using DfCell = Cell<T, Desc, 0u, 0u, 1u>;
-
-template<typename Desc>
+using DfCell = Cell<T, D2Q9, 0u, 0u, 1u>;
 using CellInfo = Cell<UInt8, D2Q9, 1u, 0u, 0u>;
 
 /**
  * Basic type for simulation
  */
-template<typename Desc>
 using CellStruct = Struct<
-	Member<DfCell<Desc>, "dfs">,
-	Member<DfCell<Desc>, "dfs_old">,
-	Member<CellInfo<Desc>, "info">
+	Member<DfCell, "dfs">,
+	Member<DfCell, "dfs_old">,
+	Member<CellInfo, "info">
 >;
 
 template<typename T, uint64_t D>
@@ -45,12 +44,6 @@ using MacroStruct = Struct<
 	Member<T, "pressure">
 >;
 
-template<typename T>
-using GeometryStruct = Struct<
-	Member<T, "info">
->;
-
-using CavityFieldD2Q9 = CellField<D2Q9, CellStruct<D2Q9>>;
 }
 
 namespace cmpt {
@@ -70,169 +63,40 @@ struct PressureBoundaryRestrictedVelocityTo {};
  * 0 - 2 - 2
  *
  */
-template<typename FP,typename Descriptor, bool East>
-struct component<FP,Descriptor, cmpt::PressureBoundaryRestrictedVelocityTo<East>> {
-private:
-	saw::data<FP> pressure_setting_;
-	saw::data<FP> rho_setting_;
-public:
-	component(const saw::data<FP>& pressure_setting__):
-		pressure_setting_{pressure_setting__},
-		rho_setting_{pressure_setting__ * df_info<FP,Descriptor>::inv_cs2}
-	{}
-
-	template<typename CellFieldSchema>
-	void apply(saw::data<CellFieldSchema>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, uint64_t time_step){
-		using dfi = df_info<FP,Descriptor>;
-
-		bool is_even = ((time_step % 2) == 0);
-		auto& cell = field(index);
-
-		auto& info = cell.template get<"info">();
-		if(info({0u}).get() == 0u){
-			return;
-		}
-		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">();
-
-		/**
-		 * Sum all known DFs
-		 */
-		saw::data<FP> sum_df{0};
-		for(saw::data<sch::UInt64> k{0u}; k < saw::data<sch::UInt64>{Descriptor::Q}; ++k){
-			auto c_k = dfi::directions[k.get()];
-			auto& cell_n = field({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}});
-			auto& info_n = cell_n.template get<"info">();
-			auto info_n_val = info_n({0u});
-			auto k_opp = dfi::opposite_index[k.get()];
-
-			if(info_n_val.get() > 0u){
-				sum_df += dfs_old({k_opp});
-			}
-		}
-		/**
-		 * Get the sum of the unknown dfs and precalculate the direction
-		 */
-		constexpr int known_dir = East ? 1 : -1;
-		auto sum_unknown_dfs = (rho_setting_ - sum_df) * saw::data<FP>{known_dir};
-
-		for(saw::data<sch::UInt64> k{0u}; k < saw::data<sch::UInt64>{Descriptor::Q}; ++k){
-			auto c_k = dfi::directions[k.get()];
-			auto& cell_n = field({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}});
-			auto& info_n = cell_n.template get<"info">();
-			auto info_n_val = info_n({0u});
-			// auto k_opp = dfi::opposite_index[k.get()];
-
-			if(info_n_val.get() > 0u){
-				sum_unknown_dfs += dfs_old({k}) * c_k[0u];
-			}
-		}
-
-		auto vel_x = sum_unknown_dfs / rho_setting_;
-
-		if constexpr (East) {
-			dfs_old({2u}) = dfs_old({1u}) + saw::data<FP>{2.0 / 3.0} * rho_setting_ * vel_x;
-			dfs_old({6u}) = dfs_old({5u}) + saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old({4u}) - dfs_old({3u}));
-			dfs_old({8u}) = dfs_old({7u}) + saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old({3u}) - dfs_old({4u}));
-		}else if constexpr (not East){
-			dfs_old({1u}) = dfs_old({2u}) - saw::data<FP>{2.0 / 3.0} * rho_setting_ * vel_x;
-			dfs_old({5u}) = dfs_old({6u}) - saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old({3u}) - dfs_old({4u}));
-			dfs_old({7u}) = dfs_old({8u}) - saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old({4u}) - dfs_old({3u}));
-		}
-	}
-};
-}
-}
 
-void set_geometry(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
+void set_geometry(
+	std::vector<saw::data<sch::Float64>, SyclDeviceAlloc<sch::Float64>>& dfs,
+	std::vector<saw::data<sch::Float64>, SyclDeviceAlloc<sch::Float64>>& dfs_old,
+	std::vector<saw::data<sch::UInt8>, SyclDeviceAlloc<sch::UInt8>>& info,
+	const saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>>& meta,
+	acpp::sycl::queue& sycl_q
+){
 	using namespace kel::lbm;
 
-	auto meta = latt.meta();
-
-	/**
-	 * Set ghost
-	 */
-	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
-
-		info({0u}).set(0u);
-
-	}, {{0u,0u}}, meta);
-
-	/**
-	 * Set wall
-	 */
-	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
-
-		info({0u}).set(2u);
-
-	}, {{0u,0u}}, meta, {{1u,1u}});
-
-	/**
-	 * Set fluid
-	 */
-	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
-
-		info({0u}).set(1u);
-
-	}, {{0u,0u}}, meta, {{2u,2u}});
-
-	/**
-	 * Set inflow
-	 */
-	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
-
-		info({0u}).set(3u);
-
-	}, {{1u,0u}}, {{2u,meta.at({1u})}}, {{0u,2u}});
-
-	/**
-	 * Set outflow
-	 */
-	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
-
-		info({0u}).set(4u);
-
-	}, {{meta.at({0u})-2u,0u}}, {{meta.at({0u})-1u, meta.at({1u})}}, {{0u,2u}});
 }
 
-void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
+void set_initial_conditions(
+	std::vector<saw::data<sch::Float64>, SyclDeviceAlloc<sch::Float64>>& dfs,
+	std::vector<saw::data<sch::Float64>, SyclDeviceAlloc<sch::Float64>>& dfs_old,
+	std::vector<saw::data<sch::UInt8>, SyclDeviceAlloc<sch::UInt8>>& info,
+	const saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>>& meta,
+	acpp::sycl::queue& sycl_q
+){
 	using namespace kel::lbm;
 
 	saw::data<sch::T> rho{1.0};
 	saw::data<sch::FixedArray<sch::T,sch::D2Q9::D>> vel{{0.0,0.0}};
 	auto eq = equilibrium<sch::T,sch::D2Q9>(rho, vel);
 
-	auto meta = latt.meta();
-
-	/**
-	 * Set distribution
-	 */
-	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(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::D2Q9::Q}; ++k){
-			dfs(k) = eq.at(k);
-			dfs_old(k) = eq.at(k);
-		}
-
-	}, {{0u,0u}}, meta);
 }
 
 void lbm_step(
-	saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt,
-	uint64_t time_step
+	std::vector<saw::data<sch::Float64>, SyclDeviceAlloc<sch::Float64>>& dfs,
+	std::vector<saw::data<sch::Float64>, SyclDeviceAlloc<sch::Float64>>& dfs_old,
+	std::vector<saw::data<sch::UInt8>, SyclDeviceAlloc<sch::UInt8>>& info,
+	const saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>>& meta,
+	uint64_t time_step,
+	acpp::sycl::queue& sycl_q
 ){
 	using namespace kel::lbm;
 	using dfi = df_info<sch::T,sch::D2Q9>;
@@ -246,55 +110,8 @@ void lbm_step(
 	component<sch::T, sch::D2Q9, cmpt::PressureBoundaryRestrictedVelocityTo<true>> inlet{1.01 * dfi::cs2};
 	component<sch::T, sch::D2Q9, cmpt::PressureBoundaryRestrictedVelocityTo<false>> outlet{1.0 * dfi::cs2};
 
-	auto meta = latt.meta();
-
-	/**
-	 * Collision
-	 */
-	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">();
-
-		switch(info({0u}).get()){
-			case 1u: {
-				coll.apply(latt, index, time_step);
-				break;
-			}
-			case 2u: {
-				bb.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::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">();
-
-		if(info_new({0u}).get() > 0u){
-			for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){
-				auto dir = dfi::directions[dfi::opposite_index[k]];
-				auto& cell_dir_old = latt({{index.at({0u})+dir[0],index.at({1u})+dir[1]}});
-
-				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);
+}
+}
 }
 
 template<typename T, typename Desc>
@@ -307,100 +124,51 @@ saw::error_or<void> kel_main(int argc, char** argv){
 	if(eo_lbm_dir.is_error()){
 		return std::move(eo_lbm_dir.get_error());
 	}
-
-	return saw::make_void();
-}
-
-int main(int argc, char** argv){
-	using namespace kel::lbm;
-	using dfi = df_info<sch::T,sch::D2Q9>;
-	auto eo_lbm_dir = output_directory();
-	if(eo_lbm_dir.is_error()){
-		return -1;
-	}
 	auto& lbm_dir = eo_lbm_dir.get_value();
-	auto out_dir = lbm_dir / "poiseulle_particles_channel_2d";
-
-	std::string_view cfg_file_name = "config.json";
-	if(argc > 1){
-		cfg_file_name = argv[1];
-	}
-
-	auto eo_conf = load_lbm_config<sch::Float64,sch::Descriptor<2u,9u>>(cfg_file_name);
-	if(eo_conf.is_error()){
-		auto& err = eo_conf.get_error();
-		std::cerr<<"[Error]: "<<err.get_category();
-		auto err_msg = err.get_message();
-		if(!err_msg.empty()){
-			std::cerr<<" - "<<err_msg;
-		}
-		std::cerr<<std::endl;
-
-		return err.get_id();
-	}
+	auto out_dir = lbm_dir / "poiseulle_channel_2d_gpu";
 
-	auto& conf = eo_conf.get_value();
+	// Create Dir TODO
 
-	converter<sch::Float64> conv {
-		{conf.template get<"delta_x">()},
-		{conf.template get<"delta_t">()}
+	// 
+	lbm::converter<sch::Float64> conv {
+		// delta_x
+		{{1.0}},
+		// delta_t
+		{{1.0}}
 	};
 
-	print_lbm_meta<sch::Float64,sch::Descriptor<2u,9u>>(conv, {conf.template get<"kinematic_viscosity">()});
+	uint64_t x_d = 1024u;
+	uint64_t y_d = 128u;
 
-	saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>> dim{{1024u, 128u}};
-	saw::data<sch::CavityFieldD2Q9, saw::encode::Native> lattice{dim};
-	auto meta = lattice.meta();
+	acpp::sycl::queue sycl_q;
+	SyclHostAlloc<sch::CellStruct> sycl_host_alloc{sycl_q};
+	SyclDeviceAlloc<sch::CellStruct> sycl_dev_alloc{sycl_q};
 
-	/**
-	 * Setup geometry
-	 */
-	set_geometry(lattice);
+	std::vector<saw::data<sch::CellStruct>, SyclHostAlloc<sch::CellStruct>> cells{x_d * y_d,sycl_alloc};
+	std::vector<saw::data<sch::CellStruct>, SyclDeviceAlloc<sch::CellStruct>> cells{x_d * y_d,sycl_alloc};
 
 	{
-
-		saw::data<sch::Array<sch::GeometryStruct<sch::T>, sch::D2Q9::D>> geo{dim};
-
-		iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-			auto& cell = lattice(index);
-			auto& info = cell.template get<"info">();
-
-			geo(index).template get<"info">().set(info({0u}).get());
-		}, {{0u,0u}}, dim);
-
-		std::string vtk_f_name{"tmp/geometry.vtk"};
-		write_vtk_file(vtk_f_name, geo);
+		auto eov = set_geometry(cells);
+		if(eov.is_error()){
+			return eov;
+		}
 	}
 
-	/**
-	 * Setup DFs
-	 */
-	set_initial_conditions(lattice);
-
-	saw::data<sch::Array<sch::MacroStruct<sch::T,sch::D2Q9::D>,sch::D2Q9::D>> macros{dim};
 
-	for(uint64_t i = 0u; i < 4096u*16u; ++i){
-		bool even_step = ((i % 2u) == 0u);
-
-		{
-			// Stream
-			iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-				auto& cell = lattice(index);
-				auto& dfs = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">();
-
-				auto& rho = macros.at(index).template get<"pressure">();
-				auto& vel = macros.at(index).template get<"velocity">();
-				compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel);
-				rho = rho * saw::data<sch::T>{dfi::cs2};
+	return saw::make_void();
+}
 
-			}, {{0u,0u}}, meta);
-			std::string vtk_f_name{"poiseulle_2d_"};
-			vtk_f_name += std::to_string(i) + ".vtk";
-			write_vtk_file(out_dir / vtk_f_name, macros);
+int main(int argc, char** argv){
+	auto eov = kel_main(argc, argv);
+	if(eov.is_error()){
+		auto& err = eov.get_error();
+		std::cerr<<"[Error] "<<err.get_category();
+		auto err_msg = err.get_message();
+		if(err_msg.size() > 0u){
+			std::cerr<<" - "<<err_msg;
 		}
-
-		lbm_step(lattice, i);
+		std::cerr<<std::endl;
+		return err.get_id();
 	}
-
 	return 0;
 }