diff options
| -rw-r--r-- | c++/collision.hpp | 17 | ||||
| -rw-r--r-- | c++/examples/cavity_2d.cpp | 12 | ||||
| -rw-r--r-- | c++/write_vtk.hpp | 112 |
3 files changed, 138 insertions, 3 deletions
diff --git a/c++/collision.hpp b/c++/collision.hpp index 87187c2..ba35be1 100644 --- a/c++/collision.hpp +++ b/c++/collision.hpp @@ -14,15 +14,30 @@ class component<T, Descriptor, cmpt::BGK> { public: using Component = cmpt::BGK; + /** + * Thoughts regarding automating order setup + */ using access = cmpt::access_tuple< cmpt::access<"dfs", 1, true, true, true> >; + /** + * 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"; - void apply(saw::data<sch::CellField>& field){ + /** + * Raw setup + */ + void apply(saw::data<sch::CellField>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, uint64_t time_step){ + bool is_even = ((time_step % 2) == 0); + auto& cell = field.at(index); + + auto& dfs_old = is_even ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + auto& dfs = (!is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + } }; diff --git a/c++/examples/cavity_2d.cpp b/c++/examples/cavity_2d.cpp index 88a7702..0a8e9b1 100644 --- a/c++/examples/cavity_2d.cpp +++ b/c++/examples/cavity_2d.cpp @@ -71,6 +71,9 @@ void compute_rho_u ( vel[1] /= rho; } +/** + * Unsure if feasible. I would have a rho normalization on the dfs and then I would use the const rho_computation + */ template<typename Desc> void compute_const_rho_u ( saw::data<sch::DfCell<Desc>>& dfs, @@ -80,15 +83,20 @@ void compute_const_rho_u ( { using dfi = df_info<sch::T, Desc>; + saw::native_data_type<sch::T>::type real_rho = 0; std::fill(vel.begin(), vel.end(), 0); for(size_t i = 0; i < Desc::Q; ++i){ + real_rho += dfs(i).get(); vel[0] += dfi::directions[i][0] * dfs(i).get(); vel[1] += dfi::directions[i][1] * dfs(i).get(); } + for(size_t i = 0; i < Desc::Q; ++i){ + dfs(i).set(dfs(i).get() * (rho/real_rho)); + } - vel[0] /= rho; - vel[1] /= rho; + vel[0] *= real_rho / (rho*rho) + vel[1] *= real_rho / (rho*rho) } /** diff --git a/c++/write_vtk.hpp b/c++/write_vtk.hpp new file mode 100644 index 0000000..7bde854 --- /dev/null +++ b/c++/write_vtk.hpp @@ -0,0 +1,112 @@ +#pragma once + +#include <forstio/error.hpp> + +namespace kel { +namespace lbm { +namespace impl { + +template<typename... MemberTypes, saw::string_literal... MemberNames> +struct lbm_vtk_writer<sch::Struct<MemberTypes,MemberNames>...> { + + saw::error_or<void> apply(){ + return saw::make_void(); + } +}; + +template<typename CellFieldSchema> +struct lbm_vtk_writer { +}; + +template<typename Desc, typename... StructT, saw::string_literal... StructN> +struct lbm_vtk_writer<CellField<Desc,Struct<Member<StructT,StructN>...>>> { + template<uint64_t i> + saw::error_or<void> write_i(std::ofstream& vtk_file, const saw::data<CellField<Desc,Struct<Member<StructT,StructN>...>>>& field){ + + // vtk_file<< + + return saw::make_void(); + } + + template<uint64_t i> + saw::error_or<void> iterate_i(std::ofstream& vtk_file, const saw::data<CellField<Desc,Struct<Member<StructT,StructN>...>>>& field){ + { + auto eov = write_i<i>(vtk_file, field); + if(eov.is_error()){ + return eov; + } + } + if constexpr ( (i+1u) < sizeof...(StructT) ){ + return iterate_i<i+1u>(vtk_file, field); + } + return saw::make_void(); + } + + saw::error_or<void> apply(std::ofstream& vtk_file, + const saw::data<CellField<Desc,Struct<Member<StructT,StructN>...>>>& field){ + + auto meta = field.meta(); + // DIMENSIONS + { + vtk_file << "DIMENSIONS "; + for(saw::data<sch::UInt64> i{0u}; i < Desc::D; ++i){ + vtk_file << " " << meta.at(i).get(); + } + for(saw::data<sch::UInt64> i{Desc::D}; i < 3u; ++i){ + vtk_file << " 1"; + } + + vtk_file << "\n"; + } + + // HEADER TO BODY + { + vtk_file << "\n"; + } + + if constexpr (sizeof...(StructT) > 0u){ + return iterate_i<0u>(vtk_file, field); + } + + return saw::make_void(); + } +}; +} + +template<typename Struct> +saw::error_or<void> write_vtk(const std::string_view file_name, const saw::data<Struct>& field){ + + std::string vtk_file_name{file_name}; + std::ofstream vtk_file{vtk_file_name}; + + if(!vtk_file.is_open()){ + return saw::make_error<saw::err::critical>("Could not open file."); + } + + vtk_file + <<"# vtk DataFile Version 3.0\n" + <<"LBM File\n" + <<"ASCII\n" + <<"DATASET STRUCTURED POINTS\n" + <<"SPACING 1.0 1.0 1.0\n" + <<"ORIGIN 0.0 0.0 0.0\n" + ; + + + auto eov = lbm_vtk_writer::apply(vtk_file, field); + return eov; + /* + vtk_file <<"# vtk DataFile Version 3.0\n"; + vtk_file <<"Velocity Cavity2D example\n"; + vtk_file <<"ASCII\n"; + vtk_file <<"DATASET STRUCTURED_POINTS\n"; + vtk_file <<"DIMENSIONS "<< dim_x <<" "<<dim_y<<" 1\n"; + vtk_file <<"SPACING 1.0 1.0 1.0\n"; + vtk_file <<"ORIGIN 0.0 0.0 0.0\n"; + vtk_file <<"POINT_DATA "<<(dim_x*dim_y)<<"\n"; + vtk_file <<"VECTORS Velocity float\n"; + */ +:w +} +} +} |