Moving project structure for more less compilation

2 files changed, 236 insertions, 86 deletions

diff --git a/examples/cavity_2d_gpu/.nix/derivation.nix b/examples/cavity_2d_gpu/.nix/derivation.nix
index f8d7314..73b81b7 100644
--- a/examples/cavity_2d_gpu/.nix/derivation.nix
+++ b/examples/cavity_2d_gpu/.nix/derivation.nix
@@ -24,6 +24,8 @@ stdenv.mkDerivation {
     forstio.async
     forstio.codec
     forstio.codec-unit
+    forstio.remote
+    forstio.remote-sycl
 		forstio.codec-json
 		adaptive-cpp
 		kel-lbm
diff --git a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
index cc0a811..c4836d0 100644
--- a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
+++ b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp
@@ -1,7 +1,6 @@
 #include <kel/lbm/lbm.hpp>
 
 #include <forstio/codec/data.hpp>
-// #include <forstio/remote/
 
 #include <AdaptiveCpp/sycl/sycl.hpp>
 
@@ -9,6 +8,83 @@
 #include <fstream>
 #include <cmath>
 
+namespace saw {
+namespace encode {
+template<typename InnerEnc>
+struct Sycl {
+};
+}
+template<typename Schema>
+class data<Schema, encode::Sycl<encode::Native>> {
+private:
+	cl::sycl::buffer<data<Schema, encode::Native>> data_;
+	data<schema::UInt64, encode::Native> size_;
+public:
+	data(const data<Schema, encode::Native>& data__):
+		data_{&data__, 1u},
+		size_{data__.size()}
+	{}
+
+	auto& get_handle() {
+		return data_;
+	}
+
+	const auto& get_handle() const {
+		return data_;
+	}
+
+	data<schema::UInt64, encode::Native> size() const {
+		return size_;
+	}
+
+	template<cl::sycl::access::mode AccessMode>
+	auto access(cl::sycl::handler& h){
+		return data_.template get_access<AccessMode>(h);
+	}
+
+	template<cl::sycl::access::mode AccessMode>
+	auto access(cl::sycl::handler& h) const {
+		return data_.template get_access<AccessMode>(h);
+	}
+};
+
+template<typename Sch, uint64_t Dim>
+class data<schema::Array<Sch, Dim>, encode::Sycl<encode::Native>> {
+public:
+	using Schema = schema::Array<Sch,Dim>;
+private:
+	cl::sycl::buffer<data<Sch, encode::Native>> 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> size() const {
+		return size_;
+	}
+
+	template<cl::sycl::access::mode AccessMode>
+	auto access(cl::sycl::handler& h){
+		return data_.template get_access<AccessMode>(h);
+	}
+
+	template<cl::sycl::access::mode AccessMode>
+	auto access(cl::sycl::handler& h) const {
+		return data_.template get_access<AccessMode>(h);
+	}
+};
+}
+
 namespace kel {
 namespace lbm {
 namespace sch {
@@ -36,10 +112,13 @@ using CellInfo = Cell<UInt8, Desc, 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">
+using CellStruct = CellFieldStruct<
+	Desc,
+	Struct<
+		Member<CellField<Desc,DfCell<Desc>>, "dfs">,
+		Member<CellField<Desc,DfCell<Desc>>, "dfs_old">,
+		Member<CellField<Desc,CellInfo<Desc>>, "info">
+	>
 >;
 
 template<typename T, uint64_t D>
@@ -48,7 +127,7 @@ using MacroStruct = Struct<
 	Member<T, "pressure">
 >;
 
-using CavityFieldD2Q9 = CellField<D2Q9, CellStruct<D2Q9>>;
+using CavityFieldD2Q9 = CellStruct<D2Q9>;
 }
 
 
@@ -104,22 +183,20 @@ public:
 		*/
 	}
 };
-}
-}
 
 constexpr size_t dim_size = 2;
 constexpr size_t dim_x = 128;
 constexpr size_t dim_y = 128;
 
