Streamlining and thoughts about changes in particle dynamicsdev

5 files changed, 108 insertions, 17 deletions

diff --git a/c++/descriptor.hpp b/c++/descriptor.hpp
index 0c6707a..d04d217 100644
--- a/c++/descriptor.hpp
+++ b/c++/descriptor.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <forstio/codec/data.hpp>
+#include <forstio/codec/data_math.hpp>
 #include <forstio/codec/schema_factory.hpp>
 
 namespace kel {
diff --git a/c++/macroscopic.hpp b/c++/macroscopic.hpp
index 8a25248..51368e9 100644
--- a/c++/macroscopic.hpp
+++ b/c++/macroscopic.hpp
@@ -59,5 +59,34 @@ void compute_rho_u (
 		vel().at({i}) = vel().at({i}) / rho();
 	}
 }
+
+/**
+ * Calculate the macroscopic variables rho and u in Lattice Units.
+ */
+template<typename T, typename Desc>
+void compute_rho_u (
+		const saw::data<sch::Cell<T, Desc, 0, 0, 1>>& dfs,
+		saw::ref<saw::data<T>> rho,
+		saw::ref<saw::data<sch::Vector<T,Desc::D>>> vel
+	)
+{
+	using dfi = df_info<T, Desc>;
+
+	rho().set(0);
+	for(uint64_t i = 0; i < Desc::D; ++i){
+		vel().at({{i}}).set(0);
+	}
+
+	for(size_t j = 0; j < Desc::Q; ++j){
+		rho() = rho() + dfs(j);
+		for(size_t i = 0; i < Desc::D; ++i){
+			vel().at({{i}}) = vel().at({{i}}) + saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[j][i])} * dfs(j);
+		}
+	}
+
+	for(size_t i = 0; i < Desc::D; ++i){
+		vel().at({{i}}) = vel().at({{i}}) / rho();
+	}
+}
 }
 }
diff --git a/c++/util.hpp b/c++/util.hpp
index 441ffc8..c5b6a05 100644
--- a/c++/util.hpp
+++ b/c++/util.hpp
@@ -4,6 +4,8 @@
 #include <forstio/codec/data.hpp>
 
 #include <iostream>
