Dangling changes

Fixed ForcedGuo collision

9 files changed, 54 insertions, 36 deletions

diff --git a/c++/boundary.hpp b/c++/boundary.hpp
index 3cf6f3f..bc4859c 100644
--- a/c++/boundary.hpp
+++ b/c++/boundary.hpp
@@ -51,7 +51,7 @@ public:
 		auto& cell = field(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">();
+		// auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
 
 		using dfi = df_info<T,Descriptor>;
 
diff --git a/c++/collision.hpp b/c++/collision.hpp
index df9a734..cd6cf14 100644
--- a/c++/collision.hpp
+++ b/c++/collision.hpp
@@ -17,11 +17,11 @@ template<typename T, typename Descriptor>
 class component<T, Descriptor, cmpt::BGK> {
 private:
 	typename saw::native_data_type<T>::type relaxation_;
-	typename saw::native_data_type<T>::type frequency_;
+	saw::data<T> frequency_;
 public:
 	component(typename saw::native_data_type<T>::type relaxation__):
 		relaxation_{relaxation__},
-		frequency_{1 / relaxation_}
+		frequency_{typename saw::native_data_type<T>::type(1) / relaxation_}
 	{}
 
 	using Component = cmpt::BGK;
@@ -57,7 +57,8 @@ public:
 		auto eq = equilibrium<T,Descriptor>(rho,vel);
 
 		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
-			dfs_old({i}).set(dfs_old({i}).get() + (1.0 / relaxation_) * (eq.at(i).get() - dfs_old({i}).get()));
+			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()));
 		}
 	}
 };
@@ -66,11 +67,11 @@ template<typename T, typename Descriptor>
 class component<T, Descriptor, cmpt::BGKGuo> {
 private:
 	typename saw::native_data_type<T>::type relaxation_;
-	typename saw::native_data_type<T>::type frequency_;
+	saw::data<T> frequency_;
 public:
 	component(typename saw::native_data_type<T>::type relaxation__):
 		relaxation_{relaxation__},
-		frequency_{1 / relaxation_}
+		frequency_{typename saw::native_data_type<T>::type(1) / relaxation_}
 	{}
 
 	using Component = cmpt::BGKGuo;
@@ -89,7 +90,7 @@ public:
 
 	template<typename CellFieldSchema>
 	void apply(saw::data<CellFieldSchema>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, uint64_t time_step){
-		bool is_even = ((time_step & 2) == 0);
+		bool is_even = ((time_step % 2) == 0);
 		auto& cell = field(index);
 
 		auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
@@ -105,22 +106,21 @@ public:
 		auto& force = cell.template get<"force">();
 
 		for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+			saw::data<T> ci_min_u{0};
 			saw::data<T> ci_dot_u{0};
-			saw::data<T> ci_dot_F{0};
-			saw::data<T> u_dot_F{0};
 			for(uint64_t d = 0u; d < Descriptor::D; ++d){
-				ci_dot_u = ci_dot_u + vel.at({d}) * dfi::directions[i][d];
-				ci_dot_F = ci_dot_F + force({d}) * dfi::directions[i][d];
-				u_dot_F = u_dot_F + vel.at({d}) * force({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));
 
-			// auto df_old = dfs_old({i});
-			saw::data<T> F_i = saw::data<T>{dfi::weights[i]} * (ci_dot_F / dfi::cs2 + ci_dot_u * ci_dot_F / (dfi::cs2 * dfi::cs2));
-			if(F_i.get() != 0){
-				std::cerr<<"AAAH";
+			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}) + saw::data<T>{frequency_} * (eq.at(i) - dfs_old({i}) ) - F_i * saw::data<T>{1.0f - 0.5f * frequency_};
+
+			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_);
 		}
 	}
 };
diff --git a/c++/particle/particle.hpp b/c++/particle/particle.hpp
index a55b06e..ae7769a 100644
--- a/c++/particle/particle.hpp
+++ b/c++/particle/particle.hpp
@@ -27,7 +27,7 @@ using ParticleMask = Struct<
 template<typename T, uint64_t D>
 using Particle = Struct<
 	Member<ParticleRigidBody<T,D>, "rigid_body">,
-	Member<ParticleMask<Float32,D>, "mask">,
+	Member<ParticleMask<T,D>, "mask">,
 	Member<T, "size">
 >;
 }
