End of day, attempts at utilization

4 files changed, 458 insertions, 6 deletions

diff --git a/c++/util.hpp b/c++/util.hpp
index a215cc1..1a54541 100644
--- a/c++/util.hpp
+++ b/c++/util.hpp
@@ -2,11 +2,22 @@
 
 namespace kel {
 namespace lbm {
-template<typename T, typename Descriptor>
+template<typename T, typename Descriptor, typename CellFieldSchema>
 struct is_neighbour {
 	template<typename CellFieldSchema>
-	static bool operator()(saw::data<CellFieldSchema>& latt,  saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>& index, saw::data<sch::UInt64> range) {
+	static bool operator()(saw::data<CellFieldSchema>& latt,  saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>& index) {
 		using dfi = df_info<T,Descriptor>;
+
+		if(index.at({0u}).get() == 0u or index.at({1u}).get() == 0u){
+			return false;
+		}
+
+		for(saw::data<sch::UInt64> k{0u}; k.get() < Descriptor::Q; ++k){
+			// TODO
+			saw::data<sch::FixedArray<sch::UInt64,2u>>
+		}
+
+		auto& cell = latt(index);
 	}
 };
 }
diff --git a/examples/SConscript b/examples/SConscript
index e6eaaf4..13335c8 100644
--- a/examples/SConscript
+++ b/examples/SConscript
@@ -28,6 +28,11 @@ examples_objects = [];
 examples_env.add_source_files(examples_objects, ['cavity_2d.cpp'], shared=False);
 examples_env.cavity_2d = examples_env.Program('#bin/cavity_2d', [env.library_static, examples_objects]);
 
+# Cavity3D
+examples_objects = [];
+examples_env.add_source_files(examples_objects, ['cavity_3d.cpp'], shared=False);
+examples_env.cavity_3d = examples_env.Program('#bin/cavity_3d', [env.library_static, examples_objects]);
+
 # Channel Throat 2D
 #examples_objects = [];
 #examples_env.add_source_files(examples_objects, ['particle_ibm.cpp'], shared=False);
@@ -51,7 +56,8 @@ examples_env.poiseulle_particles_channel_2d = examples_env.Program('#bin/poiseul
 # Set Alias
 env.examples = [
 	examples_env.meta_2d,
-	examples_env.cavity_2d,
+#	examples_env.cavity_2d,
+#	examples_env.cavity_3d,
 #	examples_env.particle_ibm_2d
 	examples_env.poiseulle_2d,
 	examples_env.poiseulle_channel_2d,
diff --git a/examples/cavity_3d.cpp b/examples/cavity_3d.cpp
index e69de29..64356be 100644
--- a/examples/cavity_3d.cpp
+++ b/examples/cavity_3d.cpp
@@ -0,0 +1,436 @@
+#include "../c++/descriptor.hpp"
+#include "../c++/macroscopic.hpp"
+#include "../c++/lbm.hpp"
+#include "../c++/component.hpp"
+#include "../c++/collision.hpp"
+#include "../c++/boundary.hpp"
+
+#include <forstio/codec/data.hpp>
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+
+namespace kel {
+namespace lbm {
+namespace sch {
+using namespace saw::schema;
+
+/**
+ * Basic distribution function
+ * Base type
+ * D
+ * Q
+ * Scalar factor
+ * D factor
+ * Q factor
+ */
+using T = Float32;
+using D3Q27 = Descriptor<3u,27u>;
+
+template<typename Desc>
+using DfCell = Cell<T, Desc, 0u, 0u, 1u>;
+
+template<typename Desc>
+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">
+>;
+
+template<typename T, uint64_t D>
+using MacroStruct = Struct<
+	Member<FixedArray<T,D>, "velocity">,
+	Member<T, "pressure">
+>;
+
+using CavityFieldD3Q27 = CellField<D3Q27, CellStruct<D3Q27>>;
+}
+
+
+/*
+template<typename T, typename Encode>
+class df_cell_view;
+*/
+/**
+ * Minor helper for the AA-Pull Pattern, so I can use only one lattice
+ *
+ * Am I sure I want to use AA this way?
+ * Esoteric Twist technically reduces the needed memory access footprint
+ */
+/*
+template<typename Desc, size_t SN, size_t DN, size_t QN, typename Encode>
+class df_cell_view<sch::Cell<sch::T, Desc, SN, DN, QN>, Encode> {
+public:
+	using Schema = sch::Cell<sch::T,Desc,SN,DN,QN>;
+private:
+		std::array<std::decay_t<typename saw::native_data_type<sch::T>::type>*, QN> view_;
+public:
+		df_cell_view(const std::array<std::decay_t<typename saw::native_data_type<sch::T>::type>*, QN>& view):
+				view_{view}
+		{}
+};
+*/
+namespace cmpt {
+struct MovingWall {};
+}
+
+/**
+ * Full-Way moving wall Bounce back, something is not right here.
+ * Technically it should reflect properly.
+ */
+template<typename Desc>
+class component<sch::T, Desc, cmpt::MovingWall> {
+public:
+	std::array<typename saw::native_data_type<sch::T>::type, Desc::D> lid_vel;
+
+public:
+	void apply(
+		saw::data<sch::DfCell<Desc>>& dfs
+	){
+		using dfi = df_info<sch::T,Desc>;
+
+		// Technically use .copy()
+		/*
+		auto dfs_cpy = dfs;
+
+		for(uint64_t i = 0u; i < Desc::Q; ++i){
+			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;
+		}
+		*/
+	}
+};
+}
+}
+
+constexpr size_t dim_size = 3;
+constexpr size_t dim_x = 128;
+constexpr size_t dim_y = 128;
+constexpr size_t dim_z = 128;
+
+
+template<typename Func>
+void apply_for_cells(Func&& func, saw::data<kel::lbm::sch::CavityFieldD3Q27>& dat){
+	for(std::size_t i = 0; i < dat.meta().at({1u}).get(); ++i){
+		for(std::size_t j = 0; j < dat.meta().at({0u}).get(); ++j){
+			saw::data<saw::schema::UInt64> di{i};
+			saw::data<saw::schema::UInt64> dj{j};
+			auto& cell_v = dat({{dj,di}});
+			func(cell_v, j, i);
+		}
+	}
+}
+
+void set_geometry(saw::data<kel::lbm::sch::CavityFieldD3Q27>& latt){
+	using namespace kel::lbm;
+	apply_for_cells([](auto& cell, std::size_t i, std::size_t j){
+		uint8_t val = 0;
+		if(j == 1){
+			val = 2u;
+		}
+		if(i == 1 || (i+2) == dim_x || (j+2) == dim_y){
+			val = 3u;
+		}
+		if(i > 1 && (i+2) < dim_x && j > 1 && (j+2) < dim_y){
+			val = 1u;
+		}
+		if(i == 0 || j == 0 || (i+1) == dim_x || (j+1) == dim_y){
+			val = 0u;
+		}
+		cell.template get<"info">()(0u).set(val);
+	}, latt);
+}
+
+void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
+	using namespace kel::lbm;
+
+	saw::data<sch::T> rho{1.0};
+	{
+		saw::data<sch::FixedArray<sch::T,sch::D3Q27::D>> vel{{0.0,0.0}};
+		auto eq = equilibrium<sch::T,sch::D3Q27>(rho, vel);
+
+		apply_for_cells([&eq](auto& cell, std::size_t i, std::size_t j){
+			(void) i;
+			(void) j;
+			auto& dfs = cell.template get<"dfs">();
+			auto& dfs_old = cell.template get<"dfs_old">();
+			auto info = cell.template get<"info">()(0u).get();
+			for(uint64_t k = 0; k < sch::D3Q27::Q; ++k){
+				dfs(k) = eq.at({k});
+				dfs_old(k) = eq.at({k});
+			}
+		}, latt);
+	}
+
+	{
+		saw::data<sch::FixedArray<sch::T,sch::D3Q27::D>> vel{{0.1,0.0}};
+		auto eq = equilibrium<sch::T,sch::D3Q27>(rho, vel);
+
+		apply_for_cells([&eq](auto& cell, std::size_t i, std::size_t j){
+			(void) i;
+			(void) j;
+			auto& dfs = cell.template get<"dfs">();
+			auto& dfs_old = cell.template get<"dfs_old">();
+			auto info = cell.template get<"info">()(0u).get();
+			if(info == 2u){
+				for(uint64_t k = 0; k < sch::D3Q27::Q; ++k){
+					dfs(k) = eq.at({k});
+					dfs_old(k) = eq.at({k});
+				}
+			}
+		}, latt);
+	}
+}
+
+void lbm_step(
+	saw::data<kel::lbm::sch::CavityFieldD3Q27>& latt,
+	bool even_step,
+	uint64_t time_step
+){
+	using namespace kel::lbm;
+  using dfi = df_info<sch::T,sch::D3Q27>;
+
+  /**
+	* 1. Relaxation parameter \tau
+	*/
+	component<sch::T, sch::D3Q27, cmpt::BGK> coll{0.5384};
+	component<sch::T, sch::D3Q27, cmpt::BounceBack> bb;
+
+	component<sch::T, sch::D3Q27, cmpt::MovingWall> bb_lid;
+	bb_lid.lid_vel = {0.1,0.0,0.0};
+
+	auto meta = latt.meta();
+
+	// Collide
+	for(saw::data<sch::UInt64> i{0u}; i < meta.at(0u); ++i){
+		for(saw::data<sch::UInt64> j{0u}; j < meta.at(1u); ++j ){
+			auto& cell = latt({{i,j}});
+			auto& info = cell.template get<"info">();
+
+			switch(info({0u}).get()){
+				case 1u: {
+					coll.apply(latt, {{i,j}}, time_step);
+					break;
+				}
+				case 2u: {
+					auto& df = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
+					bb_lid.apply(df);
+					break;
+				}
+				case 3u: {
+					bb.apply(latt, {{i,j}}, time_step);
+					break;
+				}
+				default:
+					break;
+			}
+		}
+	}
+
+	// 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() != 2u){
+				for(uint64_t k = 0u; k < sch::D3Q27::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() == 2u ){
+						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});
+					}
+				}
+			}
+		}
+	}
+}
+
+int main(){
+	using namespace kel::lbm;
+
+	saw::data<sch::FixedArray<sch::UInt64,sch::D3Q27::D>> dim{{dim_x, dim_y}};
+
+	saw::data<sch::CavityFieldD3Q27, saw::encode::Native> lattice{dim};
+
+	converter<sch::T> conv{
+		{0.1},
+		{0.1}
+	};
+
+	print_lbm_meta<sch::T, sch::D3Q27>(conv, {1e-3});
+
+	//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};
+	//}
+
+	/**
+	 * Set meta information describing what this cell is
+	 */
+	set_geometry(lattice);
+
+	/**
+	 *
+	 */
+	set_initial_conditions(lattice);
+
+	/**
+	 * Timeloop
+	 */
+
+	/**
+	 * Print basic setup info
+	 */
+	apply_for_cells([](auto& cell, std::size_t i, std::size_t j){
+			// Not needed
+			(void) j;
+			std::cout<<(static_cast<uint32_t>(cell.template get<"info">()({0}).get()));
+			if( (i+1) < dim_x){
+				std::cout<<" ";
+			}else{
+				std::cout<<"\n";
+			}
+	}, lattice);
+
+	uint64_t lattice_steps = 512000u;
+	bool even_step = true;
+
+	uint64_t print_every = 256u;
+	uint64_t file_no = 0u;
+
+	saw::data<sch::Array<sch::MacroStruct<sch::T,sch::D3Q27::D>,sch::D3Q27::D>> macros{dim};
+
+	for(uint64_t step = 0; step < lattice_steps; ++step){
+
+		if(step % print_every == 0u){
+			std::cout<<"\n";
+			typename saw::native_data_type<sch::T>::type sum = 0.0;
+
+
+			apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){
+				auto& dfs = cell.template get<"dfs">();
+				typename saw::native_data_type<sch::T>::type rho;
+				std::array<typename saw::native_data_type<sch::T>::type, sch::D3Q27::D> vel;
+				compute_rho_u<sch::T,sch::D3Q27>(dfs,rho,vel);
+
+				if(i > 1 && (i+2) < dim_x && j > 1 && (j+2) < dim_y){
+					sum += rho;
+				}
+				auto angle = std::atan2(vel[1],vel[0]);
+				double vel_mag = vel[1] * vel[1] + vel[0] * vel[0];
+
+				auto dir_char = [&]() -> std::string_view {
+					/**
+					 * Flipped y due to print order
+					 */
+					constexpr auto pi = M_PI;
+					return "■";
+					if(vel_mag < 1e-4){
+						return "•";
+					}
+					if(angle > 7.0 / 8.0 * pi){
+						return "←";
+					}
+					if(angle > 5.0 / 8.0 * pi){
+						return "↙";
+					}
+					if(angle > 3.0 / 8.0 * pi){
+						return "↓";
+					}
+					if(angle > 1.0 / 8.0 * pi){
+						return "↘";
+					}
+					if(angle > -1.0 / 8.0 * pi){
+						return "→";
+					}
+					if(angle > -3.0 / 8.0 * pi){
+						return "↗";
+					}
+					if(angle > -5.0 / 8.0 * pi){
+						return "↑";
+					}
+					if(angle > -7.0 / 8.0 * pi){
+						return "↖";
+					}
+					return "←";
+				}();
+
+				std::array<uint32_t,3u> rgb_start{64,64,255};
+				std::array<uint32_t,3u> rgb_stop{255,64,64};
+				std::array<uint32_t,3u> rgb_middle{255,255,255};
+
+
+				auto col_interpol = [&](){
+					std::array<uint32_t, 3u> rgb_interpol = rgb_start;
+					double vel_mag_cut = std::min(1.0,std::max(vel_mag/(0.07*0.07),0.0));
+
+					if(vel_mag_cut < 0.5){
+						uint32_t vel_mag_i = static_cast<uint32_t>(2.0 * vel_mag_cut * 256);
+						for(uint8_t i = 0u; i < 3u; ++i){
+							rgb_interpol[i] = (rgb_middle[i] * vel_mag_i + (256-vel_mag_i) * rgb_start[i]) / 256;
+						}
+					}else{
+						uint32_t vel_mag_i = static_cast<uint32_t>((2.0*(vel_mag_cut-0.5)) * 256);
+						for(uint8_t i = 0u; i < 3u; ++i){
+							rgb_interpol[i] = (rgb_stop[i] * vel_mag_i + (256-vel_mag_i) * rgb_middle[i]) / 256;
+						}
+					}
+					return rgb_interpol;
+				};
+				auto rgb_interpol = col_interpol();
+
+				std::cout<<"\x1b[38;2;"<<rgb_interpol[0]<<";"<<rgb_interpol[1]<<";"<<rgb_interpol[2]<<"m"<<dir_char;
+				if( (i+1) < dim_x){
+					std::cout<<" ";
+				}else{
+					std::cout<<"\n";
+				}
+			}, lattice);
+			std::cout<<"\x1b[38;2;255;255;255m";
+
+			/// std::cout<<"Average rho: "<<(sum / ((dim_x-4)*(dim_y-4)))<<std::endl;
+			std::cout.flush();
+
+			apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){
+				auto dfs = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
+
+				auto& rho = macros.at({{i,j}}).template get<"pressure">();
+				auto& vel = macros.at({{i,j}}).template get<"velocity">();
+				compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel);
+			}, lattice);
+
+			std::string vtk_f_name{"tmp/cavity_3d_"};
+			vtk_f_name += std::to_string(file_no) + ".vtk";
+			write_vtk_file(vtk_f_name, macros);
+
+			++file_no;
+		}
+
+		lbm_step(lattice, even_step, step);
+
+		even_step = not even_step;
+	}
+
+	/**
+	 * Flush cout
+	 */
+	std::cout<<"\n\n";
+	std::cout.flush();
+	return 0;
+}
diff --git a/examples/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d.cpp
index 9498ba7..902a19e 100644
--- a/examples/poiseulle_particles_channel_2d.cpp
+++ b/examples/poiseulle_particles_channel_2d.cpp
@@ -451,9 +451,8 @@ void couple_particles_to_lattice(
 					// If neighbour is wall, then add force pushing the particle away
 					if(n_info.get() <= 1u){
 						// add to p_acc
-						p_acc.at({{0u}}) = p_acc.at({{0u}}) + saw::data<sch::Int32>{10 * 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>{10 * dfi::directions[dfi::opposite_index[k.get()]][1u]}.template cast_to<sch::T>();
-						;
+						p_acc.at({{0u}}) = p_acc.at({{0u}}) - saw::data<sch::Int32>{10 * 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>{10 * dfi::directions[dfi::opposite_index[k.get()]][1u]}.template cast_to<sch::T>();
 					}
 				}
 				/// 2. Add force pushing away from wall