Merge branch 'dev'HEADmaster

11 files changed, 142 insertions, 99 deletions

diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
index a41d2b7..b2061f9 100644
--- a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
+++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp
@@ -155,7 +155,14 @@ saw::error_or<void> setup_initial_conditions(
 		auto& rbp = rb.template get<"position">();
 		rbp.at({{0u}}) = 0.25 * dim_x;
 		rbp.at({{1u}}) = 0.5 * dim_y;
-		auto& rbp_old = rbp;
+		rb.template get<"position_old">() = rbp;
+		saw::data<sch::Vector<T,Desc::D>> zero;
+		zero.at({{0u}}) = 0.0f;
+		zero.at({{1u}}) = 0.0f;
+		rb.template get<"acceleration">() = zero;
+		rb.template get<"rotation">() = {};
+		rb.template get<"rotation_old">() = {};
+		rb.template get<"angular_acceleration">() = {};
 	}
 
 	return saw::make_void();
diff --git a/modules/core/c++/descriptor/common.hpp b/modules/core/c++/descriptor/common.hpp
new file mode 100644
index 0000000..f821510
--- /dev/null
+++ b/modules/core/c++/descriptor/common.hpp
@@ -0,0 +1,9 @@
+#pragma once
+
+#include "../common.hpp"
+
+namespace kel {
+namespace lbm {
+
+}
+}
diff --git a/modules/core/c++/descriptor/d1q3.hpp b/modules/core/c++/descriptor/d1q3.hpp
new file mode 100644
index 0000000..7fc8591
--- /dev/null
+++ b/modules/core/c++/descriptor/d1q3.hpp
@@ -0,0 +1,9 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+
+}
+}
diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp
index 6fce4ca..799d2b5 100644
--- a/modules/core/c++/hlbm.hpp
+++ b/modules/core/c++/hlbm.hpp
@@ -79,6 +79,7 @@ public:
 
 		auto& N = particle_N_f.at(index);
 		auto& D = particle_D_f.at(index);