@@ -84,6 +84,7 @@ public:
 
 	template<typename LbmLattice>
 	void update_particle_border(saw::data<LbmLattice>& latt){
+		(void) latt;
 		for(auto& iter : particles_){
 			auto& par = iter;
 
diff --git a/c++/statistics.hpp b/c++/statistics.hpp
new file mode 100644
index 0000000..c07ccb7
--- /dev/null
+++ b/c++/statistics.hpp
@@ -0,0 +1,11 @@
+#pragma once
+
+namespace kel {
+namespace lbm {
+template<typename T>
+class statistics {
+private:
+public:
+};
+}
+}
diff --git a/c++/write_vtk.hpp b/c++/write_vtk.hpp
index eb2545a..c9ecc99 100644
--- a/c++/write_vtk.hpp
+++ b/c++/write_vtk.hpp
@@ -81,7 +81,7 @@ struct lbm_vtk_writer<sch::Array<sch::Struct<sch::Member<StructT,StructN>...> ,
 	template<uint64_t i>
 	static saw::error_or<void> write_i(std::ostream& vtk_file, const saw::data<sch::Array<sch::Struct<sch::Member<StructT,StructN>...>,Dim>>& field){
 
-		auto meta = field.get_dims();
+		// auto meta = field.get_dims();
 		saw::data<sch::FixedArray<sch::UInt64,Dim>> index;
 		for(saw::data<sch::UInt64> it{0}; it.get() < Dim; ++it){
 			index.at({0u}).set(0u);
diff --git a/examples/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu.cpp
index 19563e2..b5c9b7f 100644
--- a/examples/cavity_2d_gpu.cpp
+++ b/examples/cavity_2d_gpu.cpp
@@ -189,7 +189,8 @@ void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
 
 void lbm_step(
 	saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt,
-	uint64_t time_step
+	uint64_t time_step,
+	sycl::queue& sycl_q
 ){
 	using namespace kel::lbm;
 	using dfi = df_info<sch::T,sch::D2Q9>;
@@ -260,8 +261,6 @@ void lbm_step(
 int main(){
 	using namespace kel::lbm;
 
-	saw::remote<saw::rmt::Sycl> sycl_rmt;
-
 	saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>> dim{{dim_x, dim_y}};
 
 	saw::data<sch::CavityFieldD2Q9, saw::encode::Native> lattice{dim};
@@ -273,6 +272,13 @@ int main(){
 
 	print_lbm_meta<sch::T, sch::D2Q9>(conv, {1e-3});
 
+	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 / "cavity_gpu_2d";
+
 	//auto& df_field = lattices.at(0).template get<"dfs">();
 	//for(uint64_t i = 0; i < df_field.get_dim_size<0u>(); ++i){
 	//	lattices.at(i) = {dim_x, dim_y};
@@ -288,6 +294,8 @@ int main(){
 	 */
 	set_initial_conditions(lattice);
 
+	sycl::queue sycl_q{sycl::default_selector_v, sycl::property::queue::in_order{}};
+
 	/**
 	 * Timeloop
 	 */
@@ -307,7 +315,7 @@ int main(){
 			write_vtk_file(vtk_f_name, macros);
 		}
 
-		// lbm_step(lattice, i);
+		lbm_step(lattice, i, sycl_q);
 	}
 	return 0;
 }
diff --git a/examples/poiseulle_2d.cpp b/examples/poiseulle_2d.cpp
index 8e1ba94..8c31cb8 100644
--- a/examples/poiseulle_2d.cpp
+++ b/examples/poiseulle_2d.cpp
@@ -93,7 +93,7 @@ public:
 			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">();
+		// auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
 
 		/**
 		 * Sum all known DFs
@@ -121,7 +121,7 @@ public:
 			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()];
+			// auto k_opp = dfi::opposite_index[k.get()];
 
 			if(info_n_val.get() > 0u){
 				sum_unknown_dfs += dfs_old({k}) * c_k[0u];
@@ -255,9 +255,6 @@ void lbm_step(
 		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: {
 				coll.apply(latt, index, time_step);
diff --git a/examples/poiseulle_channel_2d.cpp b/examples/poiseulle_channel_2d.cpp
index d34703a..ffc7201 100644
--- a/examples/poiseulle_channel_2d.cpp
+++ b/examples/poiseulle_channel_2d.cpp
@@ -93,7 +93,7 @@ public:
 			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">();
+		//auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
 
 		/**
 		 * Sum all known DFs
@@ -121,7 +121,7 @@ public:
 			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()];
+			//auto k_opp = dfi::opposite_index[k.get()];
 
 			if(info_n_val.get() > 0u){
 				sum_unknown_dfs += dfs_old({k}) * c_k[0u];
@@ -294,8 +294,8 @@ void lbm_step(
 		auto& cell = latt(index);
 		auto& info = cell.template get<"info">();
 
-		auto& dfs = cell.template get<"dfs">();
-		auto& dfs_old = cell.template get<"dfs_old">();
+		//auto& dfs = cell.template get<"dfs">();
+		//auto& dfs_old = cell.template get<"dfs_old">();
 
 		switch(info({0u}).get()){
 			case 1u: {
diff --git a/examples/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d.cpp
index fed78b2..bc5d241 100644
--- a/examples/poiseulle_particles_channel_2d.cpp
+++ b/examples/poiseulle_particles_channel_2d.cpp
@@ -24,6 +24,7 @@ using namespace saw::schema;
  * Q factor
  */
 using T = Float32;
+// using T = MixedPrecision<Float64, Float32>;
 using D2Q5 = Descriptor<2u,5u>;
 using D2Q9 = Descriptor<2u,9u>;
 
@@ -349,6 +350,7 @@ void couple_particles_to_lattice(
 		saw::data<kel::lbm::sch::T> time_step
 ){
 	using namespace kel::lbm;
+	(void) time_step;
 
 	auto meta = latt.meta();
 	using dfi = df_info<sch::T,sch::D2Q9>;
@@ -398,7 +400,7 @@ void couple_particles_to_lattice(
 
 		// Fluid to Particle Coupling
 		// Prepare force sum
-		saw::data<sch::FixedArray<sch::T,2u>> forces;
+		saw::data<sch::Vector<sch::T,2u>> forces;
 
 		iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){
 
@@ -500,8 +502,8 @@ void lbm_step(
 		auto& cell = latt(index);
 		auto& info = cell.template get<"info">();
 
-		auto& dfs = cell.template get<"dfs">();
-		auto& dfs_old = cell.template get<"dfs_old">();
+		//auto& dfs = cell.template get<"dfs">();
+		//auto& dfs_old = cell.template get<"dfs_old">();
 
 		switch(info({0u}).get()){
 			case 1u: {
@@ -595,7 +597,6 @@ int main(int argc, char** argv){
 	particle_system<sch::T, 2u> particle_sys;
 	add_particles(particle_sys);
 
-
 	/**
 	 * Setup geometry
 	 */