diff options
Diffstat (limited to 'lib/core')
| -rw-r--r-- | lib/core/c++/collision.hpp | 22 | ||||
| -rw-r--r-- | lib/core/c++/equilibrium.hpp | 8 | ||||
| -rw-r--r-- | lib/core/c++/macroscopic.hpp | 93 | ||||
| -rw-r--r-- | lib/core/c++/stream.hpp | 27 | ||||
| -rw-r--r-- | lib/core/c++/write_vtk.hpp | 18 | ||||
| -rw-r--r-- | lib/core/tests/equilibrium.cpp | 7 |
6 files changed, 79 insertions, 96 deletions
diff --git a/lib/core/c++/collision.hpp b/lib/core/c++/collision.hpp index ed26a08..a05e263 100644 --- a/lib/core/c++/collision.hpp +++ b/lib/core/c++/collision.hpp @@ -61,6 +61,28 @@ public: dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i})); } } + + 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 { + bool is_even = ((time_step.get() % 2) == 0); + + // auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + + saw::data<sch::Scalar<T>>& rho = rho_f.at(index); + saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index); + + compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel); + + auto eq = equilibrium<T,Descriptor>(rho,vel); + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i})); + } + } }; template<typename T, typename Descriptor, typename Encode> diff --git a/lib/core/c++/equilibrium.hpp b/lib/core/c++/equilibrium.hpp index bb55d00..7d1324e 100644 --- a/lib/core/c++/equilibrium.hpp +++ b/lib/core/c++/equilibrium.hpp @@ -5,7 +5,7 @@ namespace kel { namespace lbm { template<typename T, typename Descriptor> -saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<T> rho, saw::data<sch::FixedArray<T,Descriptor::D>> vel){ +saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<sch::Scalar<T>> rho, saw::data<sch::Vector<T,Descriptor::D>> vel){ using dfi = df_info<T, Descriptor>; saw::data<sch::FixedArray<T,Descriptor::Q>> eq; @@ -17,7 +17,7 @@ saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<T> rho, saw:: // Velocity * Velocity meaning || vel ||_2^2 or <vel,vel>_2 saw::data<T> vel_vel{0.0}; for(uint64_t j = 0u; j < Descriptor::D; ++j){ - vel_vel = vel_vel + vel.at(j) * vel.at(j); + vel_vel = vel_vel + vel.at({{j}}) * vel.at({{j}}); } /** @@ -27,13 +27,13 @@ saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<T> rho, saw:: saw::data<T> vel_c{}; for(uint64_t j = 0u; j < Descriptor::D; ++j){ // <vel,c_i>_2 - vel_c = vel_c + (vel.at(j) * saw::data<T>{static_cast<saw::native_data_type<T>::type>(dfi::directions[i][j])}); + vel_c = vel_c + (vel.at({{j}}) * saw::data<T>{static_cast<saw::native_data_type<T>::type>(dfi::directions[i][j])}); } auto vel_c_cs2 = vel_c * saw::data<T>{dfi::inv_cs2}; eq.at(i).set( - dfi::weights[i] * rho.get() * + dfi::weights[i] * rho.at({}).get() * ( 1.0 + vel_c_cs2.get() diff --git a/lib/core/c++/macroscopic.hpp b/lib/core/c++/macroscopic.hpp index 1532873..19ab3e9 100644 --- a/lib/core/c++/macroscopic.hpp +++ b/lib/core/c++/macroscopic.hpp @@ -4,113 +4,32 @@ namespace kel { namespace lbm { -template<typename T, typename Desc> -void compute_rho_u ( - const saw::data<sch::FixedArray<T,Desc::Q>>& dfs, - saw::data<T>& rho, - saw::data<sch::FixedArray<T,Desc::D>>& vel - ) -{ - using dfi = df_info<T, Desc>; - - rho = 0; - for(size_t j = 0; j < Desc::D; ++j){ - vel.at({{j}}).set(0); - } - - for(size_t j = 0; j < Desc::Q; ++j){ - rho = rho + dfs.at({j}); - for(size_t i = 0; i < Desc::D; ++i){ - vel.at({{i}}) = vel.at({{i}}) + dfi::directions[j][i] * dfs.at({j}).get(); - } - } - - for(size_t i = 0; i < Desc::D; ++i){ - vel.at({i}) = vel.at({i}) / rho; - } -} -/** - * Calculate the macroscopic variables rho and u in Lattice Units. - */ -template<typename T, typename Desc> -void compute_rho_u ( - const saw::data<sch::Cell<T, Desc, 0, 0, 1>>& dfs, - typename saw::native_data_type<T>::type& rho, - std::array<typename saw::native_data_type<T>::type, 2>& vel - ) -{ - using dfi = df_info<T, Desc>; - - rho = 0; - std::fill(vel.begin(), vel.end(), 0); - - for(size_t j = 0; j < Desc::Q; ++j){ - rho += dfs(j).get(); - for(size_t i = 0; i < Desc::D; ++i){ - vel[i] += dfi::directions[j][i] * dfs(j).get(); - } - } - - for(size_t i = 0; i < Desc::D; ++i){ - vel[i] /= rho; - } -} - /** * Calculate the macroscopic variables rho and u in Lattice Units. */ template<typename T, typename Desc> void compute_rho_u ( - const saw::data<sch::Cell<T, Desc, 0, 0, 1>>& dfs, - saw::ref<saw::data<T>> rho, - saw::ref<saw::data<sch::FixedArray<T,Desc::D>>> vel - ) -{ - using dfi = df_info<T, Desc>; - - rho().set(0); - for(uint64_t i = 0; i < Desc::D; ++i){ - vel().at({i}).set(0); - } - - for(size_t j = 0; j < Desc::Q; ++j){ - rho() = rho() + dfs(j); - for(size_t i = 0; i < Desc::D; ++i){ - vel().at({i}) = vel().at({i}) + saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[j][i])} * dfs(j); - } - } - - for(size_t i = 0; i < Desc::D; ++i){ - vel().at({i}) = vel().at({i}) / rho(); - } -} - -/** - * Calculate the macroscopic variables rho and u in Lattice Units. - */ -template<typename T, typename Desc> -void compute_rho_u ( - const saw::data<sch::Cell<T, Desc, 0, 0, 1>>& dfs, - saw::ref<saw::data<T>> rho, + const saw::data<sch::FixedArray<T,Desc::Q>>& dfs, + saw::ref<saw::data<sch::Scalar<T>>> rho, saw::ref<saw::data<sch::Vector<T,Desc::D>>> vel ) { using dfi = df_info<T, Desc>; - rho().set(0); + rho().at({}).set(0); for(uint64_t i = 0; i < Desc::D; ++i){ vel().at({{i}}).set(0); } for(size_t j = 0; j < Desc::Q; ++j){ - rho() = rho() + dfs(j); + rho().at({}) = rho().at({}) + dfs.at({j}); for(size_t i = 0; i < Desc::D; ++i){ - vel().at({{i}}) = vel().at({{i}}) + saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[j][i])} * dfs(j); + vel().at({{i}}) = vel().at({{i}}) + saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[j][i])} * dfs.at({j}); } } for(size_t i = 0; i < Desc::D; ++i){ - vel().at({{i}}) = vel().at({{i}}) / rho(); + vel().at({{i}}) = vel().at({{i}}) / rho().at({}); } } } diff --git a/lib/core/c++/stream.hpp b/lib/core/c++/stream.hpp index 8c0342b..dc7cfb3 100644 --- a/lib/core/c++/stream.hpp +++ b/lib/core/c++/stream.hpp @@ -17,19 +17,44 @@ public: static constexpr saw::string_literal before = ""; template<typename CellFieldSchema> - void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step) const { + void apply(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + using dfi = df_info<T,Descriptor>; + auto g_index = index; + + for(uint64_t i = 0u; i < Descriptor::D; ++i){ + ++g_index.at({{i}}); + } + 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_new_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); auto info_f = field.template get<"info">(); auto info_meta = info_f.get_dims(); + /* bool border = false; for(uint64_t i = 0u; i < Descriptor::D; ++i){ auto ind_i = index.at({i}); border |= (ind_i.get()) == 0u or (ind_i == info_meta.at({i})); } + + if (not border){ + } + */ + + auto& dfs_new = dfs_new_f.ghost_at(g_index); + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + auto dir = dfi::directions[dfi::opposite_index[i]]; + auto g_index_nb = g_index; + for(uint64_t z = 0u; z < Descriptor::D; ++z){ + g_index_nb.at({z}).set(g_index.at({z}).get() + dir[z]); + } + + dfs_new.at({i}) = dfs_old_f.ghost_at(g_index_nb).at({i}); + } } }; } diff --git a/lib/core/c++/write_vtk.hpp b/lib/core/c++/write_vtk.hpp index e84cf9d..3c60a66 100644 --- a/lib/core/c++/write_vtk.hpp +++ b/lib/core/c++/write_vtk.hpp @@ -100,6 +100,8 @@ struct lbm_vtk_writer<sch::Chunk<T,Ghost,D...>> { } } + vtk_file<<"\n"; + return saw::make_void(); } }; @@ -129,6 +131,21 @@ struct lbm_vtk_writer<sch::Vector<T,D>> { } }; +template<typename T> +struct lbm_vtk_writer<sch::Scalar<T>> { + static saw::error_or<void> apply_header(std::ostream& vtk_file, std::string_view name){ + vtk_file<<"SCALARS "<<name<<" float\n"; + vtk_file<<"LOOKUP_TABLE default\n"; + return saw::make_void(); + } + static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::Scalar<T>>& field){ + // vtk_file<<"VECTORS "<<name<<" float\n"; + vtk_file<<field.at({}).get(); + vtk_file<<"\n"; + return saw::make_void(); + } +}; + template<typename... MemberT, saw::string_literal... Keys, uint64_t... Ghost, uint64_t... Dims> struct lbm_vtk_writer<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>> final { template<uint64_t i> @@ -176,7 +193,6 @@ struct lbm_vtk_writer<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>, static_assert(saw::ct_multiply<uint64_t,Dims...>::value > 0u, "Invalid Dim size resulting in length 0u"); for(saw::data<sch::UInt64> i{0u}; i.get() < sizeof...(Dims); ++i){ - pd_size = pd_size * meta.at(i); vtk_file << " " << meta.at(i).get(); } for(saw::data<sch::UInt64> i{sizeof...(Dims)}; i.get() < 3u; ++i){ diff --git a/lib/core/tests/equilibrium.cpp b/lib/core/tests/equilibrium.cpp index 9201e55..20a1f08 100644 --- a/lib/core/tests/equilibrium.cpp +++ b/lib/core/tests/equilibrium.cpp @@ -11,10 +11,11 @@ void check_equilibrium(){ using dfi = lbm::df_info<lbm::sch::Float64,Descriptor>; - saw::data<lbm::sch::Float64> rho{1.0}; - saw::data<lbm::sch::FixedArray<lbm::sch::Float64,Descriptor::D>> vel; + saw::data<lbm::sch::Scalar<lbm::sch::Float64>> rho; + rho.at({}).set(1.0); + saw::data<lbm::sch::Vector<lbm::sch::Float64,Descriptor::D>> vel; for(saw::data<lbm::sch::UInt64> i{0u}; i.get() < Descriptor::D; ++i){ - vel.at(i) = {0.0}; + vel.at({{i}}) = {0.0}; } auto eq = lbm::equilibrium<lbm::sch::Float64,Descriptor>(rho,vel); |