+
 		// Convex combination of velocities
 		vel = vel * porosity + [&]() -> saw::data<sch::Vector<T,Desc::D>> {
 			return (D.at({}).get() > 0.0 ? N * flip_porosity / D : N);
@@ -113,13 +114,20 @@ public:
 		using dfi = df_info<T,Desc>;
 		bool is_even = ((time_step.get() % 2) == 0);
 		
+		saw::data<sch::Scalar<T>> one;
+		one.at({}) = 1.0;
+		
 		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
 
 		auto& part_spheroid_group = part_group;
-		auto& mvel = macros.template get<"velocity">();
+		auto& mvel_f = macros.template get<"velocity">();
 		auto& mpor_f = macros.template get<"porosity">();
+		
+		auto& particle_N_f = field.template get<"particle_N">();
+		auto& particle_D_f = field.template get<"particle_D">();
+
 		{
-			auto& parts = part_spheroid_group.template get<"particles">();
+			auto parts = part_spheroid_group.template get<"particles">();
 			auto parts_size = parts.meta().at({0u});
 
 			auto& p_coll = part_spheroid_group.template get<"collision">().at({});
@@ -132,7 +140,7 @@ public:
 			saw::data<sch::FixedArray<sch::UInt64,Desc::D>> start;
 			saw::data<sch::FixedArray<sch::UInt64,Desc::D>> stop;
 
-			auto eo_aabb = particle_aabb<typename ParticleSchema::ValueType>::calculate(part_spheroid_group,{{0u}},mvel.meta());
+			auto eo_aabb = particle_aabb<typename ParticleSchema::ValueType>::calculate(part_spheroid_group,{{0u}},mvel_f.meta());
 			if(eo_aabb.is_error()){
 				return;
 			}
@@ -143,21 +151,22 @@ public:
 			stop = aabb.template get<"b">();
 
 			saw::data<sch::Vector<T,Desc::D>> force_p;
-
-			iterator<Desc::D>::apply([&](const auto& index){
+			
+			iterator<Desc::D>::apply([&](const auto& index_f){
 				// ask for the d_k value here.
 				// For every value im iterating over I need sth
 				// std::cout<<"Pos: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl;
 			
-				auto& dfs = dfs_old_f.at(index);
+				auto& dfs = dfs_old_f.at(index_f);
 				
 				saw::data<sch::Vector<T,Desc::D>> momentum;
+
 				for(uint64_t i = 0u; i < Desc::Q; ++i){
 					saw::data<sch::Vector<T,Desc::D>> e_i;
 					saw::data<sch::FixedArray<sch::UInt64,Desc::D>> n_ind_i;
 					for(uint64_t k{0u}; k < Desc::D; ++k){
 						e_i.at({{k}}) = (dfi::directions[i])[k];
-						n_ind_i.at({k}) = (dfi::directions[i])[k];
+						n_ind_i.at({k}) = index_f.at({k}) + (dfi::directions[i])[k];
 					}
 
 					uint64_t i_opp = dfi::opposite_index[i];
@@ -170,21 +179,26 @@ public:
 					momentum = momentum + ei_dfs;
 				}
 
-				auto& mpor = mpor_f.at(index);
-				auto rel_dist = saw::math::vectorize_data(index).template cast_to<T>() - pirb_pos;
+				auto& mpor = mpor_f.at(index_f);
+				auto rel_dist = saw::math::vectorize_data(index_f).template cast_to<T>() - pirb_pos;
 				saw::data<sch::Scalar<T>> eps;
 				eps.at({}) = 1.5f;
 				mpor = particle_porosity<T,Desc::D,1u,por::ParticleSpheroid<T>>::calculate(rel_dist,p_rad,eps);
 				force_p = force_p + momentum * mpor;
 
+				auto& N = particle_N_f.at(index_f);
+				auto& D = particle_D_f.at(index_f);
+				auto flip_porosity = one - mpor;
+				D = D + flip_porosity;
+				N = N + mvel_f.at(index_f) * flip_porosity;
 			},start,stop);
-
+			
 			auto& pirb_acc = pirb.template get<"acceleration">();
 			pirb_acc = force_p;
 
 			saw::data<sch::Scalar<T>> ts;
 			ts.at({}) = 1u;
-			verlet_step_lambda<T,Desc::D,1u,coll::Spheroid<T>>(part_spheroid_group,{{{0u}}},ts);
+			verlet_step_lambda<T,Desc::D>(pi,ts);
 
 			// Check
 		}
diff --git a/modules/core/c++/particle/particle.hpp b/modules/core/c++/particle/particle.hpp
index 29b08e6..f759c8d 100644
--- a/modules/core/c++/particle/particle.hpp
+++ b/modules/core/c++/particle/particle.hpp
@@ -137,7 +137,7 @@ saw::data<sch::Particle<T,D, sch::ParticleCollisionSpheroid<T>>> create_spheroid
 */
 
 template<typename T,uint64_t D>
-constexpr auto verlet_step_lambda_old = [](saw::data<sch::Particle<T,D>>& particle, saw::data<sch::Scalar<T>> time_step_delta){
+constexpr auto verlet_step_lambda = [](saw::data<sch::Particle<T,D>>& particle, saw::data<sch::Scalar<T>> time_step_delta){
 	auto& body = particle.template get<"rigid_body">();
 
 	auto& pos = body.template get<"position">();
@@ -158,51 +158,14 @@ constexpr auto verlet_step_lambda_old = [](saw::data<sch::Particle<T,D>>& partic
 	two.at({}).set(2.0);
 	pos_new = pos * two - pos_old + pos_acc * tsd_squared;
 
-	// Angular
-	saw::data<typename sch::impl::rotation_type_helper<T,D>::Schema> rot_new;
-	rot_new = rot * two - rot_old + rot_acc * tsd_squared;
-
-	// Swap - Could be std::swap?
+	// Swap - Could be std::swap? No, this is device code
 	pos_old = pos;
 	pos = pos_new;
 
-	rot_old = rot;
-	rot = rot_new;
-};
-
-template<typename T,uint64_t D,uint64_t PC, typename Coll>
-constexpr auto verlet_step_lambda = [](saw::data<sch::ParticleGroup<T,D,PC,Coll>>& pgrp, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::Scalar<T>> time_step_delta){
-	auto& parts = pgrp.template get<"particles">();
-	auto& particle = parts.at(index);
-
-	auto& body = particle.template get<"rigid_body">();
-
-	auto& pos = body.template get<"position">();
-	auto& pos_old = body.template get<"position_old">();
-
-	auto& pos_acc = body.template get<"acceleration">();
-	
-	auto& rot = body.template get<"rotation">();
-	auto& rot_old = body.template get<"rotation_old">();
-
-	auto& rot_acc = body.template get<"angular_acceleration">();
-
-	auto tsd_squared = time_step_delta * time_step_delta;
-
-	saw::data<sch::Vector<T,D>> pos_new;
-	// Actual step
-	saw::data<sch::Scalar<T>> two;
-	two.at({}).set(2.0);
-	pos_new = pos * two - pos_old + pos_acc * tsd_squared;
-
 	// Angular
 	saw::data<typename sch::impl::rotation_type_helper<T,D>::Schema> rot_new;
 	rot_new = rot * two - rot_old + rot_acc * tsd_squared;
-
-	// Swap - Could be std::swap?
-	pos_old = pos;
-	pos = pos_new;
-
+	
 	rot_old = rot;
 	rot = rot_new;
 };
diff --git a/modules/core/tests/particles.cpp b/modules/core/tests/particles.cpp
index 592c367..e060097 100644
--- a/modules/core/tests/particles.cpp
+++ b/modules/core/tests/particles.cpp
@@ -24,7 +24,7 @@ SAW_TEST("Verlet step 2D - Planar"){
 
 	saw::data<sch::Scalar<sch::T>> dt;
 	dt.at({}).set(0.5);
-	lbm::verlet_step_lambda_old<sch::T,2u>(particle,dt);
+	lbm::verlet_step_lambda<sch::T,2u>(particle,dt);
 
 	SAW_EXPECT(pos.at({{0}}).get() == 0.25, std::string{"Incorrect Pos X: "} + std::to_string(pos.at({{0}}).get()));
 	SAW_EXPECT(pos.at({{1}}).get() == 0.0, std::string{"Incorrect Pos Y: "} + std::to_string(pos.at({{1}}).get()));
diff --git a/modules/core/tests/schema.cpp b/modules/core/tests/schema.cpp
new file mode 100644
index 0000000..b705af0
--- /dev/null
+++ b/modules/core/tests/schema.cpp
@@ -0,0 +1,19 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/abstract/schema.hpp"
+
+namespace {
+namespace ts {
+using namespace kel::lbm::sch;
+
+using TPtr = Ptr<UInt16>;
+using TVal = UInt16;
+}
+
+SAW_TEST("Abstract/Schema/Ptr"){
+	using namespace kel::lbm;
+
+	// static_assert(std::is_same_v<typename schema_decay<ts::TPtr>::Schema, ts::TVal>, "schema decay failed");
+}
+
+}
diff --git a/poc/triangle_stream/common.hpp b/poc/triangle_stream/common.hpp
new file mode 100644
index 0000000..927f8f8
--- /dev/null
+++ b/poc/triangle_stream/common.hpp
@@ -0,0 +1,46 @@
+#pragma once
+
+#include <kel/lbm/lbm.hpp>
+
+namespace kel {
+namespace lbm {
+
+constexpr uint64_t dim_y = 256ul;
+constexpr uint64_t dim_x = dim_y * 20ul;
+
+namespace sch {
+using namespace saw::schema;
+
+using InfoChunk = Chunk<UInt8, 0u, dim_x, dim_y>;
+
+template<typename T, typename Desc>
+using DfChunk = Chunk<FixedArray<T,Desc::Q>, 1u, dim_x, dim_y>;
+
+template<typename T, typename Desc>
+using ScalarChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
+
+template<typename T, typename Desc>
+using VectorChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>;
+
+template<typename T, typename Desc>
+using ChunkStruct = Struct<
+	Member<InfoChunk, "info">,
+	Member<DfChunk<T,Desc>, "dfs">,
+	Member<DfChunk<T,Desc>, "dfs_old">
+>;
+
+template<typename T, typename Desc>
+using VelChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>;
+
+template<typename T>
+using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
+
+template<typename T, typename Desc>
+using MacroStruct = Struct<
+	Member<VelChunk<T,Desc>, "velocity">,
+	Member<RhoChunk<T>, "density">
+>;
+}
+
+}
+}
diff --git a/poc/triangle_stream/init.hpp b/poc/triangle_stream/init.hpp
new file mode 100644
index 0000000..532f08c
--- /dev/null
+++ b/poc/triangle_stream/init.hpp
@@ -0,0 +1,13 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+
+template<typename T, typename Desc>
+void init(){
+}
+
+}
+}
diff --git a/poc/triangle_stream/sim.cpp b/poc/triangle_stream/sim.cpp
index 36986e4..b2121ee 100644
--- a/poc/triangle_stream/sim.cpp
+++ b/poc/triangle_stream/sim.cpp
@@ -1,53 +1,8 @@
-#include <kel/lbm/lbm.hpp>
+#include "common.hpp"
 
 namespace kel {
 namespace lbm {
 
-constexpr uint64_t dim_y = 256ul;
-constexpr uint64_t dim_x = dim_y * 20ul;
-
-namespace sch {
-using namespace saw::schema;
-
-using InfoChunk = Chunk<UInt8, 0u, dim_x, dim_y>;
-
-template<typename T, typename Desc>
-using DfChunk = Chunk<FixedArray<T,Desc::Q>, 1u, dim_x, dim_y>;
-
-template<typename T, typename Desc>
-using ScalarChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
-
-template<typename T, typename Desc>
-using VectorChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>;
-
-template<typename T, typename Desc>
-using ChunkStruct = Struct<
-	Member<InfoChunk, "info">,
-	Member<DfChunk<T,Desc>, "dfs">,
-	Member<DfChunk<T,Desc>, "dfs_old">
->;
-
-template<typename T, typename Desc>
-using VelChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>;
-
-template<typename T>
-using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
-
-template<typename T, typename Desc>
-using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
-	Member<RhoChunk<T>, "density">
->;
-
-//template<typename T, typename Desc>
-//using ParticleArray = Array<
-//	Particle<T,Desc::D>
-//>;
-}
-
-template<typename T, typename Desc>
-void init(){
-}
 
 template<typename T, typename Desc>
 saw::error_or<void> kmain(int argc, char** argv){
diff --git a/poc/triangle_stream/step.hpp b/poc/triangle_stream/step.hpp
new file mode 100644
index 0000000..e10c4db
--- /dev/null
+++ b/poc/triangle_stream/step.hpp
@@ -0,0 +1,8 @@
+#pragma once
+
+#include "schema.hpp"
+
+namespace kel {
+namespace lbm {
+}
+}