+#include <iomanip>
+#include <sys/ioctl.h>
 
 namespace kel {
 namespace lbm {
@@ -46,18 +48,37 @@ public:
 };
 */
 
-void print_progress_bar(uint32_t progress, uint32_t progress_target){
-	std::cout
-		<<"\r"
-		<<"Progress: "
-		<<((100 * progress) / progress_target)
-		<<"% [";
+void print_progress_bar(const uint32_t progress, const uint32_t progress_target){
+	std::cout<<"\r";
+	//	<<"Progress: "
+	//	<<((100 * progress) / progress_target)
+	//	<<"% [";
 
-	uint64_t i{0u};
+	uint64_t perc_prog = (100ul*progress) / progress_target;
+	constexpr uint64_t max_perc_progress = 100u;
+	perc_prog = std::min(perc_prog, max_perc_progress);
 
-	uint64_t max_display = 64u;
+	std::string_view prog_str{"Progress: "};
+	// Progress string + (100 width and perc char) + ([] brackets) + space + pointer
+	uint64_t i{prog_str.size() + 4u + 2u + 1u + 1u};
+
+	std::cout<<prog_str;
+	std::cout<<std::setw(3u)<<perc_prog<<"%";
+	std::cout<<" ";
+	std::cout<<"[";
+
+	uint64_t max_display = []() -> uint64_t{
+		struct winsize w;
+		if(ioctl(STDOUT_FILENO, TIOCGWINSZ,&w) == -1){
+			// Standardized Terminal size
+			return 80u;
+		}
+		return w.ws_col;
+	}();
+	max_display = std::max(max_display,i) - i;
 	uint64_t progress_display = (max_display * progress) / progress_target;
-	for(; i < progress_display; ++i){
+
+	for(i = 0u; i < progress_display; ++i){
 		std::cout<<"#";
 	}
 	for(; i < max_display; ++i){
@@ -65,6 +86,7 @@ void print_progress_bar(uint32_t progress, uint32_t progress_target){
 	}
 
 	std::cout<<"]";
+
 	std::cout<<std::flush;
 }
 }
diff --git a/c++/write_vtk.hpp b/c++/write_vtk.hpp
index c9ecc99..0647db5 100644
--- a/c++/write_vtk.hpp
+++ b/c++/write_vtk.hpp
@@ -3,6 +3,7 @@
 #include <forstio/error.hpp>
 
 #include <forstio/codec/data.hpp>
+#include <forstio/codec/data_math.hpp>
 
 #include "descriptor.hpp"
 
@@ -55,6 +56,31 @@ struct lbm_vtk_writer<sch::FixedArray<T,D>> {
 	}
 };
 
+template<typename T, uint64_t D>
+struct lbm_vtk_writer<sch::Vector<T,D>> {
+	static saw::error_or<void> apply_header(std::ostream& vtk_file, std::string_view name){
+		vtk_file<<"VECTORS "<<name<<" float\n";
+		return saw::make_void();
+	}
+	static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::Vector<T,D>>& field){
+		static_assert(D > 0, "Non-dimensionality is bad for velocity.");
+		static_assert(D <= 3, "4th dimension as well. Mostly due to vtk.");
+
+		// vtk_file<<"VECTORS "<<name<<" float\n";
+		for(uint64_t i = 0u; i < D; ++i){
+			if(i > 0){
+				vtk_file<<" ";
+			}
+			vtk_file<<field.at({{i}}).get();
+		}
+		for(uint64_t i = D; i < 3; ++i){
+			vtk_file<<" 0";
+		}
+		vtk_file<<"\n";
+		return saw::make_void();
+	}
+};
+
 template<uint64_t Dim, typename... StructT, saw::string_literal... StructN>
 struct lbm_vtk_writer<sch::Array<sch::Struct<sch::Member<StructT,StructN>...> , Dim>> {
 
diff --git a/examples/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d.cpp
index 9f7475d..4673a2a 100644
--- a/examples/poiseulle_particles_channel_2d.cpp
+++ b/examples/poiseulle_particles_channel_2d.cpp
@@ -54,10 +54,11 @@ using CellStruct = Struct<
 
 template<typename T, uint64_t D>
 using MacroStruct = Struct<
-	Member<FixedArray<T,D>, "velocity">,
+	Member<Vector<T,D>, "velocity">,
 	Member<T, "pressure">,
 	Member<UInt16, "particle">,
 	Member<FixedArray<T,D>, "force">
+	//Member<Vector<T,D>, "debug_acceleration">
 >;
 
 template<typename T>
@@ -348,13 +349,16 @@ void couple_particles_to_lattice(
 
 			// p_pos - p_pos_old
 			// auto p_vel_rel = closest_u - p_vel;
-			saw::data<sch::Vector<sch::T, 2u>> p_vel_rel;
-			p_vel_rel.at({{0u}}) = closest_u.at({0u}) - p_vel.at({{0u}});
-			p_vel_rel.at({{1u}}) = closest_u.at({1u}) - p_vel.at({{1u}});
+			saw::data<sch::Vector<sch::T, 2u>> p_vel_rel = closest_u - p_vel;
 
+			// Add relative velocity to particle acceleration (Since our timestep is 1 we don't need to do anything else)
+			p_acc = p_acc + p_vel_rel;
+
+			/*
 			for(saw::data<sch::UInt64> i{0u}; i.get() < 2u; ++i){
 				p_acc.at({{i}}) = p_acc.at({{i}}) + p_vel_rel.at({{i}});
 			}
+			*/
 
 			// Add forces to put away from walls
 			/// 1. Check if neighbour is wall
@@ -366,10 +370,11 @@ void couple_particles_to_lattice(
 
 				auto& n_cell = latt(n_p_cell_pos);
 				auto& n_info = n_cell.template get<"info">()({0u});
-				auto& n_macro_cell_particle = macros.at(n_p_cell_pos).template get<"particle">();
+				auto& n_macro_cell = macros.at(n_p_cell_pos);
+				auto& n_macro_cell_particle = n_macro_cell.template get<"particle">();
 
 				// If neighbour is wall, then add force pushing the particle away
-				if(n_info.get() <= 1u or (n_macro_cell_particle.get() != i.get() and n_macro_cell_particle.get() > 0u) ) {
+				if(n_info.get() <= 1u or (n_macro_cell_particle.get() != (i.get()+1u) and n_macro_cell_particle.get() > 0u) ) {
 					// add to p_acc
 
 					saw::data<sch::T> dist_dot{
@@ -379,8 +384,16 @@ void couple_particles_to_lattice(
 
 					if(dist_dot.get() > 0){
 						// TODO add if particle is close
-						p_acc.at({{0u}}) = p_acc.at({{0u}}) + saw::data<sch::Int32>{100 * dfi::directions[dfi::opposite_index[k.get()]][0u]}.template cast_to<sch::T>();
-						p_acc.at({{1u}}) = p_acc.at({{1u}}) + saw::data<sch::Int32>{100 * dfi::directions[dfi::opposite_index[k.get()]][1u]}.template cast_to<sch::T>();
+						//auto& n_debug_acc = n_macro_cell.template get<"debug_acceleration">();
+
+						saw::data<sch::Vector<sch::T,sch::D2Q9::D>> push_force;
+						{
+							push_force.at({{0u}}) = saw::data<sch::Int32>{1000 * dfi::directions[dfi::opposite_index[k.get()]][0u]}.template cast_to<sch::T>();
+							push_force.at({{1u}}) = saw::data<sch::Int32>{1000 * dfi::directions[dfi::opposite_index[k.get()]][1u]}.template cast_to<sch::T>();
+						}
+
+						//n_debug_acc = n_debug_acc + push_force;
+						p_acc = p_acc + push_force;
 					}
 				}
 			}