-void set_geometry(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
+void set_geometry(saw::data<sch::CavityFieldD2Q9>& lattice){
 	using namespace kel::lbm;
-	auto meta = latt.meta();
+	auto meta = lattice.meta();
+	auto& info_field = lattice.template get<"info">();
 	/**
 	 * Set ghost
 	 */
 	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
+		auto& info = info_field.at(index);
 
 		info({0u}).set(0u);
 
@@ -129,8 +206,7 @@ void set_geometry(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 	 * Set wall
 	 */
 	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
+		auto& info = info_field.at(index);
 
 		info({0u}).set(2u);
 
@@ -140,8 +216,7 @@ void set_geometry(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 	 * Set fluid
 	 */
 	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
+		auto& info = info_field.at(index);
 
 		info({0u}).set(1u);
 
@@ -151,30 +226,32 @@ void set_geometry(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 	 * Set top lid
 	 */
 	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
+		auto& info = info_field.at(index);
 
 		info({0u}).set(3u);
 
 	}, {{0u,1u}}, {{meta.at({0u}), 2u}}, {{2u,0u}});
 }
 
-void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
+void set_initial_conditions(
+	saw::data<sch::CavityFieldD2Q9>& lattice
+){
 	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();
+	auto meta = lattice.meta();
+	auto& dfs_field = lattice.template get<"dfs">();
+	auto& dfs_old_field = lattice.template get<"dfs_old">();
 
 	/**
 	 * 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">();
+		auto& dfs = dfs_field.at(index);
+		auto& dfs_old = dfs_old_field.at(index);
 
 		for(saw::data<sch::UInt64> k = 0; k < saw::data<sch::UInt64>{sch::D2Q9::Q}; ++k){
 			dfs(k) = eq.at(k);
@@ -184,84 +261,150 @@ void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 	}, {{0u,0u}}, meta);
 }
 
+struct sycl_component_context {
+	acpp::sycl::handler* handler = nullptr;
+};
+
 void lbm_step(
-	acpp::sycl::buffer<saw::data<sch::CavityFieldD2Q9>>& latt,
+	saw::data<sch::CellField<sch::D2Q9,sch::CellInfo<sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>>& info,
+	saw::data<sch::CellField<sch::D2Q9,kel::lbm::sch::DfCell  <sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>>& dfs,
+	saw::data<sch::CellField<sch::D2Q9,kel::lbm::sch::DfCell  <sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>>& dfs_old,
 	bool even_step,
 	uint64_t time_step,
 	acpp::sycl::queue& sycl_q
 ){
-	using namespace kel::lbm;
-	using dfi = df_info<sch::T,sch::D2Q9>;
+  using namespace kel::lbm;
+  using namespace acpp;
+
+  using dfi = df_info<sch::T, sch::D2Q9>;
+
+  //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;
+  bb_lid.lid_vel = {0.1, 0.0};
+
+  const size_t Nx = latt.template get_dim_size<0>().get();
+  const size_t Ny = latt.template get_dim_size<1>().get();
+
+  // Submit collision kernel
+  sycl_q.submit([&](sycl::handler& cgh) {
+    // Accessor for latt with read/write
+    auto info_acc = info.template get_access<sycl::access::mode::read>(cgh);
+    auto dfs_acc = dfs.template get_access<sycl::access::mode::read_write>(cgh);
+    auto dfs_old_acc = dfs_old.template get_access<sycl::access::mode::read_write>(cgh);
+
+    cgh.parallel_for(
+      sycl::range<2>{Nx, Ny}, [=](sycl::id<2> idx) {
+        size_t i = idx[0];
+        size_t j = idx[1];
+
+        // Read cell info
+        auto& info = info_acc[idx];
+
+        switch (info({0u}).get()) {
+          case 1u: {
+            // bb.apply(latt_acc, {i, j}, time_step);
+				// bool is_even = ((time_step % 2) == 0);
+
+				auto& dfs_old = (is_even) ? dfs_old_acc[idx] : dfs_acc[idx];
+				// auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
+
+				saw::data<T> rho;
+				saw::data<sch::FixedArray<T,Descriptor::D>> vel;
+				compute_rho_u<T,Descriptor>(dfs_old,rho,vel);
+				auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+				for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+					dfs_old({i}) = dfs_old({i}) + frequency_ * (eq.at(i) - dfs_old({i}));
+					// dfs_old({i}).set(dfs_old({i}).get() + (1.0 / relaxation_) * (eq.at(i).get() - dfs_old({i}).get()));
+				}
+            break;
+          }
+          case 2u: {
+            // coll.apply(latt_acc, {i, j}, time_step);
+				// bool is_even = ((time_step % 2) == 0);
 
-  /**
-	* 1. Relaxation parameter \tau
-	*/
-	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;
-	bb_lid.lid_vel = {0.1,0.0};
+				auto& dfs_old = (is_even) ? dfs_old_acc[idx] : dfs_acc[idx];
+				// auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
 
-	auto meta = latt.meta();
+				saw::data<T> rho;
+				saw::data<sch::FixedArray<T,Descriptor::D>> vel;
+				compute_rho_u<T,Descriptor>(dfs_old,rho,vel);
+				auto eq = equilibrium<T,Descriptor>(rho,vel);
 
-	/**
-	 * Collision
-	 */
-	iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
-		auto& cell = latt(index);
-		auto& info = cell.template get<"info">();
-
-		auto& dfs = cell.template get<"dfs">();
-		auto& dfs_old = cell.template get<"dfs_old">();
-
-		switch(info({0u}).get()){
-			case 1u: {
-				bb.apply(latt, index, time_step);
-				break;
-			}
-			case 2u: {
-				coll.apply(latt, index, time_step);
-				break;
-			}
-			default:
-				break;
-		}
+				using dfi = df_info<T,Descriptor>;
 
-	}, {{0u,0u}}, meta);
+				auto& force = cell.template get<"force">();
 
-	// Stream
-	for(uint64_t i = 1u; (i+1u) < latt.template get_dim_size<0>().get(); ++i){
-		for(uint64_t j = 1u; (j+1u) < latt.template get_dim_size<1>().get(); ++j){
-			auto& cell = latt({{i,j}});
-			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 && info_new({0u}).get() != 3u){
-				for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){
-					auto dir = dfi::directions[dfi::opposite_index[k]];
-					auto& cell_dir_old = latt({{i+dir[0],j+dir[1]}});
-
-					auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">();
-					auto& info_old = cell_dir_old.template get<"info">();
-
-					if( info_old({0}).get() == 3u ){
-						auto& df_old_loc = even_step ? latt({{i,j}}).template get<"dfs_old">() : latt({{i,j}}).template get<"dfs">();
-						df_new({k}) = df_old_loc({dfi::opposite_index.at(k)}) - 2.0 * dfi::inv_cs2 * dfi::weights.at(k) * 1.0 * ( bb_lid.lid_vel[0] * dir[0] + bb_lid.lid_vel[1] * dir[1]);
-						// dfs({dfi::opposite_index.at(i)}) = dfs_cpy({i}) - 2.0 * dfi::weights[i] * 1.0 * ( lid_vel[0] * dfi::directions[i][0] + lid_vel[1] * dfi::directions[i][1]) * dfi::inv_cs2;
-					} else {
-						df_new({k}) = df_old({k});
+				for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+					// saw::data<T> ci_min_u{0};
+					saw::data<T> ci_dot_u{0};
+					for(uint64_t d = 0u; d < Descriptor::D; ++d){
+						ci_dot_u = ci_dot_u + vel.at({d}) * saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])};
 					}
+					saw::data<T> F_i{0};// = saw::data<T>{dfi::weights[i]} * (ci_dot_F * dfi::inv_cs2 + ci_dot_u * ci_dot_F * (dfi::inv_cs2 * dfi::inv_cs2));
+
+					for(uint64_t d = 0u; d < Descriptor::D; ++d){
+						F_i = F_i +
+							saw::data<T>{dfi::weights[i]} * ((saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])}
+							- vel.at({d}) ) * dfi::inv_cs2 + ci_dot_u * saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])} * dfi::inv_cs2  * dfi::inv_cs2 ) * force({d});
+					}
+
+
+					dfs_old({i}) = dfs_old({i}) + frequency_ * (eq.at(i) - dfs_old({i}) ) + F_i * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_);
 				}
