Reworking from AoS to SoA for data handling

8 files changed, 109 insertions, 21 deletions

diff --git a/c++/converter.hpp b/c++/converter.hpp
index bc307f5..e07847f 100644
--- a/c++/converter.hpp
+++ b/c++/converter.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include "lbm_unit.h"
+#include "lbm_unit.hpp"
 
 namespace kel {
 namespace lbm {
diff --git a/c++/descriptor.hpp b/c++/descriptor.hpp
index 5ee9918..376e733 100644
--- a/c++/descriptor.hpp
+++ b/c++/descriptor.hpp
@@ -22,14 +22,14 @@ struct Cell {
 	static constexpr uint64_t Size = SC + Desc::D * DC + Desc::Q * QC;
 };
 
-template<typename Desc, typename CellStruct>
-struct Field;
+template<typename Desc, typename Cell>
+struct CellFields;
 
-template<typename Desc, typename... CellMembers>
-struct Field<
+template<typename Desc, typename... CellFieldTypes, saw::string_literal... CellFieldNames>
+struct CellFields<
 	Desc,
 	Struct<
-		CellMembers...
+		Member<Array<CellFieldTypes,Desc::D>, CellFieldNames>...
 	>
 >;
 
@@ -38,6 +38,37 @@ struct Field<
 template<typename T, typename Desc>
 class df_info{};
 
+template<typename T> 
+class df_info<T,sch::Descriptor<1,3>> {
+public:
+	using Descriptor = sch::Descriptor<1,3>;
+
+	static constexpr uint64_t D = 1u;
+	static constexpr uint64_t Q = 3u;
+
+	static constexpr std::array<std::array<int32_t,D>,Q> directions = {{
+		{ 0},
+		{-1},
+		{ 1}
+	}};
+
+	static constexpr std::array<typename saw::native_data_type<T>::type, Q> weights = {
+		2./3.,
+		1./6.,
+		1./6.
+	};
+
+	static constexpr std::array<uint64_t,Q> opposite_index = {
+		0,2,1
+	};
+	
+	static constexpr typename saw::native_data_type<T>::type inv_cs2 = 3.0;
+	static constexpr typename saw::native_data_type<T>::type cs2 = 1./3.;
+};
+
+/**
+ * D2Q5 Descriptor
+ */
 template<typename T>
 class df_info<T,sch::Descriptor<2, 5>> {
 public:
@@ -132,9 +163,9 @@ struct meta_schema<kel::lbm::sch::Cell<T,Desc,S,D,Q>> {
 };
 
 template<typename Desc, typename CellT>
-struct meta_schema<kel::lbm::sch::Field<Desc, CellT>> {
+struct meta_schema<kel::lbm::sch::CellFields<Desc, CellT>> {
   using MetaSchema = schema::FixedArray<schema::UInt64,Desc::D>;
-	using Schema = kel::lbm::sch::Field<Desc, CellT>;
+	using Schema = kel::lbm::sch::CellFields<Desc, CellT>;
 };
 
 template<typename Sch, typename Desc, uint64_t S, uint64_t D, uint64_t Q, typename Encode>
@@ -156,16 +187,18 @@ public:
 };
 
 template<typename Desc, typename CellT, typename Encode>
-class data<kel::lbm::sch::Field<Desc, CellT>, Encode> final {
+class data<kel::lbm::sch::CellFields<Desc, CellT>, Encode> final {
 public:
-	using Schema = kel::lbm::sch::Field<Desc,CellT>;
+	using Schema = kel::lbm::sch::CellFields<Desc,CellT>;
 	using MetaSchema = typename meta_schema<Schema>::MetaSchema;
 private:
-	data<schema::Array<CellT,Desc::D>, Encode> inner_;
+	data<CellT, Encode> inner_;
+	data<MetaSchema, Encode> inner_meta_;
 public:
 	data() = default;
-	data(const data<schema::FixedArray<schema::UInt64, Desc::D>>& inner_meta__):
-		inner_{inner_meta__}
+	data(const data<MetaSchema,Encode>& inner_meta__):
+		inner_{inner_meta__},
+		inner_meta_{inner_meta__}
 	{}
 
 	template<uint64_t i>
diff --git a/c++/equilibrium.hpp b/c++/equilibrium.hpp
index 8342ed4..326fc9e 100644
--- a/c++/equilibrium.hpp
+++ b/c++/equilibrium.hpp
@@ -9,7 +9,11 @@ saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<sch::T> rho,
 	using dfi = df_info<T, Descriptor>;
 
 	saw::data<sch::FixedArray<T,Descriptor::Q>> eq;
-	saw::data<T> vel_vel{};
+	//   ^
+	//  0.0
+	//  / \
+	//  | |
+	saw::data<T> vel_vel{0.0};
 	for(uint64_t j = 0u; j < Descriptor::D; ++j){
 		vel_vel = vel_vel + vel.at(j) * vel.at(j);
 	}
diff --git a/c++/examples/cavity_2d.cpp b/c++/examples/cavity_2d.cpp
index 9e91c6c..f93bc64 100644
--- a/c++/examples/cavity_2d.cpp
+++ b/c++/examples/cavity_2d.cpp
@@ -37,12 +37,12 @@ using CellInfo = Cell<UInt8, D2Q9, 1, 0, 0>;
  */
 template<typename Desc>
 using CellStruct = Struct<
-	Member<DfCell<Desc>, "dfs">,
-	Member<CellInfo<Desc>, "info">
+	Member<Array<DfCell<Desc>,Desc::D>, "dfs">,
+	Member<Array<CellInfo<Desc>,Desc::D>, "info">
 >;
 
 
-using CavityFieldD2Q9 = Field<D2Q9, CellStruct<D2Q9>>;
+using CavityFieldD2Q9 = CellFields<D2Q9, CellStruct<D2Q9>>;
 }
 
 /**
@@ -350,9 +350,8 @@ void lbm_step(
 	// Stream
 	for(uint64_t i = 1; (i+1) < old_latt.template get_dim_size<0>().get(); ++i){
 		for(uint64_t j = 1; (j+1) < old_latt.template get_dim_size<1>().get(); ++j){
-			auto& cell_new = new_latt({{i,j}});
-			auto& df_new = cell_new.template get<"dfs">();
-			auto& info_new = cell_new.template get<"info">();
+			auto& df_new = new_latt.template get<"dfs">()({{i,j}});
+			auto& info_new = new_latt.template get<"info">()({{i,j}});
 
 			if(info_new({0u}).get() > 0u && info_new({0u}).get() != 2u){
 				for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){
diff --git a/c++/lbm_unit.hpp b/c++/lbm_unit.hpp
index c4d09f6..dc1f8e9 100644
--- a/c++/lbm_unit.hpp
+++ b/c++/lbm_unit.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include <kel/unit/unit.h>
+#include <forstio/codec/unit/unit.hpp>
 
 #include <string_view>
 
diff --git a/c++/particle/particle.hpp b/c++/particle/particle.hpp
index fd0ff59..58c028c 100644
--- a/c++/particle/particle.hpp
+++ b/c++/particle/particle.hpp
@@ -20,6 +20,23 @@ template<typename T, uint64_t D>
 using ParticleMask = Struct<
 	Member<Array<T,D>, "mask">
 >;
+
+template<typename T, uint64_t D>
+using Particle = Struct<
+	Member<ParticleRigidBody<T,D>, "rigid_body">
+>;
 }
+
+template<typename T, uint64_t D>
+class particle_system {
+private:
+	saw::data<sch::Array<sch::Particle<T,D>>> particles_;
+public:
+
+	void step(T time_step){
+		for(auto& iter : particles_){
+		}
+	}
+};
 }
 }
diff --git a/tests/descriptor.cpp b/tests/descriptor.cpp
index 5e23ebc..a8337e6 100644
--- a/tests/descriptor.cpp
+++ b/tests/descriptor.cpp
@@ -19,6 +19,11 @@ void check_opposite_dirs(){
 	}
 }
 
+SAW_TEST("Opposites and Dirs D1Q3"){
+	using namespace kel;
+	check_opposite_dirs<lbm::sch::Descriptor<1,3>>();
+}
+
 SAW_TEST("Opposites and Dirs D2Q5"){
 	using namespace kel;
 	check_opposite_dirs<lbm::sch::Descriptor<2,5>>();
diff --git a/tests/equilibrium.cpp b/tests/equilibrium.cpp
new file mode 100644
index 0000000..88542d2
--- /dev/null
+++ b/tests/equilibrium.cpp
@@ -0,0 +1,30 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/equilibrium.hpp"
+
+
+namespace {
+
+template<typename Descriptor>
+void check_equilibrium(){
+	using namespace kel;
+
+	using dfi = df_info<lbm::sch::Float64,Descriptor>;
+
+	saw::data<lbm::sch::Float64> rho{1.0};
+	saw::data<lbm::sch::FixedArray<lbm::sch::Float64,>> vel;
+	for(saw::data<lbm::sch::UInt64> i{0u}; i < {Descriptor::D}; ++i){
+		vel.at(i) = {0.0};
+	}
+	auto eq = lbm::equilibrium<lbm::sch::Float64,Descriptor>(rho,vel);
+
+	for(saw::data<lbm::sch::UInt64> i{0u}; i < {Descriptor::Q}; ++i){
+		SAW_CHECK(eq(i) == dfi::weights[i.get()], "No velocity and normalized rho should be exactly the weights");
+	}
+}
+
+SAW_TEST("Equilibrium At rest D1Q3"){
+	using namespace kel;
+	check_equilibrium<lbm::sch::Descriptor<1,3>>();
+}
+}