summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--examples/one_particle_sedimentation_2d/common.hpp3
-rw-r--r--examples/one_particle_sedimentation_2d/sim.cpp4
-rw-r--r--examples/one_particle_sedimentation_2d/step.hpp33
-rw-r--r--examples/schaefer_turek_durst_krause_rannbacher/common.hpp4
-rw-r--r--examples/schaefer_turek_durst_krause_rannbacher/sim.cpp6
-rw-r--r--examples/stdkr_2d_bgk/common.hpp4
-rw-r--r--modules/core/c++/fplbm.hpp22
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">());
generated by cgit v1.2.3 (git 2.25.1) at 2026-08-30 08:07:14 +0000