-			}
-		}
-	}
+            break;
+          }
+          default:
+            // Do nothing
+            break;
+        }
+      });
+  });
+
+  // Submit streaming kernel
+  sycl_q.submit([&](sycl::handler& cgh) {
+    auto info_acc = info.template get_access<sycl::access::mode::read_write>(cgh);
+
+    cgh.parallel_for(
+      sycl::range<2>{Nx - 2, Ny - 2}, [=](sycl::id<2> idx) {
+        size_t i = idx[0] + 1; // avoid boundary
+        size_t j = idx[1] + 1;
+
+        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 && info_new({0u}).get() != 3u) {
+          for (uint64_t k = 0u; k < sch::D2Q9::Q; ++k) {
+            auto dir = dfi::directions[dfi::opposite_index[k]];
+            auto& cell_dir_old = latt_acc({i + dir[0], j + dir[1]});
+
+            auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">();
+            auto& info_old = cell_dir_old.template get<"info">();
+
+            if (info_old({0}).get() == 3u) {
+              auto& df_old_loc = even_step ? latt_acc({i, j}).template get<"dfs_old">() : latt_acc({i, j}).template get<"dfs">();
+              df_new({k}) = df_old_loc({dfi::opposite_index.at(k)}) - 2.0 * dfi::inv_cs2 * dfi::weights.at(k) * 1.0 * (bb_lid.lid_vel[0] * dir[0] + bb_lid.lid_vel[1] * dir[1]);
+            } else {
+              df_new({k}) = df_old({k});
+            }
+          }
+        }
+      });
+  });
+
+  sycl_q.wait();
+}
+}
 }
 
 int main(){
 	using namespace kel::lbm;
 
-	saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>> dim{{dim_x, dim_y}};
+	saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>> dim{{dim_x,info_field.at dim_y}};
+	const dim_size = dim_x * dim_y;
 
-	saw::data<sch::CavityFieldD2Q9, saw::encode::Native> lattice{dim};
+	saw::data<sch::CavityFieldD2Q9>> lattice{dim};
 
 	converter<sch::T> conv{
 		{0.1},
@@ -287,13 +430,18 @@ int main(){
 	 */
 	set_initial_conditions(lattice);
 
+	auto& info_field = lattice.template get<"info">();
+	auto& dfs_field = lattice.template get<"dfs">();
+	auto& dfs_old_field = lattice.template get<"dfs_old">();
+
 	acpp::sycl::queue sycl_q{acpp::sycl::default_selector_v, acpp::sycl::property::queue::in_order{}};
-	acpp::sycl::buffer<saw::data<sch::CavityFieldD2Q9>> df_sycl{&lattice, 1u};
+	saw::data<sch::CellField<Desc,sch::CellInfo<sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>> sycl_info_field{info_field};
+	saw::data<sch::CellField<Desc,sch::DfCell<sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>> sycl_dfs_field{dfs_field};
+	saw::data<sch::CellField<Desc,sch::DfCell<sch::D2Q9>>, saw::encode::Sycl<saw::encode::Native>> sycl_dfs_old_field{dfs_old_field};
 
 	/**
 	 * Timeloop
 	 */
-
 	uint64_t lattice_steps = 512000u;
 	bool even_step = true;
 
@@ -310,7 +458,7 @@ int main(){
 			write_vtk_file(vtk_f_name, macros);
 		}
 
-		lbm_step(df_sycl, (i%2u == 0u), i, sycl_q);
+		lbm_step(sycl_info_field, sycl_dfs_field, sycl_dfs_old_field, (i%2u == 0u), i, sycl_q);
 	}
 	return 0;
 }