Danglingdev

5 files changed, 30 insertions, 156 deletions

diff --git a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp
index 328f6e2..cc7649a 100644
--- a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp
+++ b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp
@@ -6,6 +6,7 @@
 #include <forstio/remote/filesystem/easy.hpp>
 #include <forstio/codec/json/json.hpp>
 #include <forstio/codec/simple.hpp>
+#include <forstio/codec/args.hpp>
 
 namespace kel {
 namespace lbm {
@@ -18,6 +19,14 @@ constexpr uint64_t particle_amount = 1ul;
 namespace sch {
 using namespace saw::schema;
 
+using LbmArgs = Args<
+	Struct<
+		Member<String, "name">,
+		Member<UInt64, "sub_steps">
+	>,
+	Tuple<>
+>;
+
 using InfoChunk = Chunk<UInt8, 0u, dim_x, dim_y>;
 
 template<typename T, typename Desc>
@@ -185,7 +194,7 @@ saw::error_or<void> step(
 	}).wait();
 
 	q.submit([&](acpp::sycl::handler& h){
-		component<T,Desc,cmpt::FpLbmOneParticleImplicit,encode::Sycl<saw::encode::Native>> fplbm_one_part;
+		component<T,Desc,cmpt::FpLbmOneParticle,encode::Sycl<saw::encode::Native>> fplbm_one_part;
 
 		h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){
 			saw::data<sch::FixedArray<sch::UInt64,1u>> index;
@@ -193,7 +202,7 @@ saw::error_or<void> step(
 				index.at({{i}}).set(idx[i]);
 			}
 
-			fplbm_one_part.apply(fields,macros,particles,index,t_i,{16u});
+			fplbm_one_part.apply(fields,macros,particles,index,t_i,{64u});
 		});
 	}).wait();
 
diff --git a/examples/poiseulle_3d_gpu/sim.cpp b/examples/poiseulle_3d_gpu/sim.cpp
index 6916ef2..167234d 100644
--- a/examples/poiseulle_3d_gpu/sim.cpp
+++ b/examples/poiseulle_3d_gpu/sim.cpp
@@ -39,22 +39,17 @@ using ChunkStruct = Struct<
 >;
 
 template<typename T, typename Desc>
-using VelChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y, dim_z>;
+using MomChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y, dim_z>;
 
 template<typename T>
 using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y, dim_z>;
 
 template<typename T, typename Desc>
 using MacroStruct = Struct<
-	Member<VelChunk<T,Desc>, "velocity">,
+	Member<MomChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<ScalarChunk<T,Desc>, "porosity">
 >;
-
-//template<typename T, typename Desc>
-//using ParticleArray = Array<
-//	Particle<T,Desc::D>
-//>;
 }
 
 template<typename T, typename Desc>
