diff options
| -rw-r--r-- | c++/collision.hpp | 70 | ||||
| -rw-r--r-- | examples/cavity_2d_gpu.cpp | 1 | ||||
| -rw-r--r-- | examples/poiseulle_particles_channel_2d.cpp | 6 |
3 files changed, 72 insertions, 5 deletions
diff --git a/c++/collision.hpp b/c++/collision.hpp index d9eb2c0..df9a734 100644 --- a/c++/collision.hpp +++ b/c++/collision.hpp @@ -7,6 +7,7 @@ namespace kel { namespace lbm { namespace cmpt { struct BGK {}; +struct BGKGuo {}; } /** @@ -16,9 +17,11 @@ template<typename T, typename Descriptor> class component<T, Descriptor, cmpt::BGK> { private: typename saw::native_data_type<T>::type relaxation_; + typename saw::native_data_type<T>::type frequency_; public: component(typename saw::native_data_type<T>::type relaxation__): - relaxation_{relaxation__} + relaxation_{relaxation__}, + frequency_{1 / relaxation_} {} using Component = cmpt::BGK; @@ -46,7 +49,7 @@ public: auto& cell = field(index); auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); - auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + // auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); saw::data<T> rho; saw::data<sch::FixedArray<T,Descriptor::D>> vel; @@ -58,5 +61,68 @@ public: } } }; + +template<typename T, typename Descriptor> +class component<T, Descriptor, cmpt::BGKGuo> { +private: + typename saw::native_data_type<T>::type relaxation_; + typename saw::native_data_type<T>::type frequency_; +public: + component(typename saw::native_data_type<T>::type relaxation__): + relaxation_{relaxation__}, + frequency_{1 / relaxation_} + {} + + using Component = cmpt::BGKGuo; + + using access = cmpt::access_tuple< + cmpt::access<"dfs", 1, true, true, true>, + cmpt::access<"force", 0, true, false> + >; + + /** + * Thoughts regarding automating order setup + */ + static constexpr saw::string_literal name = "collision"; + static constexpr saw::string_literal after = ""; + static constexpr saw::string_literal before = "streaming"; + + template<typename CellFieldSchema> + void apply(saw::data<CellFieldSchema>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, uint64_t time_step){ + bool is_even = ((time_step & 2) == 0); + auto& cell = field(index); + + auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + // auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + + saw::data<T> rho; + saw::data<sch::FixedArray<T,Descriptor::D>> vel; + compute_rho_u<T,Descriptor>(dfs_old,rho,vel); + auto eq = equilibrium<T,Descriptor>(rho,vel); + + using dfi = df_info<T,Descriptor>; + + auto& force = cell.template get<"force">(); + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + saw::data<T> ci_dot_u{0}; + saw::data<T> ci_dot_F{0}; + saw::data<T> u_dot_F{0}; + for(uint64_t d = 0u; d < Descriptor::D; ++d){ + ci_dot_u = ci_dot_u + vel.at({d}) * dfi::directions[i][d]; + ci_dot_F = ci_dot_F + force({d}) * dfi::directions[i][d]; + u_dot_F = u_dot_F + vel.at({d}) * force({d}); + } + + // auto df_old = dfs_old({i}); + saw::data<T> F_i = saw::data<T>{dfi::weights[i]} * (ci_dot_F / dfi::cs2 + ci_dot_u * ci_dot_F / (dfi::cs2 * dfi::cs2)); + if(F_i.get() != 0){ + std::cerr<<"AAAH"; + } + + dfs_old({i}) = dfs_old({i}) + saw::data<T>{frequency_} * (eq.at(i) - dfs_old({i}) ) - F_i * saw::data<T>{1.0f - 0.5f * frequency_}; + } + } +}; } } diff --git a/examples/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu.cpp index b24ca38..19563e2 100644 --- a/examples/cavity_2d_gpu.cpp +++ b/examples/cavity_2d_gpu.cpp @@ -199,7 +199,6 @@ void lbm_step( */ component<sch::T, sch::D2Q9, cmpt::BGK> coll{0.59}; component<sch::T, sch::D2Q9, cmpt::BounceBack> bb; - component<sch::T, sch::D2Q9, cmpt::MovingWall> bb_lid; bb_lid.lid_vel = {0.1,0.0}; diff --git a/examples/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d.cpp index e495043..fed78b2 100644 --- a/examples/poiseulle_particles_channel_2d.cpp +++ b/examples/poiseulle_particles_channel_2d.cpp @@ -643,7 +643,9 @@ int main(int argc, char** argv){ compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel); rho = rho * saw::data<sch::T>{dfi::cs2}; - force = + for(uint64_t d = 0u; d < sch::D2Q9::D; ++d){ + force.at({d}) = cell.template get<"force">()({d}); + } part_mask.set(0u); }, {{0u,0u}}, meta); @@ -683,7 +685,7 @@ int main(int argc, char** argv){ // Cast down to get lower corner. // Before casting shift by 0.5 for closest pick saw::data<sch::FixedArray<sch::UInt64,2u>> p_cell_pos {{ - static_cast<uint64_t>(p_pos_lie.at(saw::data<T>{{0u}}).get()+0.5), + static_cast<uint64_t>(p_pos_lie.at({{0u}}).get()+0.5), static_cast<uint64_t>(p_pos_lie.at({{1u}}).get()+0.5) }}; |