Merge branch 'dev'HEADmaster

7 files changed, 30 insertions, 46 deletions

diff --git a/examples/one_particle_sedimentation_2d/common.hpp b/examples/one_particle_sedimentation_2d/common.hpp
index 954557d..18d6584 100644
--- a/examples/one_particle_sedimentation_2d/common.hpp
+++ b/examples/one_particle_sedimentation_2d/common.hpp
@@ -62,7 +62,8 @@ using MacroStruct = Struct<
 	Member<VelChunk<T,Desc>, "momentum">,
 	Member<RhoChunk<T>, "density">,
 	Member<ScalarChunk<T,Desc>, "porosity">,
-	Member<VelChunk<T,Desc>, "force">
+	Member<VelChunk<T,Desc>, "force">,
+	Member<VelChunk<T,Desc>, "external_force">
 >;
 
 template<typename T, typename Desc>
diff --git a/examples/one_particle_sedimentation_2d/sim.cpp b/examples/one_particle_sedimentation_2d/sim.cpp
index b0b374c..a722fc3 100644
--- a/examples/one_particle_sedimentation_2d/sim.cpp
+++ b/examples/one_particle_sedimentation_2d/sim.cpp
@@ -48,10 +48,10 @@ saw::error_or<void> lbm_main(const saw::data<args::LbmArgs>& args){
 		// delta_x
 		{{0.02 / dim_x}},
 		// delta_t
-		{{1.067e-4}}
+		{{1.5e-6}}
 	};
 	
-	print_lbm_meta<T,Desc>(conv,{1e-3},{1e-4},{0.5 * dim_y});
+	print_lbm_meta<T,Desc>(conv,{1e-3},{1e-4},{0.02});
 
 	// saw::data<sch::FixedArray<sch::UInt64,Desc::D>> meta{{dim_x,dim_y}};
 	auto lbm_data_ptr = saw::heap<saw::data<sch::ChunkStruct<T,Desc>>>();