@@ -106,7 +101,7 @@ saw::error_or<void> setup_initial_conditions(
 	//
 	auto& df_f = fields.template get<"dfs_old">();
 	auto& rho_f = macros.template get<"density">();
-	auto& vel_f = macros.template get<"velocity">();
+	auto& mom_f = macros.template get<"momentum">();
 	auto& por_f = macros.template get<"porosity">();
 	
 	iterator<Desc::D>::apply(
@@ -115,8 +110,8 @@ saw::error_or<void> setup_initial_conditions(
 			auto& rho = rho_f.at(index);
 			por_f.at(index).at({}) = {1};
 			rho.at({}) = {1};
-			auto& vel = vel_f.at(index);
-			auto eq = equilibrium<T,Desc>(rho,vel);
+			auto& mom = mom_f.at(index);
+			auto eq = equilibrium<T,Desc>(rho,mom);
 
 			df = eq;
 		},
@@ -129,9 +124,9 @@ saw::error_or<void> setup_initial_conditions(
 			auto& df = df_f.at(index);
 			auto& rho = rho_f.at(index);
 			rho.at({}) = {1};
-			auto& vel = vel_f.at(index);
-			vel.at({{0u}}) = 0.01;
-			auto eq = equilibrium<T,Desc>(rho,vel);
+			auto& mom = mom_f.at(index);
+			mom.at({{0u}}) = 0.01;
+			auto eq = equilibrium<T,Desc>(rho,mom);
 
 			df = eq;
 		},
@@ -165,13 +160,6 @@ saw::error_or<void> step(
 		// component<T,Desc,cmpt::HLBM,encode::Sycl<saw::encode::Native>> collision{0.6};
 		component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb;
 
-		saw::data<sch::Scalar<T>> rho_b;
-		rho_b.at({}) = 1.0;
-		saw::data<sch::Vector<T,Desc::D>> vel_b;
-		vel_b.at({{0u}}) = 0.01;
-
-		component<T,Desc,cmpt::Equilibrium,encode::Sycl<saw::encode::Native>> equi{rho_b,vel_b};
-
 		component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> flow_in{
 			[&](){
 				uint64_t target_t_i = 256u;
@@ -201,13 +189,11 @@ saw::error_or<void> step(
 					collision.apply(fields,macros,index,t_i);
 				break;
 				case 3u:
-					equi.apply(fields,index,t_i);
-					// flow_in.apply(fields,index,t_i);
+					flow_in.apply(fields,index,t_i);
 					collision.apply(fields,macros,index,t_i);
 				break;
 				case 4u:
-					equi.apply(fields,index,t_i);
-					// flow_out.apply(fields,index,t_i);
+					flow_out.apply(fields,index,t_i);
 					collision.apply(fields,macros,index,t_i);
 				break;
 				default:
diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp
index 1f38f51..37074b0 100644
--- a/modules/core/c++/fplbm.hpp
+++ b/modules/core/c++/fplbm.hpp
@@ -80,20 +80,25 @@ public:
 		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
 		auto& dfs = dfs_old_f.at(index);
 
+		saw::data<sch::Scalar<T>> half;
+		half.at({}) = 0.5f;
+
 		auto& rho_f = macros.template get<"density">();
 		saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
 
 		auto& vel_f = macros.template get<"velocity">();
-		saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);
+		saw::data<sch::Vector<T,Descriptor::D>> vel = vel_f.at(index);
+
+		auto& force_f = macros.template get<"force">();
+		auto& force = force_f.at(index);
+
+		vel = vel + force * half;
 
 		// compute_rho_u<T,Descriptor>(dfs,rho,vel);
 		auto eq = equilibrium<T,Descriptor>(rho,vel);
 
 		using dfi = df_info<T,Descriptor>;
 
-		auto& force_f = macros.template get<"force">();
-		auto& force = force_f.at(index);
-
 		auto& por_f = macros.template get<"porosity">();
 		auto& por = por_f.at(index);
 
@@ -214,93 +219,6 @@ public:
 };
 
 template<typename T, typename Descriptor, typename Encode>
-class component<T, Descriptor, cmpt::FpLbmOneParticleImplicit, Encode> final {
-public:
-	using Component = cmpt::FpLbmOneParticle;
-private:
-public:
-	component() = default;
-
-	template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema>
-	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& pg, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step, saw::data<sch::UInt64> sub_steps) const {
-	// void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, saw::data<sch::UInt64> time_step){
-		bool is_even = ((time_step.get() % 2) == 0);
-
-		//auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
-		// auto& dfs = dfs_old_f.at(index);
-
-		auto& rho_f = macros.template get<"density">();
-		auto& vel_f = macros.template get<"velocity">();
-		auto& por_f = macros.template get<"porosity">();
-		auto& force_f = macros.template get<"force">();
-
-		/**
-		 * The other methods all work with a flipped porosity.
-		 * For compat reasons this is also flipped
-		 */
-		// TODO - Change to tuple later
-		auto parts = pg.template get<"particles">();
-
-		auto parts_size = parts.meta().at({0u});
-		auto& pi = parts.at(index);		
-
-		auto& pirb = pi.template get<"rigid_body">();
-		auto& pirb_pos = pirb.template get<"position">();
-		auto& pirb_pos_old = pirb.template get<"position_old">();
-		
-
-		auto& p_coll = pg.template get<"collision">().at({});
-		auto& p_rad = p_coll.template get<"radius">();
-
-		auto eo_aabb = particle_aabb<typename ParticleSchema::ValueType>::calculate(pg,{{0u}},vel_f.meta());
-		if(eo_aabb.is_error()){
-			return;
-		}
-		auto& aabb = eo_aabb.get_value();
-
-		saw::data<sch::Scalar<T>> two;
-		two.at({}).set(2.0f);
-
-		saw::data<sch::Scalar<T>> one;
-		one.at({}) = 1.0;
-		
-		saw::data<sch::Scalar<T>> eps;
-		eps.at({}) = 1.5f;
-		
-		saw::data<sch::Scalar<T>> sss;
-		sss.at({}) = one.at({}) / sub_steps.template cast_to<T>();
-		auto vel_s = (pirb_pos-pirb_pos_old) / sss;
-
-		saw::data<sch::Vector<T,Descriptor::D>> force_p{};
-
-		iterator<Descriptor::D>::apply([&](const auto& index_f) -> void {
-			auto& force = force_f.at(index_f);
-			auto& vel = vel_f.at(index_f);
-			auto& por = por_f.at(index_f);
-			auto& rho = rho_f.at(index_f);
-
-			saw::data<sch::Vector<T,Descriptor::D>> rel_dist = saw::math::vectorize_data(index_f).template cast_to<T>() - pirb_pos;
-
-			por = particle_porosity<T,Descriptor::D,1u,por::ParticleSpheroid<T>>::calculate(rel_dist,p_rad,eps);
-			auto flip_por = one - por;
-
-			// vel_s is technically time the density of the particle?
-
-			force = ( vel_s * rho - vel * rho ) * two * flip_por / (one + flip_por);
-
-			force_p = force_p - force;
-		}, aabb.template get<"a">(), aabb.template get<"b">());
-
-		auto& pirb_acc = pirb.template get<"acceleration">();
-		pirb_acc = force_p / pg.template get<"total_mass">().at({});
-
-		for(saw::data<sch::UInt64> i{0u}; i < sub_steps; ++i){
-			verlet_step_lambda<T,Descriptor::D>(pi,sss);
-		}
-	}
-};
-
-template<typename T, typename Descriptor, typename Encode>
 class component<T, Descriptor, cmpt::FpLbmOneParticleNoVelocity, Encode> final {
 private:
 	saw::data<T> relaxation_;
diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp
index fc19728..cc21dcf 100644
--- a/modules/core/c++/hlbm.hpp
+++ b/modules/core/c++/hlbm.hpp
@@ -234,43 +234,5 @@ public:
 		}
 	}
 };
-
-template<typename T, typename Desc, typename Encode>
-class component<T, Desc, cmpt::HlbmOneParticleMomentumExchange, Encode> final {
-public:
-	template<typename CellFieldSchema, typename MacroFieldSchema>
-	void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Desc::D>> index, saw::data<sch::UInt64> time_step) const {
-		//
-		using dfi = df_info<T,Desc>;
-		bool is_even = ((time_step.get() % 2) == 0);
-		
-		auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
-		auto& dfs = dfs_old_f.at(index);
-
-		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];
-			}
-
-			uint64_t i_opp = dfi::opposite_index[i];
-
-			saw::data<T> dfs_added = dfs.at({i}) - dfs_old_f.at(n_ind_i).at({i_opp});
-			saw::data<sch::Scalar<T>> dfs_added_v;
-			dfs_added_v.at({}) = dfs_added;
-			auto ei_dfs = e_i * dfs_added_v;
-
-			momentum = momentum + ei_dfs;
-		}
-		auto& force_f = macros.template get<"force">();
-		
-		// Set Force
-		force_f.at(index) = momentum * macros.template get<"porosity">().at(index);
-	}
-};
 }
 }
diff --git a/modules/core/c++/lbm.hpp b/modules/core/c++/lbm.hpp
index 23d30b1..4e73537 100644
--- a/modules/core/c++/lbm.hpp
+++ b/modules/core/c++/lbm.hpp
@@ -1,6 +1,5 @@
 #pragma once
 
-#include "args.hpp"
 #include "schema.hpp"
 #include "descriptor.hpp"
 #include "boundary.hpp"