Experimental validation

6 files changed, 81 insertions, 72 deletions

diff --git a/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix b/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix
index d5361c2..958955f 100644
--- a/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix
+++ b/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix
@@ -12,7 +12,6 @@
 }:
 
 let
-	slip_txt = if slip then "slip" else "noslip";
 in stdenv.mkDerivation {
 	pname = "${pname}-examples-stdkr_2d_gpu_${particle_coupling}";
 	inherit version;
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/SConscript b/examples/schaefer_turek_durst_krause_rannbacher/SConscript
index 4483d58..b33b7ae 100644
--- a/examples/schaefer_turek_durst_krause_rannbacher/SConscript
+++ b/examples/schaefer_turek_durst_krause_rannbacher/SConscript
@@ -23,7 +23,7 @@ env.headers += examples_env.headers;
 # Cavity2D
 examples_objects = [];
 examples_env.add_source_files(examples_objects, ['sim.cpp'], shared=False);
-examples_env.poiseulle_2d_gpu = examples_env.Program('#bin/poiseulle_particles_2d_hlbm_gpu', [examples_objects]);
+examples_env.poiseulle_2d_gpu = examples_env.Program('#bin/stdkr_2d', [examples_objects]);
 
 # Set Alias
 env.examples = [
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/common.hpp b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp
index 2d51614..ec6f80f 100644
--- a/examples/schaefer_turek_durst_krause_rannbacher/common.hpp
+++ b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp
@@ -15,11 +15,14 @@ namespace lbm {
 
 /**
  * length 2.2m
+ * height 0.41m
  * viscosity 10^-3
+ * particle diameter 0.1m
+ * position at (0.16m,0.2m)
  */
 
-constexpr uint64_t dim_y = 256ul;
-constexpr uint64_t dim_x = dim_y * 20ul;
+constexpr uint64_t dim_y = 512ul;
+constexpr uint64_t dim_x = static_cast<uint64_t>(dim_y * (220.0 / 41.0));
 
 constexpr uint64_t particle_amount = 1ul;
 
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/init.hpp b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp
index b7fdeb0..b0e4020 100644
--- a/examples/schaefer_turek_durst_krause_rannbacher/init.hpp
+++ b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp
@@ -114,7 +114,7 @@ saw::error_or<void> init(
 
 	{
 		saw::data<sch::Scalar<T>> radius_p;
-		radius_p.at({}).set(2);
+		radius_p.at({}).set(conv.meter_si_to_lbm({{0.16}}).handle().get());
 		saw::data<sch::Scalar<T>> dense_p;
 		dense_p.at({}).set(1);
 		particles = create_spheroid_particle_group<T,Desc::D,1u>(radius_p, dense_p, {{16u}});
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp
index c7f1014..6113f6c 100644
--- a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp
+++ b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp
@@ -46,12 +46,12 @@ saw::error_or<void> lbm_main(const saw::data<args::LbmArgs>& args){
 
 	converter<T> conv {
 		// delta_x
-		{{1.0}},
+		{{0.41 / dim_y}},
 		// delta_t
-		{{1.0}}
+		{{1.067e-4}}
 	};
 	
-	print_lbm_meta<T,Desc>(conv,{0.1},{1e-4},{0.4 * dim_y});
+	print_lbm_meta<T,Desc>(conv,{1e-3},{1e-4},{0.5 * dim_y});
 
 	// 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>>>();
@@ -138,7 +138,7 @@ saw::error_or<void> lbm_main(const saw::data<args::LbmArgs>& args){
 				}
 			}
 			{
-				auto eov = write_csv_file(out_dir,"m",i.get(), *lbm_macro_data_ptr);
+				auto eov = write_csv_file(out_dir,"stdkr_2d",i.get(), *lbm_macro_data_ptr);
 				if(eov.is_error()){
 					return eov;
 				}
@@ -178,7 +178,7 @@ saw::error_or<void> lbm_main(const saw::data<args::LbmArgs>& args){
 		}
 	}
 	{
-		auto eov = write_vtk_file(out_dir,"m",time_steps.get(), *lbm_macro_data_ptr);
+		auto eov = write_vtk_file(out_dir,"stdkr_2d",time_steps.get(), *lbm_macro_data_ptr);
 		if(eov.is_error()){
 			return eov;
 		}
@@ -197,13 +197,17 @@ saw::error_or<void> k_main(int argc, char** argv){
 	saw::data<args::LbmArgs> args;
 	{
 		auto eov = saw::parse_args(args,argc,argv);
-		return eov;
+		if(eov.is_error()){
+			return eov;
+		}
 	}
 	auto& an = args.template get<"args">();
 	auto& coupling = an.template get<"coupling">();
 
 	if(coupling.stl_view() == "hlbm"){
 		return lbm_main<FloatT,DescT,method::Hlbm>(args);
+	} else if(coupling.stl_view() == "fplbm"){
+		return lbm_main<FloatT,DescT,method::FpLbm>(args);
 	}
 
 	return saw::make_error<saw::err::critical>("Invalid coupling");
diff --git a/examples/schaefer_turek_durst_krause_rannbacher/step.hpp b/examples/schaefer_turek_durst_krause_rannbacher/step.hpp
index 37eb4e7..1a70f6a 100644
--- a/examples/schaefer_turek_durst_krause_rannbacher/step.hpp
+++ b/examples/schaefer_turek_durst_krause_rannbacher/step.hpp
@@ -14,81 +14,84 @@ saw::error_or<void> step(
 		saw::data<sch::UInt64> t_i,
 		device& dev
 ){
+	static constexpr saw::data<sch::UInt64> sub_steps{0u};
+
 	auto& q = dev.get_handle();
 	auto& info_f = fields.template get<"info">();
 	auto& porous_f = macros.template get<"porosity">();
 	
-	q.submit([&](acpp::sycl::handler& h){
-		component<T,Desc,cmpt::HlbmReset,encode::Sycl<saw::encode::Native>> hlbm_reset;
-		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;
-			for(uint64_t i = 0u; i < Desc::D; ++i){
-				index.at({{i}}).set(idx[i]);
-			}
+	if constexpr ( std::is_same_v<Coll, method::Hlbm> ){
+		q.submit([&](acpp::sycl::handler& h){
+			component<T,Desc,cmpt::HlbmReset,encode::Sycl<saw::encode::Native>> hlbm_reset;
+			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;
+				for(uint64_t i = 0u; i < Desc::D; ++i){
+					index.at({{i}}).set(idx[i]);
+				}
 
-			hlbm_reset.apply(fields,macros,index,t_i);
-		});
-	}).wait();
+				hlbm_reset.apply(fields,macros,index,t_i);
+			});
+		}).wait();
+		q.submit([&](acpp::sycl::handler& h){
+			component<T,Desc,cmpt::HlbmOneParticle,encode::Sycl<saw::encode::Native>> hlbm_one_part;
 
-	q.submit([&](acpp::sycl::handler& h){
-		component<T,Desc,cmpt::HlbmOneParticle,encode::Sycl<saw::encode::Native>> hlbm_one_part;
+			h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){
+				saw::data<sch::FixedArray<sch::UInt64,1u>> index;
+				for(uint64_t i = 0u; i < 1u; ++i){
+					index.at({{i}}).set(idx[i]);
+				}
 
-		h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){
-			saw::data<sch::FixedArray<sch::UInt64,1u>> index;
-			for(uint64_t i = 0u; i < 1u; ++i){
-				index.at({{i}}).set(idx[i]);
-			}
+				hlbm_one_part.apply(fields,macros,particles,index,t_i,{0u});
+			});
+		}).wait();
 
-			hlbm_one_part.apply(fields,macros,particles,index,t_i,{16u});
-		});
-	}).wait();
+		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;
 
-	// auto coll_ev = 
-	q.submit([&](acpp::sycl::handler& h){
-		component<T,Desc,cmpt::Hlbm,encode::Sycl<saw::encode::Native>> collision{0.8};
-		component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb;
+			component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> flow_in{
+				[&](){
+					uint64_t target_t_i = 16u;
+					if(t_i.get() < target_t_i){
+						return 1.0 + (0.01 / target_t_i) * t_i.get();
+					}
+					return 1.001;
+				}()
+			};
+			component<T,Desc,cmpt::ZouHeHorizontal<false>,encode::Sycl<saw::encode::Native>> flow_out{1.0};
 
-		component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> flow_in{
-			[&](){
-				uint64_t target_t_i = 16u;
-				if(t_i.get() < target_t_i){
-					return 1.0 + (0.01 / target_t_i) * t_i.get();
+			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;
+				for(uint64_t i = 0u; i < Desc::D; ++i){
+					index.at({{i}}).set(idx[i]);
 				}
-				return 1.01;
-			}()
-		};
-		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;
-			for(uint64_t i = 0u; i < Desc::D; ++i){
-				index.at({{i}}).set(idx[i]);
-			}
 
-			auto info = info_f.at(index);
-			
-			switch(info.get()){
-				case 0u:
-				break;
-				case 1u:
-					bb.apply(fields,index,t_i);
-				break;
-				case 2u:
-					collision.apply(fields,macros,index,t_i);
+				auto info = info_f.at(index);
+				
+				switch(info.get()){
+					case 0u:
 					break;
-				case 3u:
-					flow_in.apply(fields,index,t_i);
-					collision.apply(fields,macros,index,t_i);
+					case 1u:
+						bb.apply(fields,index,t_i);
 					break;
-				case 4u:
-					flow_out.apply(fields,index,t_i);
-					collision.apply(fields,macros,index,t_i);
-				break;
-				default:
-				break;
-			}
-		});
-	}).wait();
+					case 2u:
+						collision.apply(fields,macros,index,t_i);
+						break;
+					case 3u:
+						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;
+				}
+			});
+		}).wait();
+	}else if constexpr ( std::is_same_v<Coll, method::FpLbm> ){
+	}
 
 	return saw::make_void();
 }