diff --git a/examples/one_particle_sedimentation_2d/step.hpp b/examples/one_particle_sedimentation_2d/step.hpp
index 20b7f4b..86064ae 100644
--- a/examples/one_particle_sedimentation_2d/step.hpp
+++ b/examples/one_particle_sedimentation_2d/step.hpp
@@ -19,8 +19,7 @@ saw::error_or<void> step(
 	auto& q = dev.get_handle();
 	auto& info_f = fields.template get<"info">();
 	auto& porous_f = macros.template get<"porosity">();
-
-
+	auto& ext_force_f = macros.template get<"external_force">();
 	
 	if constexpr ( std::is_same_v<Coll, method::Hlbm> ){
 		q.submit([&](acpp::sycl::handler& h){
@@ -31,6 +30,7 @@ saw::error_or<void> step(
 					index.at({{i}}).set(idx[i]);
 				}
 				hlbm_reset.apply(fields,macros,index,t_i);
+
 			});
 		}).wait();
 		q.submit([&](acpp::sycl::handler& h){
@@ -49,7 +49,6 @@ saw::error_or<void> step(
 		q.submit([&](acpp::sycl::handler& h){
 			component<T,Desc,cmpt::Hlbm,encode::Sycl<saw::encode::Native>> collision{1.0};
 			component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb;
-			component<T,Desc,cmpt::ZouHeHorizontal<false>,encode::Sycl<saw::encode::Native>> flow_out{1.0};
 
 			h.parallel_for(acpp::sycl::range<Desc::D>{dim_x,dim_y}, [=](acpp::sycl::id<Desc::D> idx){
 				saw::data<sch::FixedArray<sch::UInt64,Desc::D>> index;
@@ -67,27 +66,6 @@ saw::error_or<void> step(
 					break;
 					case 2u:
 						collision.apply(fields,macros,index,t_i);
-						break;
-					case 3u:
-						{
-							component<T,Desc,cmpt::ZouHeVelocityX<true>,encode::Sycl<saw::encode::Native>> flow_in{
-								[&]() -> saw::data<sch::Vector<T,Desc::D>> {
-									saw::data<sch::Vector<T,Desc::D>> vel;
-									{
-										auto y_si = conv.meter_lbm_to_si(index.at({{1u}}).template cast_to<T>()).handle();
-										vel.at({{0u}}) = conv.velocity_si_to_lbm({{static_cast<typename saw::native_data_type<T>::type>((1.2 / 0.41) * y_si.get() - (1.2 / (0.41*0.41)) * y_si.get() * y_si.get())}}).handle();
-										vel.at({{1u}}) = 0.0f;
-									}
-									return vel;
-								}()
-							};
-							flow_in.apply(fields,index,t_i);
-						}
-						collision.apply(fields,macros,index,t_i);
-						break;
-					case 4u:
-						flow_out.apply(fields,index,t_i);
-						collision.apply(fields,macros,index,t_i);
 					break;
 					default:
 					break;
@@ -104,11 +82,16 @@ saw::error_or<void> step(
 				}
 				
 				fplbm_reset.apply(fields,macros,index,t_i);
+
+				// TODO apply external force
+				// auto& ext_force = ext_force_f.at(index);
 			});
 		}).wait();
 
 		q.submit([&](acpp::sycl::handler& h){
-			component<T,Desc,cmpt::FpLbmOneParticle,encode::Sycl<saw::encode::Native>> fplbm_one_part;
+			saw::data<sch::Vector<T,Desc::D>> gf;
+			gf.at({{1u}}) = conv.acceleration_si_to_lbm({-9.81}).handle();
+			component<T,Desc,cmpt::FpLbmOneParticle,encode::Sycl<saw::encode::Native>> fplbm_one_part{gf};
 
 			h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){
 				saw::data<sch::FixedArray<sch::UInt64,1u>> index;
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/common.hpp b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp
index be9a538..a85cd50 100644
--- a/examples/schaefer_turek_durst_krause_rannbacher/common.hpp
+++ b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp
@@ -21,8 +21,8 @@ namespace lbm {
  * position at (0.16m,0.2m)
  */
 
-constexpr uint64_t dim_y = 512ul;
-constexpr uint64_t dim_x = static_cast<uint64_t>(dim_y * (220.0 / 41.0));
+constexpr uint64_t dim_y = 41u;
+constexpr uint64_t dim_x = 220u;
 
 constexpr uint64_t particle_amount = 1ul;
 
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp
index 5a92681..a152789 100644
--- a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp
+++ b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp
@@ -87,12 +87,6 @@ saw::error_or<void> lbm_main(const saw::data<args::LbmArgs>& args){
 			return eov;
 		}
 	}
-	{
-		auto eov = write_vtk_file(out_dir,"initial_state",0u,*lbm_data_ptr);
-		if(eov.is_error()){
-			return eov;
-		}
-	}
 
 	saw::data<sch::ChunkStruct<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_data{sycl_q};
 	saw::data<sch::MacroStruct<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_macro_data{sycl_q};
diff --git a/examples/stdkr_2d_bgk/common.hpp b/examples/stdkr_2d_bgk/common.hpp
index eff9d84..0b597c0 100644
--- a/examples/stdkr_2d_bgk/common.hpp
+++ b/examples/stdkr_2d_bgk/common.hpp
@@ -24,8 +24,8 @@ struct Bgk {};
  * position at (0.16m,0.2m)
  */
 
-constexpr uint64_t dim_y = 512ul;
-constexpr uint64_t dim_x = static_cast<uint64_t>(dim_y * (220.0 / 41.0));
+constexpr uint64_t dim_y = 41u;
+constexpr uint64_t dim_x = 220u;
 
 constexpr uint64_t particle_amount = 1ul;
 
diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp
index 714a300..5ad692b 100644
--- a/modules/core/c++/fplbm.hpp
+++ b/modules/core/c++/fplbm.hpp
@@ -136,8 +136,14 @@ class component<T, Descriptor, cmpt::FpLbmOneParticle, Encode> final {
 public:
 	using Component = cmpt::FpLbmOneParticle;
 private:
+	saw::data<sch::Vector<T,Descriptor::D>> global_force_;
 public:
-	component() = default;
+	component():
+		global_force_{}
+	{}
+	component(saw::data<sch::Vector<T,Descriptor::D>> global_force__):
+		global_force_{global_force__}
+	{}
 
 	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 {
@@ -148,7 +154,7 @@ public:
 		// auto& dfs = dfs_old_f.at(index);
 
 		auto& rho_f = macros.template get<"density">();
-		auto& vel_f = macros.template get<"momentum">();
+		auto& mom_f = macros.template get<"momentum">();
 		auto& por_f = macros.template get<"porosity">();
 		auto& force_f = macros.template get<"force">();
 
@@ -170,7 +176,7 @@ public:
 		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());
+		auto eo_aabb = particle_aabb<typename ParticleSchema::ValueType>::calculate(pg,{{0u}},mom_f.meta());
 		if(eo_aabb.is_error()){
 			return;
 		}
@@ -183,7 +189,7 @@ public:
 		two.at({}).set(2.0f);
 
 		saw::data<sch::Scalar<T>> one;
-		one.at({}) = 1.0;
+		one.at({}) = 1.0f;
 		
 		saw::data<sch::Scalar<T>> eps;
 		eps.at({}) = 1.5f;
@@ -192,11 +198,11 @@ public:
 		sss.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{};
+		saw::data<sch::Vector<T,Descriptor::D>> force_p = global_force_;
 
 		iterator<Descriptor::D>::apply([&](const auto& index_f) -> void {
 			auto& force = force_f.at(index_f);
-			auto& vel = vel_f.at(index_f);
+			auto& mom = mom_f.at(index_f);
 			auto& por = por_f.at(index_f);
 			auto& rho = rho_f.at(index_f);
 
@@ -207,8 +213,8 @@ public:
 
 			// vel_s is technically time the density of the particle?
 
-			force = ( vel_s * rho - vel * rho ) * flip_por;
-			// force = ( vel_s * rho - vel * rho ) /* (two / three) */ * flip_por;
+			force = ( vel_s - mom) * flip_por;
+			// force = (vel_s - mom) * flip_por / (flip_por / two + one);
 
 			force_p = force_p - force;
 		}, aabb.template get<"a">(), aabb.template get<"b">());