diff options
| -rw-r--r-- | examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp | 65 | ||||
| -rw-r--r-- | modules/core/c++/fplbm.hpp | 107 |
2 files changed, 100 insertions, 72 deletions
diff --git a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp index 3164945..42495c6 100644 --- a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp +++ b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp @@ -48,7 +48,8 @@ template<typename T, typename Desc> using MacroStruct = Struct< Member<VelChunk<T,Desc>, "velocity">, Member<RhoChunk<T>, "density">, - Member<ScalarChunk<T,Desc>, "porosity"> + Member<ScalarChunk<T,Desc>, "porosity">, + Member<VelChunk<T,Desc>, "force"> >; //template<typename T, typename Desc> @@ -61,13 +62,12 @@ using ParticleSpheroidGroup = ParticleGroup<T,Desc::D,1u,coll::Spheroid<T>>; } template<typename T, typename Desc> -saw::error_or<void> setup_initial_conditions( +saw::error_or<void> init( saw::data<sch::ChunkStruct<T,Desc>>& fields, saw::data<sch::MacroStruct<T,Desc>>& macros, saw::data<sch::ParticleSpheroidGroup<T,Desc>>& particles ){ auto& info_f = fields.template get<"info">(); - auto& porous_f = macros.template get<"porosity">(); // Set everything as walls iterator<Desc::D>::apply( [&](auto& index){ @@ -107,42 +107,34 @@ saw::error_or<void> setup_initial_conditions( {{0u,1u}} ); // - auto& df_f = fields.template get<"dfs_old">(); + auto& dfs_old_f = fields.template get<"dfs_old">(); + auto& dfs_f = fields.template get<"dfs">(); 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">(); iterator<Desc::D>::apply( [&](auto& index){ - auto& df = df_f.at(index); + auto& dfs = dfs_f.at(index); 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); - - df = eq; - }, - {},// 0-index - df_f.get_dims() - ); - - iterator<Desc::D>::apply( - [&](auto& index){ - auto& df = df_f.at(index); - auto& rho = rho_f.at(index); - rho.at({}) = {1}; - auto& vel = vel_f.at(index); - if(info_f.at(index).get() == 2u){ - vel.at({{0u}}) = 0.0; + for(uint64_t i{0u}; i < Desc::D; ++i){ + vel.at({{i}}).set(0); } auto eq = equilibrium<T,Desc>(rho,vel); - df = eq; + dfs = eq; + dfs_old_f.at(index) = dfs; + auto& force = force_f.at(index); + for(uint64_t i{0u}; i < Desc::D; ++i){ + force.at({{i}}).set(0); + } }, {},// 0-index - df_f.get_dims(), - {{1u,1u}} + df_f.get_dims() ); saw::data<sch::Scalar<T>> rad; @@ -181,19 +173,19 @@ saw::error_or<void> step( auto& porous_f = macros.template get<"porosity">(); q.submit([&](acpp::sycl::handler& h){ - component<T,Desc,cmpt::FplbmReset,encode::Sycl<saw::encode::Native>> hlbm_reset; + component<T,Desc,cmpt::FpLbmReset,encode::Sycl<saw::encode::Native>> fplbm_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); + fplbm_reset.apply(fields,macros,index,t_i); }); }).wait(); q.submit([&](acpp::sycl::handler& h){ - component<T,Desc,cmpt::FplbmOneParticle,encode::Sycl<saw::encode::Native>> hlbm_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; @@ -201,13 +193,13 @@ saw::error_or<void> step( index.at({{i}}).set(idx[i]); } - hlbm_one_part.apply(fields,macros,particles,index,t_i,{16u}); + fplbm_one_part.apply(fields,macros,particles,index,t_i,{1u}); }); }).wait(); // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ - component<T,Desc,cmpt::Fplbm,encode::Sycl<saw::encode::Native>> collision{0.8}; + component<T,Desc,cmpt::FpLbm,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{ @@ -221,7 +213,6 @@ saw::error_or<void> step( }; 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){ @@ -253,14 +244,6 @@ saw::error_or<void> step( }); }).wait(); - - // Step - /* - q.submit([&](acpp::sycl::handler& h){ - // h.depends_on(collision_ev); - }).wait(); - */ - return saw::make_void(); } } @@ -278,7 +261,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){ } auto& lbm_dir = eo_lbm_dir.get_value(); - auto out_dir = lbm_dir / "moving_poiseulle_particles_2d_hlbm_gpu"; + auto out_dir = lbm_dir / "moving_poiseulle_particles_2d_fplbm_gpu"; { std::error_code ec; @@ -327,7 +310,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){ sycl_q.wait(); { - auto eov = setup_initial_conditions<T,Desc>(*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_parts_data_ptr); + auto eov = init<T,Desc>(*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_parts_data_ptr); if(eov.is_error()){ return eov; } @@ -421,6 +404,8 @@ saw::error_or<void> lbm_main(int argc, char** argv){ } }); }).wait(); + + // EPoll wait.poll(); if(print_status){ std::cout<<"Status: "<<i.get()<<" of "<<time_steps.get()<<" - "<<(i.template cast_to<sch::Float64>().get() * 100 / time_steps.get())<<"%"<<std::endl; diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index 31b66c4..3c8e173 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -1,6 +1,9 @@ #pragma once #include "common.hpp" +#include "particle/aabb.hpp" +#include "particle/porosity.hpp" +#include "iterator.hpp" namespace kel { namespace lbm { @@ -28,6 +31,24 @@ public: auto& por = por_f.at(index); por.at({}) = 1.0; + + auto& force_f = macros.template get<"force">(); + auto& force = force_f.at(index); + for(uint64_t i{0u}; i < Descriptor::D; ++i){ + force.at({{i}}) = 0.0; + } + + 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& rho = rho_f.at(index); + auto& vel_f = macros.template get<"velocity">(); + auto& vel = vel_f.at(index); + + compute_rho_u<T,Descriptor>(dfs,rho,vel); } }; @@ -52,19 +73,18 @@ 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,Descriptor::D>> index, saw::data<sch::UInt64> time_step) 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">(); - 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); - compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel); + compute_rho_u<T,Descriptor>(dfs,rho,vel); auto eq = equilibrium<T,Descriptor>(rho,vel); using dfi = df_info<T,Descriptor>; @@ -72,36 +92,34 @@ public: auto& force_f = macros.template get<"force">(); auto& force = force_f.at(index); - auto& porosity_f = macros.template get<"porosity">(); - auto& porosity = porosity_f.at(index); + auto& por_f = macros.template get<"porosity">(); + auto& por = por_f.at(index); saw::data<sch::Scalar<T>> dfi_inv_cs2; dfi_inv_cs2.at({}).set(dfi::inv_cs2); - for(uint64_t i = 0u; i < Descriptor::Q; ++i){ - // saw::data<T> ci_min_u{0}; + for(uint64_t i{0u}; i < Descriptor::Q; ++i){ saw::data<sch::Vector<T,Descriptor::D>> ci; - for(uint64_t d = 0u; d < Descriptor::D; ++d){ + for(uint64_t d{0u}; d < Descriptor::D; ++d){ ci.at({{d}}).set(static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])); } auto ci_dot_u = saw::math::dot(ci,vel); - // saw::data<sch::Vector<T,Descriptor::D>> F_i; - // F_i = f * (c_i - u * ics2 + <c_i,u> * c_i * ics2 * ics2) * w_i; saw::data<sch::Scalar<T>> w; w.at({}).set(dfi::weights[i]); - auto F_i_d = saw::math::dot(force * w, (ci - vel * dfi_inv_cs2 + ci * ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2 )); - /* - saw::data<sch::Scalar<T>> F_i_sum; - for(uint64_t d = 0u; d < Descriptor::D; ++d){ - saw::data<sch::Scalar<T>> F_i_d; - F_i_d.at({}) = F_i.at({{d}}); - F_i_sum = F_i_sum + F_i_d; - } - */ + auto term1 = (ci-vel) * dfi_inv_cs2; + auto term2 = ci * (ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2); + + auto force_projection = saw::math::dot(term1 + term2, force); + + auto F_i = w * force_projection; - dfs.at({i}) = dfs.at({i}) + frequency_ * (eq.at(i) - dfs.at({i}) ) + F_i_d.at({}) * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_); + dfs.at({i}) = dfs.at({i}) + + frequency_ + * (eq.at(i) - dfs.at({i}) ) + + F_i.at({}) + * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_); } } }; @@ -119,14 +137,12 @@ public: // 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& 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">(); - - // Temporary find a better way auto& force_f = macros.template get<"force">(); /** @@ -136,8 +152,19 @@ public: // TODO - Change to tuple later auto parts = pg.template get<"particles">(); - auto parts_size = parts.meta.at({0u}); - auto pi = parts.at(index); + 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">(); + + saw::data<sch::Scalar<T>> sss; + sss.at({}) = sub_steps.template cast_to<T>(); + auto vel_s = (pirb_pos-pirb_pos_old) / sss; + + 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()){ @@ -145,9 +172,28 @@ public: } auto& aabb = eo_aabb.get_value(); - iterator<Desc::D>::apply([&](const auto& index_f) -> void { + saw::data<sch::Scalar<T>> two; + two.at({}).set(2); - }, start, stop); + saw::data<sch::Scalar<T>> one; + one.at({}) = 1.0; + + saw::data<sch::Scalar<T>> eps; + eps.at({}) = 1.5f; + + 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; + + force = ( vel_s * rho - vel * rho ) * two * flip_por; + }, aabb.template get<"a">(), aabb.template get<"b">()); } }; @@ -208,15 +254,12 @@ public: saw::data<sch::Scalar<T>> min_two; min_two.at({}).set(-2); - // Maybe ? auto& force = force_f.at(index); force = vel * rho * min_two * flip_porosity; saw::data<sch::Scalar<T>> dfi_inv_cs2; dfi_inv_cs2.at({}).set(dfi::inv_cs2); - // auto vel = vel_f.at(index); - for(uint64_t i = 0u; i < Descriptor::Q; ++i){ // saw::data<T> ci_min_u{0}; saw::data<sch::Vector<T,Descriptor::D>> ci; |