diff options
| author | Claudius "keldu" Holeksa <mail@keldu.de> | 2026-08-14 23:43:54 +0200 |
|---|---|---|
| committer | Claudius "keldu" Holeksa <mail@keldu.de> | 2026-08-14 23:43:54 +0200 |
| commit | 0e9322fe0d024a06f23430285a055c16cf1e3eeb (patch) | |
| tree | 657a33098c96b14c0e09748c6b9d4a7e9cd16160 | |
| parent | 0d15416f60491bd4d17bd8b45a9487c98b8132b0 (diff) | |
| parent | c634aac438972b7bafb5fbbeba7962d174cc0d03 (diff) | |
| download | libs-lbm-0e9322fe0d024a06f23430285a055c16cf1e3eeb.tar.gz | |
Merge branch 'dev'
47 files changed, 1716 insertions, 205 deletions
diff --git a/default.nix b/default.nix index f560d7e..d7aeb9c 100644 --- a/default.nix +++ b/default.nix @@ -218,6 +218,11 @@ in rec { inherit kel; }; + stdkr_2d = pkgs.callPackage ./examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix { + inherit pname version stdenv forstio adaptive-cpp; + inherit kel; + }; + heterogeneous_computing = pkgs.callPackage ./examples/heterogeneous_computing/.nix/derivation.nix { inherit pname version stdenv forstio adaptive-cpp; inherit kel; diff --git a/examples/cavity_2d/cavity_2d.cpp b/examples/cavity_2d/cavity_2d.cpp index 702661a..af22796 100644 --- a/examples/cavity_2d/cavity_2d.cpp +++ b/examples/cavity_2d/cavity_2d.cpp @@ -42,7 +42,7 @@ using CellStruct = Struct< template<typename T, uint64_t D> using MacroStruct = Struct< - Member<FixedArray<T,D>, "velocity">, + Member<FixedArray<T,D>, "momentum">, Member<T, "pressure"> >; @@ -406,7 +406,7 @@ int main(){ auto dfs = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); auto& rho = macros.at({{i,j}}).template get<"pressure">(); - auto& vel = macros.at({{i,j}}).template get<"velocity">(); + auto& vel = macros.at({{i,j}}).template get<"momentum">(); compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel); }, lattice); diff --git a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp index 0fc490a..fe44baf 100644 --- a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp +++ b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp @@ -47,7 +47,7 @@ using CellStruct = CellFieldStruct< template<typename T, uint64_t D> using MacroStruct = Struct< - Member<FixedArray<T,D>, "velocity">, + Member<FixedArray<T,D>, "momentum">, Member<T, "pressure"> >; @@ -362,7 +362,7 @@ int main(){ auto dfs_field = even_step ? dfs_old : dfs; auto& rho = macros.at(index).template get<"pressure">(); - auto& vel = macros.at(index).template get<"velocity">(); + auto& vel = macros.at(index).template get<"momentum">(); compute_rho_u<sch::T,sch::D2Q9>(dfs_field[j],rho,vel); }, {{0u,0u}}, meta); std::string vtk_f_name{"tmp/cavity_2d_gpu_"}; diff --git a/examples/cavity_3d.cpp b/examples/cavity_3d.cpp index 64356be..5ae4834 100644 --- a/examples/cavity_3d.cpp +++ b/examples/cavity_3d.cpp @@ -46,7 +46,7 @@ using CellStruct = Struct< template<typename T, uint64_t D> using MacroStruct = Struct< - Member<FixedArray<T,D>, "velocity">, + Member<FixedArray<T,D>, "momentum">, Member<T, "pressure"> >; @@ -411,7 +411,7 @@ int main(){ auto dfs = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); auto& rho = macros.at({{i,j}}).template get<"pressure">(); - auto& vel = macros.at({{i,j}}).template get<"velocity">(); + auto& vel = macros.at({{i,j}}).template get<"momentum">(); compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel); }, lattice); diff --git a/examples/meta/meta.cpp b/examples/meta/meta.cpp index a5e2ae5..7232c61 100644 --- a/examples/meta/meta.cpp +++ b/examples/meta/meta.cpp @@ -75,7 +75,7 @@ saw::error_or<void> meta_main(int argc, char** argv){ meta_lbm_info<sch::Prec,sch::Descriptor<1u,3u>>(args_str); break; default: - return saw::make_error<saw::err::critical>("Second positional needs to be a valid velocity set"); + return saw::make_error<saw::err::critical>("Second positional needs to be a valid momentum set"); } } break; @@ -89,7 +89,7 @@ saw::error_or<void> meta_main(int argc, char** argv){ meta_lbm_info<sch::Prec,sch::Descriptor<2u,9u>>(args_str); break; default: - return saw::make_error<saw::err::critical>("Second positional needs to be a valid velocity set"); + return saw::make_error<saw::err::critical>("Second positional needs to be a valid momentum set"); } } break; @@ -108,7 +108,7 @@ saw::error_or<void> meta_main(int argc, char** argv){ meta_lbm_info<sch::Prec,sch::Descriptor<3u,27u>>(args_str); break; default: - return saw::make_error<saw::err::critical>("Second positional needs to be a valid velocity set"); + return saw::make_error<saw::err::critical>("Second positional needs to be a valid momentum set"); } } break; diff --git a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp index cc7649a..4c9c8bd 100644 --- a/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp +++ b/examples/moving_poiseulle_particles_2d_fplbm_gpu/sim.cpp @@ -55,7 +55,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<ScalarChunk<T,Desc>, "porosity">, Member<VelChunk<T,Desc>, "force"> @@ -119,7 +119,7 @@ saw::error_or<void> init( 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& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); auto& force_f = macros.template get<"force">(); diff --git a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp index 9e812e2..bd754bf 100644 --- a/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp +++ b/examples/moving_poiseulle_particles_2d_hlbm_gpu/sim.cpp @@ -46,7 +46,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<ScalarChunk<T,Desc>, "porosity"> >; @@ -109,7 +109,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator<Desc::D>::apply( diff --git a/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp index 2c4c9d1..2daa452 100644 --- a/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp +++ b/examples/moving_poiseulle_particles_2d_psm_gpu/sim.cpp @@ -46,7 +46,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<ScalarChunk<T,Desc>, "porosity">, Member<VelChunk<T,Desc>, "force"> @@ -110,7 +110,7 @@ saw::error_or<void> init( 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& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); auto& force_f = macros.template get<"force">(); diff --git a/examples/poiseulle_2d/poiseulle_2d.cpp b/examples/poiseulle_2d/poiseulle_2d.cpp index 3bc7a6c..f9db3df 100644 --- a/examples/poiseulle_2d/poiseulle_2d.cpp +++ b/examples/poiseulle_2d/poiseulle_2d.cpp @@ -41,7 +41,7 @@ using CellStruct = Struct< template<typename T, uint64_t D> using MacroStruct = Struct< - Member<FixedArray<T,D>, "velocity">, + Member<FixedArray<T,D>, "momentum">, Member<T, "pressure"> >; @@ -375,7 +375,7 @@ int main(int argc, char** argv){ auto& dfs = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">(); auto& rho = macros.at(index).template get<"pressure">(); - auto& vel = macros.at(index).template get<"velocity">(); + auto& vel = macros.at(index).template get<"momentum">(); compute_rho_u<sch::T,sch::D2Q9>(dfs,rho,vel); rho = rho * saw::data<sch::T>{dfi::cs2}; diff --git a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp index 1aaa5c0..9963c59 100644 --- a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp +++ b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp @@ -36,7 +36,7 @@ using CellStruct = Struct< >; using MacroStruct = Struct< - Member<FixedArray<Float64,D2Q9::D>, "velocity">, + Member<FixedArray<Float64,D2Q9::D>, "momentum">, Member<Float64, "pressure"> >; @@ -245,7 +245,7 @@ void step( auto& macro_c = macro_cells[acc_id]; saw::data<sch::T>& rho = macro_c.template get<"pressure">(); - saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"velocity">(); + saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"momentum">(); compute_rho_u<sch::T,Desc>(dfs_old,rho,vel); auto eq = equilibrium<sch::T,Desc>(rho,vel); @@ -262,7 +262,7 @@ void step( auto& macro_c = macro_cells[acc_id]; saw::data<sch::T>& rho = macro_c.template get<"pressure">(); - saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"velocity">(); + saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"momentum">(); compute_rho_u<sch::T,Desc>(dfs_old,rho,vel); auto eq = equilibrium<sch::T,Desc>(rho,vel); @@ -279,7 +279,7 @@ void step( auto& macro_c = macro_cells[acc_id]; saw::data<sch::T>& rho = macro_c.template get<"pressure">(); - saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"velocity">(); + saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"momentum">(); compute_rho_u<sch::T,Desc>(dfs_old,rho,vel); auto eq = equilibrium<sch::T,Desc>(rho,vel); diff --git a/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp b/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp index 1693c5e..1af0023 100644 --- a/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp +++ b/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp @@ -47,7 +47,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<ScalarChunk<T,Desc>, "porosity"> >; @@ -123,7 +123,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator<Desc::D>::apply( diff --git a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp index 84ce22b..9ed1646 100644 --- a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp @@ -45,7 +45,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density"> >; @@ -102,7 +102,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); iterator<Desc::D>::apply( [&](auto& index){ diff --git a/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp index c6aac92..f45e10e 100644 --- a/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp @@ -49,7 +49,7 @@ using ForceChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<VelChunk<T,Desc>, "force">, Member<ScalarChunk<T,Desc>, "porosity"> @@ -109,7 +109,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator<Desc::D>::apply( diff --git a/examples/poiseulle_particles_2d_gpu/common.hpp b/examples/poiseulle_particles_2d_gpu/common.hpp index 90b7687..b321344 100644 --- a/examples/poiseulle_particles_2d_gpu/common.hpp +++ b/examples/poiseulle_particles_2d_gpu/common.hpp @@ -49,7 +49,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<ScalarChunk<T,Desc>, "porosity"> >; diff --git a/examples/poiseulle_particles_2d_gpu/init.hpp b/examples/poiseulle_particles_2d_gpu/init.hpp index cdce25e..b7fdeb0 100644 --- a/examples/poiseulle_particles_2d_gpu/init.hpp +++ b/examples/poiseulle_particles_2d_gpu/init.hpp @@ -57,7 +57,7 @@ saw::error_or<void> init( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator<Desc::D>::apply( diff --git a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp index 6cadbd0..0061fa2 100644 --- a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp @@ -46,7 +46,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<VelChunk<T,Desc>, "force">, Member<ScalarChunk<T,Desc>, "porosity"> @@ -106,7 +106,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator<Desc::D>::apply( diff --git a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp index 640c9ab..d63515d 100644 --- a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp @@ -46,7 +46,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<VelChunk<T,Desc>, "force"> >; @@ -113,7 +113,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& force_f = macros.template get<"force">(); iterator<Desc::D>::apply( @@ -218,7 +218,7 @@ saw::error_or<void> step( q.submit([&](acpp::sycl::handler& h){ h.parallel_for(acpp::sycl::range<1u>{1u}, [=](acpp::sycl::id<1u> idx){ - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& dense_f = macros.template get<"density">(); auto& ps = particles; diff --git a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp index c2ab675..760f70a 100644 --- a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp @@ -47,7 +47,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<VelChunk<T,Desc>, "force">, Member<ScalarChunk<T,Desc>, "porosity"> @@ -124,7 +124,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator<Desc::D>::apply( diff --git a/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp index 1e6e75f..24832f2 100644 --- a/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp +++ b/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp @@ -68,7 +68,7 @@ using CellStruct = Struct< template<typename T, uint64_t D> using MacroStruct = Struct< - Member<Vector<T,D>, "velocity">, + Member<Vector<T,D>, "momentum">, Member<T, "pressure">, Member<UInt16, "particle">, Member<FixedArray<T,D>, "force"> @@ -369,7 +369,7 @@ void couple_particles_to_lattice( // Interpolate this from close U cells. // For now pick the closest U - auto& closest_u = macros.at(p_cell_pos).template get<"velocity">(); + auto& closest_u = macros.at(p_cell_pos).template get<"momentum">(); // auto p_shift = closest_u; // p_pos - p_pos_old @@ -380,7 +380,7 @@ void couple_particles_to_lattice( p_vel_rel.at({{1u}}) = p_vel_rel.at({{1u}}) / p_mass; } - // Add relative velocity to particle acceleration (Since our timestep is 1 we don't need to do anything else) + // Add relative momentum to particle acceleration (Since our timestep is 1 we don't need to do anything else) p_acc = p_acc + p_vel_rel; @@ -463,7 +463,7 @@ void couple_particles_to_lattice( } - /* So I want to include the relative velocity from both particles, but this is a bit hard considering I just assumed 0 velocity from the other party + /* So I want to include the relative momentum from both particles, but this is a bit hard considering I just assumed 0 momentum from the other party else if(n_macro_cell_particle.get() != (i.get() + 1u) and n_macro_cell_particle.get() > 0u){ // Generally compare } @@ -721,7 +721,7 @@ saw::error_or<void> real_main(int argc, char** argv){ auto& macro_cell = macros.at(index); auto& rho = macro_cell.template get<"pressure">(); - auto& vel = macro_cell.template get<"velocity">(); + auto& vel = macro_cell.template get<"momentum">(); auto& force = macro_cell.template get<"force">(); auto& part_mask = macro_cell.template get<"particle">(); diff --git a/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix b/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix new file mode 100644 index 0000000..958955f --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/.nix/derivation.nix @@ -0,0 +1,43 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel +, particle_coupling ? "fplbm" +}: + +let +in stdenv.mkDerivation { + pname = "${pname}-examples-stdkr_2d_gpu_${particle_coupling}"; + inherit version; + src = ./..; + + nativeBuildInputs = [ + scons + clang-tools + python3 + ]; + + buildInputs = [ + forstio.core + forstio.async + forstio.codec + forstio.codec-unit + forstio.io + forstio.remote + forstio.remote-filesystem + forstio.codec-json + adaptive-cpp + kel.lbm.core + kel.lbm.sycl + ]; + + preferLocalBuild = true; + + outputs = [ "out" "dev" ]; +} diff --git a/examples/schaefer_turek_durst_krause_rannbacher/SConscript b/examples/schaefer_turek_durst_krause_rannbacher/SConscript new file mode 100644 index 0000000..b33b7ae --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/SConscript @@ -0,0 +1,34 @@ +#!/bin/false + +import os +import os.path +import glob + + +Import('env') + +dir_path = Dir('.').abspath + +# Environment for base library +examples_env = env.Clone(); +examples_env['CXX'] = 'syclcc-clang'; +examples_env['CXXFLAGS'] += ['-O3']; + +examples_env.sources = sorted(glob.glob(dir_path + "/*.cpp")) +examples_env.headers = sorted(glob.glob(dir_path + "/*.hpp")) + +env.sources += examples_env.sources; +env.headers += examples_env.headers; + +# Cavity2D +examples_objects = []; +examples_env.add_source_files(examples_objects, ['sim.cpp'], shared=False); +examples_env.poiseulle_2d_gpu = examples_env.Program('#bin/stdkr_2d', [examples_objects]); + +# Set Alias +env.examples = [ + examples_env.poiseulle_2d_gpu +]; +env.Alias('examples', env.examples); +env.targets += ['examples']; +env.Install('$prefix/bin/', env.examples); diff --git a/examples/schaefer_turek_durst_krause_rannbacher/SConstruct b/examples/schaefer_turek_durst_krause_rannbacher/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/SConstruct @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 + +import sys +import os +import os.path +import glob +import re + + +if sys.version_info < (3,): + def isbasestring(s): + return isinstance(s,basestring) +else: + def isbasestring(s): + return isinstance(s, (str,bytes)) + +def add_kel_source_files(self, sources, filetype, lib_env=None, shared=False, target_post=""): + + if isbasestring(filetype): + dir_path = self.Dir('.').abspath + filetype = sorted(glob.glob(dir_path+"/"+filetype)) + + for path in filetype: + target_name = re.sub( r'(.*?)(\.cpp|\.c\+\+)', r'\1' + target_post, path ) + if shared: + target_name+='.os' + sources.append( self.SharedObject( target=target_name, source=path ) ) + else: + target_name+='.o' + sources.append( self.StaticObject( target=target_name, source=path ) ) + pass + +def isAbsolutePath(key, dirname, env): + assert os.path.isabs(dirname), "%r must have absolute path syntax" % (key,) + +env_vars = Variables( + args=ARGUMENTS +) + +env_vars.Add('prefix', + help='Installation target location of build results and headers', + default='/usr/local/', + validator=isAbsolutePath +) + +env_vars.Add('build_examples', + help='If examples should be built', + default="true" +) + +env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[], + CPPDEFINES=['SAW_UNIX'], + CXXFLAGS=[ + '-std=c++20', + '-g', + '-Wall', + '-Wextra' + ], + LIBS=[ + 'forstio-core', + 'forstio-async', + 'forstio-io' + ] +); +env.__class__.add_source_files = add_kel_source_files +env.Tool('compilation_db'); +env.cdb = env.CompilationDatabase('compile_commands.json'); + +env.objects = []; +env.sources = []; +env.headers = []; +env.targets = []; + +Export('env') +SConscript('SConscript') + +env.Alias('cdb', env.cdb); +env.Alias('all', [env.targets]); +env.Default('all'); + +env.Alias('install', '$prefix') diff --git a/examples/schaefer_turek_durst_krause_rannbacher/common.hpp b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp new file mode 100644 index 0000000..ec6f80f --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp @@ -0,0 +1,70 @@ +#pragma once + +#include <kel/lbm/sycl/lbm.hpp> +#include <kel/lbm/lbm.hpp> +#include <kel/lbm/particle.hpp> +#include <kel/lbm/particle/particle.hpp> + +#include <forstio/io/io.hpp> +#include <forstio/remote/filesystem/easy.hpp> +#include <forstio/codec/json/json.hpp> +#include <forstio/codec/simple.hpp> + +namespace kel { +namespace lbm { + +/** + * length 2.2m + * height 0.41m + * viscosity 10^-3 + * particle diameter 0.1m + * 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 particle_amount = 1ul; + +namespace sch { +using namespace saw::schema; + +using InfoChunk = Chunk<UInt8, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using DfChunk = Chunk<FixedArray<T,Desc::Q>, 1u, dim_x, dim_y>; + +template<typename T, typename Desc> +using ScalarChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using VectorChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using ChunkStruct = Struct< + Member<InfoChunk, "info">, + Member<DfChunk<T,Desc>, "dfs">, + Member<DfChunk<T,Desc>, "dfs_old">, + Member<VectorChunk<T,Desc>, "particle_N">, + Member<ScalarChunk<T,Desc>, "particle_D"> +>; + +template<typename T, typename Desc> +using VelChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>; + +template<typename T> +using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using MacroStruct = Struct< + Member<VelChunk<T,Desc>, "momentum">, + Member<RhoChunk<T>, "density">, + Member<ScalarChunk<T,Desc>, "porosity"> +>; + +template<typename T, typename Desc> +using ParticleSpheroidGroup = ParticleGroup<T,Desc::D,1u,coll::Spheroid<T>>; +} + +} +} diff --git a/examples/schaefer_turek_durst_krause_rannbacher/init.hpp b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp new file mode 100644 index 0000000..b0e4020 --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp @@ -0,0 +1,127 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +template<typename T, typename Desc, typename Coll> +saw::error_or<void> init( + const converter<T>& conv, + saw::data<sch::ChunkStruct<T,Desc>>& fields, + saw::data<sch::MacroStruct<T,Desc>>& macros, + saw::data<sch::ParticleSpheroidGroup<T,Desc>>& particles +){ + (void) conv; + + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + // Set everything as walls + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(1u); + }, + {}, + info_f.get_dims(), + {} + ); + // Fluid + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(2u); + }, + {}, + info_f.get_dims(), + {{1u,1u}} + ); + + // Inflow + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(3u); + }, + {{0u,0u}}, + {{1u,dim_y}}, + {{0u,1u}} + ); + + // Outflow + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(4u); + }, + {{dim_x-1u,0u}}, + {{dim_x, dim_y}}, + {{0u,1u}} + ); + // + auto& df_f = fields.template get<"dfs_old">(); + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"momentum">(); + auto& por_f = macros.template get<"porosity">(); + + iterator<Desc::D>::apply( + [&](auto& index){ + auto& df = df_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; + } + auto eq = equilibrium<T,Desc>(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims(), + {{1u,1u}} + ); + + iterator<Desc::D>::apply( + [&](auto& index){ + saw::data<sch::Vector<T,Desc::D>> middle, ind_vec; + middle.at({{0u}}) = dim_x * 0.25; + middle.at({{1u}}) = dim_y * 0.5; + + ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to<T>(); + ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to<T>(); + + auto dist = middle - ind_vec; + auto dist_2 = saw::math::dot(dist,dist); + if(dist_2.at({}).get() < dim_y*dim_y*0.01){ + porous_f.at(index).at({}) = 0.0; + } + }, + {},// 0-index + df_f.get_dims() + ); + + { + saw::data<sch::Scalar<T>> radius_p; + radius_p.at({}).set(conv.meter_si_to_lbm({{0.16}}).handle().get()); + saw::data<sch::Scalar<T>> dense_p; + dense_p.at({}).set(1); + particles = create_spheroid_particle_group<T,Desc::D,1u>(radius_p, dense_p, {{16u}}); + } + + return saw::make_void(); +} + +} +} diff --git a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp new file mode 100644 index 0000000..6113f6c --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp @@ -0,0 +1,230 @@ +#include "common.hpp" +#include "init.hpp" +#include "step.hpp" + +#include <forstio/codec/args.hpp> + +/** + * For deciding what parameters to use maybe? + */ +namespace args { +using namespace saw::schema; + +using LbmArgsStruct = Struct< + Member<String, "coupling"> +>; + +using LbmArgs = Args< + LbmArgsStruct, + Tuple<> +>; +} + +template<typename T, typename Desc, typename Coll> +saw::error_or<void> lbm_main(const saw::data<args::LbmArgs>& args){ + using namespace kel::lbm; + + auto& an = args.template get<"args">(); + + using dfi = df_info<T,Desc>; + + auto eo_lbm_dir = output_directory(); + if(eo_lbm_dir.is_error()){ + return std::move(eo_lbm_dir.get_error()); + } + auto& lbm_dir = eo_lbm_dir.get_value(); + + auto out_dir = lbm_dir / "stdkr_2d_gpu" / an.template get<"coupling">().stl_view(); + + { + std::error_code ec; + std::filesystem::create_directories(out_dir,ec); + if(ec != std::errc{}){ + return saw::make_error<saw::err::critical>("Could not create output directory"); + } + } + + converter<T> conv { + // delta_x + {{0.41 / dim_y}}, + // delta_t + {{1.067e-4}} + }; + + print_lbm_meta<T,Desc>(conv,{1e-3},{1e-4},{0.5 * dim_y}); + + // 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>>>(); + auto lbm_macro_data_ptr = saw::heap<saw::data<sch::MacroStruct<T,Desc>>>(); + auto lbm_particle_data_ptr = saw::heap<saw::data<sch::ParticleSpheroidGroup<T,Desc>>>(); + + std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl; + + auto eo_aio = saw::setup_async_io(); + if(eo_aio.is_error()){ + return std::move(eo_aio.get_error()); + } + auto& aio = eo_aio.get_value(); + saw::wait_scope wait{aio.event_loop}; + + bool krun = true; + bool print_status = false; + aio.event_port.on_signal(saw::Signal::Terminate).then([&](){ + krun = false; + }).detach(); + aio.event_port.on_signal(saw::Signal::User1).then([&](){ + print_status = true; + }).detach(); + + device dev; + + auto& sycl_q = dev.get_handle(); + + sycl_q.wait(); + { + auto eov = init<T,Desc,Coll>(conv,*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_data_ptr); + if(eov.is_error()){ + 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}; + saw::data<sch::ParticleSpheroidGroup<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_particle_data{sycl_q}; + sycl_q.wait(); + + { + auto eov = dev.copy_to_device(*lbm_data_ptr,lbm_sycl_data); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = dev.copy_to_device(*lbm_macro_data_ptr,lbm_sycl_macro_data); + if(eov.is_error()){ + return eov; + } + } + sycl_q.wait(); + auto lsd_view = make_view(lbm_sycl_data); + auto lsdm_view = make_view(lbm_sycl_macro_data); + auto lsdp_view = make_view(lbm_sycl_particle_data); + + saw::data<sch::UInt64> time_steps{64u*1024ul}; + + auto& info_f = lsd_view.template get<"info">(); + + for(saw::data<sch::UInt64> i{0u}; i < time_steps and krun; ++i){ + // BC + Collision + { + auto eov = step<T,Desc,Coll>(conv,lsd_view,lsdm_view,lsdp_view,i,dev); + if(eov.is_error()){ + return eov; + } + } + + sycl_q.wait(); + if(i.get() % 64u == 0u){ + { + auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_csv_file(out_dir,"stdkr_2d",i.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + } + // Stream + sycl_q.submit([&](acpp::sycl::handler& h){ + component<T,Desc,cmpt::Stream,encode::Sycl<saw::encode::Native>> stream; + + 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]); + } + + auto info = info_f.at(index).get(); + + if(info > 0u){ + stream.apply(lsd_view,index,i); + } + }); + }).wait(); + 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; + print_status = false; + } + print_progress_bar(i.get(), time_steps.get()-1u); + } + // After Loop + sycl_q.wait(); + { + auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_vtk_file(out_dir,"stdkr_2d",time_steps.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + + sycl_q.wait(); + return saw::make_void(); +} + +using FloatT = kel::lbm::sch::Float32; +using DescT = kel::lbm::sch::D2Q9; + +saw::error_or<void> k_main(int argc, char** argv){ + using namespace kel::lbm; + + saw::data<args::LbmArgs> args; + { + auto eov = saw::parse_args(args,argc,argv); + if(eov.is_error()){ + return eov; + } + } + auto& an = args.template get<"args">(); + auto& coupling = an.template get<"coupling">(); + + if(coupling.stl_view() == "hlbm"){ + return lbm_main<FloatT,DescT,method::Hlbm>(args); + } else if(coupling.stl_view() == "fplbm"){ + return lbm_main<FloatT,DescT,method::FpLbm>(args); + } + + return saw::make_error<saw::err::critical>("Invalid coupling"); +} + + +int main(int argc, char** argv){ + auto eov = k_main(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "<<err.get_category(); + auto err_msg = err.get_message(); + if(err_msg.size() > 0u){ + std::cerr<<" - "<<err_msg; + } + std::cerr<<std::endl; + return err.get_id(); + } + return 0; +} diff --git a/examples/schaefer_turek_durst_krause_rannbacher/step.hpp b/examples/schaefer_turek_durst_krause_rannbacher/step.hpp new file mode 100644 index 0000000..1a70f6a --- /dev/null +++ b/examples/schaefer_turek_durst_krause_rannbacher/step.hpp @@ -0,0 +1,100 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +template<typename T, typename Desc, typename Coll> +saw::error_or<void> step( + const converter<T>& conv, + saw::data<sch::Ptr<sch::ChunkStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& fields, + saw::data<sch::Ptr<sch::MacroStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& macros, + saw::data<sch::Ptr<sch::ParticleSpheroidGroup<T,Desc>>,encode::Sycl<saw::encode::Native>>& particles, + saw::data<sch::UInt64> t_i, + device& dev +){ + static constexpr saw::data<sch::UInt64> sub_steps{0u}; + + auto& q = dev.get_handle(); + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + + if constexpr ( std::is_same_v<Coll, method::Hlbm> ){ + q.submit([&](acpp::sycl::handler& h){ + component<T,Desc,cmpt::HlbmReset,encode::Sycl<saw::encode::Native>> hlbm_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); + }); + }).wait(); + q.submit([&](acpp::sycl::handler& h){ + component<T,Desc,cmpt::HlbmOneParticle,encode::Sycl<saw::encode::Native>> hlbm_one_part; + + h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){ + saw::data<sch::FixedArray<sch::UInt64,1u>> index; + for(uint64_t i = 0u; i < 1u; ++i){ + index.at({{i}}).set(idx[i]); + } + + hlbm_one_part.apply(fields,macros,particles,index,t_i,{0u}); + }); + }).wait(); + + 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<true>,encode::Sycl<saw::encode::Native>> flow_in{ + [&](){ + uint64_t target_t_i = 16u; + if(t_i.get() < target_t_i){ + return 1.0 + (0.01 / target_t_i) * t_i.get(); + } + return 1.001; + }() + }; + 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){ + index.at({{i}}).set(idx[i]); + } + + auto info = info_f.at(index); + + switch(info.get()){ + case 0u: + break; + case 1u: + bb.apply(fields,index,t_i); + break; + case 2u: + collision.apply(fields,macros,index,t_i); + break; + case 3u: + 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; + } + }); + }).wait(); + }else if constexpr ( std::is_same_v<Coll, method::FpLbm> ){ + } + + return saw::make_void(); +} + +} +} diff --git a/examples/settling_cubes_2d_ibm_gpu/sim.cpp b/examples/settling_cubes_2d_ibm_gpu/sim.cpp index 33712b5..d7c51bc 100644 --- a/examples/settling_cubes_2d_ibm_gpu/sim.cpp +++ b/examples/settling_cubes_2d_ibm_gpu/sim.cpp @@ -45,7 +45,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<VectorChunk<T,Desc>, "force"> >; @@ -93,7 +93,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template gStarted hearing about similar cases not long after, many of which were successful on the part of the crooks. Scary stuff. et<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); iterator<Desc::D>::apply( [&](auto& index){ @@ -163,7 +163,7 @@ saw::error_or<void> step( auto& parts = p_group.template get<"particles">(); auto& p_mask = p_group.template get<"mask">(); - auto& vels = macros.template get<"velocity">(); + auto& vels = macros.template get<"momentum">(); auto& forces = macros.template get<"force">(); auto p_meta = parts.meta(); @@ -214,7 +214,7 @@ saw::error_or<void> step( // Lagrange indicator position auto p_pos_lag = p_pos + transformed_pos; - // Pick the closest velocity and clamp it + // Pick the closest momentum and clamp it saw::data<sch::FixedArray<sch::UInt64,Desc::D>> p_cell_pos; saw::data<sch::Vector<sch::UInt64,Desc::D>> p_cell_pos_vec; for(uint64_t i{0u}; i < Desc::D; ++i){ diff --git a/examples/settling_spheres_2d_hlbm_gpu/sim.cpp b/examples/settling_spheres_2d_hlbm_gpu/sim.cpp index f7e47b6..5913a1f 100644 --- a/examples/settling_spheres_2d_hlbm_gpu/sim.cpp +++ b/examples/settling_spheres_2d_hlbm_gpu/sim.cpp @@ -49,7 +49,7 @@ using PorChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<PorChunk<T>, "porosity"> >; @@ -91,7 +91,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template gStarted hearing about similar cases not long after, many of which were successful on the part of the crooks. Scary stuff. et<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); iterator<Desc::D>::apply( [&](auto& index){ @@ -168,7 +168,7 @@ saw::error_or<void> step( auto& parts = p_group.template get<"particles">(); auto& p_mask = p_group.template get<"mask">(); - auto& vels = macros.template get<"velocity">(); + auto& vels = macros.template get<"momentum">(); auto& forces = macros.template get<"force">(); auto p_meta = parts.meta(); diff --git a/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp b/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp index 9d2ec12..bbe8cfd 100644 --- a/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp +++ b/examples/stokes_drag_particle_2d_bgk_gpu/sim.cpp @@ -49,7 +49,7 @@ using ForceChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density"> >; @@ -122,7 +122,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); iterator<Desc::D>::apply( [&](auto& index){ diff --git a/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp index 9e96219..9cf974a 100644 --- a/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp +++ b/examples/stokes_drag_particle_2d_fplbm_gpu/sim.cpp @@ -49,7 +49,7 @@ using ForceChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<VelChunk<T,Desc>, "force">, Member<ScalarChunk<T,Desc>, "porosity"> @@ -125,7 +125,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator<Desc::D>::apply( diff --git a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp index 9c28d88..b63a003 100644 --- a/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp +++ b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp @@ -46,7 +46,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<VelChunk<T,Desc>, "force">, Member<ScalarChunk<T,Desc>, "porosity"> @@ -122,7 +122,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator<Desc::D>::apply( diff --git a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp index 75cbe54..6b86c68 100644 --- a/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp +++ b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp @@ -47,7 +47,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density">, Member<VelChunk<T,Desc>, "force">, Member<ScalarChunk<T,Desc>, "porosity"> @@ -124,7 +124,7 @@ saw::error_or<void> setup_initial_conditions( // auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); iterator<Desc::D>::apply( diff --git a/modules/core/c++/boundary.hpp b/modules/core/c++/boundary.hpp index 22a6f34..e7567c9 100644 --- a/modules/core/c++/boundary.hpp +++ b/modules/core/c++/boundary.hpp @@ -133,14 +133,14 @@ template<typename FP, typename Descriptor, typename Encode> class component<FP, Descriptor, cmpt::Equilibrium, Encode> final { private: saw::data<sch::Scalar<FP>> density_; - saw::data<sch::Vector<FP,Descriptor::D>> velocity_; + saw::data<sch::Vector<FP,Descriptor::D>> momentum_; public: component( saw::data<sch::Scalar<FP>> density__, - saw::data<sch::Vector<FP,Descriptor::D>> velocity__ + saw::data<sch::Vector<FP,Descriptor::D>> momentum__ ): density_{density__}, - velocity_{velocity__} + momentum_{momentum__} {} template<typename CellFieldSchema> @@ -151,7 +151,7 @@ public: auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); - auto eq = equilibrium<FP,Descriptor>(density_,velocity_); + auto eq = equilibrium<FP,Descriptor>(density_,momentum_); dfs_old_f.at(index) = eq; } @@ -169,8 +169,10 @@ public: * 0 - 2 - 2 * */ -template<typename FP, typename Descriptor, bool Dir, typename Encode> -class component<FP, Descriptor, cmpt::ZouHeHorizontal<Dir>, Encode> final { +template<typename FP, bool Dir, typename Encode> +class component<FP, sch::Descriptor<2u,9u>, cmpt::ZouHeHorizontal<Dir>, Encode> final { +public: + using Descriptor = sch::Descriptor<2u,9u>; private: saw::data<FP> rho_setting_; public: @@ -179,7 +181,7 @@ public: {} template<typename CellFieldSchema> - void apply(const saw::data<CellFieldSchema, Encode>& field, 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,2u>> index, saw::data<sch::UInt64> time_step) const { using dfi = df_info<FP,Descriptor>; bool is_even = ((time_step.get() % 2) == 0); @@ -207,7 +209,7 @@ public: return {}; }(); - static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); + // static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); if constexpr (Dir) { dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data<FP>{2.0 / 3.0} * rho_vel_x; @@ -221,16 +223,515 @@ public: } }; +template<typename FP, bool East, typename Encode> +class component<FP, sch::Descriptor<3u,27u>, cmpt::ZouHeHorizontal<East>, Encode> final { +public: + using Descriptor = sch::Descriptor<3u,27u>; +private: + saw::data<FP> rho_setting_; +public: + component(const saw::data<FP>& rho_setting__): + rho_setting_{rho_setting__} + {} + + template<typename CellFieldSchema> + 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<FP,Descriptor>; + + bool is_even = ((time_step.get() % 2) == 0); + + auto& info_f = field.template get<"info">(); + + // 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& f = dfs_old_f.at(index); + + /* + * D3Q27 indexing + * + * x=-1 x=0 x=+1 + * + * z=+1: 19 22 25 18 21 24 20 23 26 + * + * z= 0: 1 4 7 0 3 6 2 5 8 + * + * z=-1: 10 13 16 9 12 15 11 14 17 + * + * Boundary normal is x. + * + * East: + * known: cx = +1 + * unknown: cx = -1 + * + * West: + * known: cx = -1 + * unknown: cx = +1 + */ + + + /* + * ============================================================ + * rho * ux + * ============================================================ + * + * rho = S + 2 * sum(unknown) + * + * Therefore: + * + * East: + * rho*ux = rho - S + * + * West: + * rho*ux = S - rho + */ + + saw::data<FP> S; + + if constexpr (East) { + + S = + // cx = 0 + f.at({0u}) + + f.at({3u}) + + f.at({6u}) + + f.at({9u}) + + f.at({12u}) + + f.at({15u}) + + f.at({18u}) + + f.at({21u}) + + f.at({24u}) + + // cx = +1 + + saw::data<FP>{2.0} * ( + f.at({2u}) + + f.at({5u}) + + f.at({8u}) + + f.at({11u}) + + f.at({14u}) + + f.at({17u}) + + f.at({20u}) + + f.at({23u}) + + f.at({26u}) + ); + + } else { + + S = + // cx = 0 + f.at({0u}) + + f.at({3u}) + + f.at({6u}) + + f.at({9u}) + + f.at({12u}) + + f.at({15u}) + + f.at({18u}) + + f.at({21u}) + + f.at({24u}) + + // cx = -1 + + saw::data<FP>{2.0} * ( + f.at({1u}) + + f.at({4u}) + + f.at({7u}) + + f.at({10u}) + + f.at({13u}) + + f.at({16u}) + + f.at({19u}) + + f.at({22u}) + + f.at({25u}) + ); + } + + saw::data<FP> rho_ux; + + if constexpr (East) { + rho_ux = rho_setting_ - S; + } else { + rho_ux = S - rho_setting_; + } + + + /* + * ============================================================ + * Tangential momentum: rho * uy + * ============================================================ + * + * For a pressure boundary, uy and uz are obtained from the + * known populations. + * + * rho*uy: + * + * East: + * + * x=0 contribution + * + 2*x=+1 contribution + * + * West: + * + * x=0 contribution + * + 2*x=-1 contribution + */ + + saw::data<FP> rho_uy; + saw::data<FP> rho_uz; + + if constexpr (East) { + + rho_uy = + // cx = 0 + -f.at({3u}) + + f.at({6u}) + -f.at({12u}) + + f.at({15u}) + -f.at({21u}) + + f.at({24u}) + + // cx = +1 + + saw::data<FP>{2.0} * ( + -f.at({5u}) + + f.at({8u}) + -f.at({14u}) + + f.at({17u}) + -f.at({23u}) + + f.at({26u}) + ); + + rho_uz = + // cx = 0 + -f.at({9u}) + -f.at({12u}) + -f.at({15u}) + + f.at({18u}) + + f.at({21u}) + + f.at({24u}) + + // cx = +1 + + saw::data<FP>{2.0} * ( + -f.at({11u}) + -f.at({14u}) + -f.at({17u}) + +f.at({20u}) + +f.at({23u}) + +f.at({26u}) + ); + + } else { + + rho_uy = + // cx = 0 + -f.at({3u}) + + f.at({6u}) + - f.at({12u}) + + f.at({15u}) + - f.at({21u}) + + f.at({24u}) + + // cx = -1 + + saw::data<FP>{2.0} * ( + -f.at({4u}) + + f.at({7u}) + -f.at({13u}) + + f.at({16u}) + -f.at({22u}) + + f.at({25u}) + ); + + rho_uz = + // cx = 0 + -f.at({9u}) + -f.at({12u}) + -f.at({15u}) + +f.at({18u}) + +f.at({21u}) + +f.at({24u}) + + // cx = -1 + + saw::data<FP>{2.0} * ( + -f.at({10u}) + -f.at({13u}) + -f.at({16u}) + +f.at({19u}) + +f.at({22u}) + +f.at({25u}) + ); + } + + + /* + * ============================================================ + * Velocity + * ============================================================ + */ + + saw::data<FP> ux = rho_ux / rho_setting_; + saw::data<FP> uy = rho_uy / rho_setting_; + saw::data<FP> uz = rho_uz / rho_setting_; + + + /* + * ============================================================ + * Reconstruct unknown distributions + * ============================================================ + * + * Zou/He non-equilibrium bounce-back: + * + * f_i = f_opp + * + (f_eq_i - f_eq_opp) + * + * For D3Q27: + * + * axis: + * + * 6*w*rho*u = 1/6 * rho*u + * + * face diagonal: + * + * 6*w*rho*(c.u) + * = 1/9 * rho*(c.u) + * + * body diagonal: + * + * 6*w*rho*(c.u) + * = 1/36 * rho*(c.u) + * + * The factors below therefore use: + * + * 1/6 + * 1/9 + * 1/36 + * + * ============================================================ + */ + + if constexpr (East) { + + /* + * -------------------------------------------------------- + * (-1, 0, 0) <- (+1, 0, 0) + * -------------------------------------------------------- + */ + + f.at({1u}) = + f.at({2u}) + - saw::data<FP>{1.0 / 6.0} * rho_ux; + + + /* + * -------------------------------------------------------- + * (-1, -1, 0) <- (+1, -1, 0) + * -------------------------------------------------------- + */ + + f.at({4u}) = + f.at({5u}) + - saw::data<FP>{1.0 / 9.0} + * (rho_ux + rho_uy); + + + /* + * -------------------------------------------------------- + * (-1, +1, 0) <- (+1, +1, 0) + * -------------------------------------------------------- + */ + + f.at({7u}) = + f.at({8u}) + + saw::data<FP>{1.0 / 9.0} + * (-rho_ux + rho_uy); + + + /* + * -------------------------------------------------------- + * (-1, 0, -1) <- (+1, 0, -1) + * -------------------------------------------------------- + */ + + f.at({10u}) = + f.at({11u}) + - saw::data<FP>{1.0 / 9.0} + * (rho_ux + rho_uz); + + + /* + * -------------------------------------------------------- + * (-1, -1, -1) <- (+1, -1, -1) + * -------------------------------------------------------- + */ + + f.at({13u}) = + f.at({14u}) + - saw::data<FP>{1.0 / 36.0} + * (rho_ux + rho_uy + rho_uz); + + + /* + * -------------------------------------------------------- + * (-1, +1, -1) <- (+1, +1, -1) + * -------------------------------------------------------- + */ + + f.at({16u}) = + f.at({17u}) + + saw::data<FP>{1.0 / 36.0} + * (-rho_ux + rho_uy + rho_uz); + + + /* + * -------------------------------------------------------- + * (-1, 0, +1) <- (+1, 0, +1) + * -------------------------------------------------------- + */ + + f.at({19u}) = + f.at({20u}) + + saw::data<FP>{1.0 / 9.0} + * (-rho_ux + rho_uz); + + + /* + * -------------------------------------------------------- + * (-1, -1, +1) <- (+1, -1, +1) + * -------------------------------------------------------- + */ + + f.at({22u}) = + f.at({23u}) + + saw::data<FP>{1.0 / 36.0} + * (-rho_ux - rho_uy + rho_uz); + + + /* + * -------------------------------------------------------- + * (-1, +1, +1) <- (+1, +1, +1) + * -------------------------------------------------------- + */ + + f.at({25u}) = + f.at({26u}) + + saw::data<FP>{1.0 / 36.0} + * (-rho_ux + rho_uy + rho_uz); + + } else { + + /* + * -------------------------------------------------------- + * (+1, 0, 0) <- (-1, 0, 0) + * -------------------------------------------------------- + */ + + f.at({2u}) = + f.at({1u}) + + saw::data<FP>{1.0 / 6.0} * rho_ux; + + + /* + * -------------------------------------------------------- + * (+1, -1, 0) <- (-1, -1, 0) + * -------------------------------------------------------- + */ + + f.at({5u}) = + f.at({4u}) + + saw::data<FP>{1.0 / 9.0} + * (rho_ux - rho_uy); + + + /* + * -------------------------------------------------------- + * (+1, +1, 0) <- (-1, +1, 0) + * -------------------------------------------------------- + */ + + f.at({8u}) = + f.at({7u}) + + saw::data<FP>{1.0 / 9.0} + * (rho_ux + rho_uy); + + + /* + * -------------------------------------------------------- + * (+1, 0, -1) <- (-1, 0, -1) + * -------------------------------------------------------- + */ + + f.at({11u}) = + f.at({10u}) + + saw::data<FP>{1.0 / 9.0} + * (rho_ux - rho_uz); + + + /* + * -------------------------------------------------------- + * (+1, -1, -1) <- (-1, -1, -1) + * -------------------------------------------------------- + */ + + f.at({14u}) = + f.at({13u}) + + saw::data<FP>{1.0 / 36.0} + * (rho_ux - rho_uy - rho_uz); + + + /* + * -------------------------------------------------------- + * (+1, +1, -1) <- (-1, +1, -1) + * -------------------------------------------------------- + */ + + f.at({17u}) = + f.at({16u}) + + saw::data<FP>{1.0 / 36.0} + * (rho_ux + rho_uy - rho_uz); + + + /* + * -------------------------------------------------------- + * (+1, 0, +1) <- (-1, 0, +1) + * -------------------------------------------------------- + */ + + f.at({20u}) = + f.at({19u}) + + saw::data<FP>{1.0 / 9.0} + * (rho_ux + rho_uz); + + + /* + * -------------------------------------------------------- + * (+1, -1, +1) <- (-1, -1, +1) + * -------------------------------------------------------- + */ + + f.at({23u}) = + f.at({22u}) + + saw::data<FP>{1.0 / 36.0} + * (rho_ux - rho_uy + rho_uz); + + + /* + * -------------------------------------------------------- + * (+1, +1, +1) <- (-1, +1, +1) + * -------------------------------------------------------- + */ + + f.at({26u}) = + f.at({25u}) + + saw::data<FP>{1.0 / 36.0} + * (rho_ux + rho_uy + rho_uz); + } + } +}; + template<typename T, typename Descriptor, bool East, typename Encode> class component<T, Descriptor, cmpt::ZouHeVelocityX<East>, Encode> final { private: - saw::data<sch::Vector<T,Descriptor::D>> velocity_; + saw::data<sch::Vector<T,Descriptor::D>> momentum_; public: component( - saw::data<sch::Vector<T,Descriptor::D>> velocity__ + saw::data<sch::Vector<T,Descriptor::D>> momentum__ ): - velocity_{velocity__} + momentum_{momentum__} {} template<typename CellFieldSchema> @@ -252,20 +753,20 @@ public: } saw::data<sch::Scalar<T>> rho; - rho.at({}) = (dfs.at({0u}) + dfs.at({4u}) + dfs.at({3u}) + dir_sum * 2) / (velocity_.at({{0u}}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0)}); + rho.at({}) = (dfs.at({0u}) + dfs.at({4u}) + dfs.at({3u}) + dir_sum * 2) / (momentum_.at({{0u}}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0)}); if constexpr (East) { - dfs.at({2u}) = dfs.at({1u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(2.0 / 3.0)} * rho.at({}) * velocity_.at({{0u}}); - dfs.at({6u}) = dfs.at({5u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * velocity_.at({{1u}})) * 0.5f; - dfs.at({8u}) = dfs.at({7u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * velocity_.at({{1u}})) * 0.5f; + dfs.at({2u}) = dfs.at({1u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(2.0 / 3.0)} * rho.at({}) * momentum_.at({{0u}}); + dfs.at({6u}) = dfs.at({5u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * momentum_.at({{1u}})) * 0.5f; + dfs.at({8u}) = dfs.at({7u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * momentum_.at({{1u}})) * 0.5f; }else if constexpr ( not East ){ - dfs.at({1u}) = dfs.at({2u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(2.0 / 3.0)} * rho.at({}) * velocity_.at({{0u}}); - dfs.at({5u}) = dfs.at({6u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * velocity_.at({{1u}}))*0.5f; - dfs.at({7u}) = dfs.at({8u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * velocity_.at({{0u}}) - + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * velocity_.at({{1u}}))*0.5f; + dfs.at({1u}) = dfs.at({2u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(2.0 / 3.0)} * rho.at({}) * momentum_.at({{0u}}); + dfs.at({5u}) = dfs.at({6u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({4u}) - dfs.at({3u}) + rho.at({}) * momentum_.at({{1u}}))*0.5f; + dfs.at({7u}) = dfs.at({8u}) + saw::data<T>{static_cast<saw::native_data_type<T>::type>(1.0 / 6.0)} * rho.at({}) * momentum_.at({{0u}}) + + (dfs.at({3u}) - dfs.at({4u}) + rho.at({}) * momentum_.at({{1u}}))*0.5f; } } }; diff --git a/modules/core/c++/boundary/common.hpp b/modules/core/c++/boundary/common.hpp new file mode 100644 index 0000000..3feb451 --- /dev/null +++ b/modules/core/c++/boundary/common.hpp @@ -0,0 +1,8 @@ +#pragma once + +#include "../common.hpp" + +namespace kel { +namespace lbm { +} +} diff --git a/modules/core/c++/boundary/zou_he_pressure.hpp b/modules/core/c++/boundary/zou_he_pressure.hpp new file mode 100644 index 0000000..7fc8591 --- /dev/null +++ b/modules/core/c++/boundary/zou_he_pressure.hpp @@ -0,0 +1,9 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +} +} diff --git a/modules/core/c++/collision.hpp b/modules/core/c++/collision.hpp index 023f61f..39514b8 100644 --- a/modules/core/c++/collision.hpp +++ b/modules/core/c++/collision.hpp @@ -75,7 +75,7 @@ public: 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">(); + auto& vel_f = macros.template get<"momentum">(); saw::data<sch::Scalar<T>>& rho = rho_f.at(index); saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index); @@ -135,7 +135,7 @@ public: auto& dfs = dfs_old_f.at(index); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data<sch::Scalar<T>>& rho = rho_f.at(index); diff --git a/modules/core/c++/component/psm.hpp b/modules/core/c++/component/psm.hpp index cada53c..5832760 100644 --- a/modules/core/c++/component/psm.hpp +++ b/modules/core/c++/component/psm.hpp @@ -39,7 +39,7 @@ public: auto& porous_f = macros.template get<"porosity">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data<sch::Scalar<T>>& rho = rho_f.at(index); saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index); diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index 37074b0..382d067 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -46,7 +46,7 @@ public: auto& rho_f = macros.template get<"density">(); auto& rho = rho_f.at(index); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& vel = vel_f.at(index); compute_rho_u<T,Descriptor>(dfs,rho,vel); @@ -86,7 +86,7 @@ public: auto& rho_f = macros.template get<"density">(); saw::data<sch::Scalar<T>>& rho = rho_f.at(index); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data<sch::Vector<T,Descriptor::D>> vel = vel_f.at(index); auto& force_f = macros.template get<"force">(); @@ -148,7 +148,7 @@ public: // auto& dfs = dfs_old_f.at(index); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); auto& force_f = macros.template get<"force">(); @@ -204,7 +204,7 @@ public: // vel_s is technically time the density of the particle? - force = ( vel_s * rho - vel * rho ) * two * flip_por; + force = ( vel_s * rho - vel * rho ) /* two */ * flip_por; force_p = force_p - force; }, aabb.template get<"a">(), aabb.template get<"b">()); @@ -245,7 +245,7 @@ public: auto& dfs = dfs_old_f.at(index); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); auto& por_f = macros.template get<"porosity">(); // Temporary find a better way diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp index cc21dcf..e8ea66a 100644 --- a/modules/core/c++/hlbm.hpp +++ b/modules/core/c++/hlbm.hpp @@ -70,7 +70,7 @@ public: auto& porosity_f = macros.template get<"porosity">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data<sch::Scalar<T>>& rho = rho_f.at(index); saw::data<sch::Vector<T,Desc::D>>& vel = vel_f.at(index); @@ -125,7 +125,7 @@ public: auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); auto& part_spheroid_group = part_group; - auto& mvel_f = macros.template get<"velocity">(); + auto& mvel_f = macros.template get<"momentum">(); auto& mrho_f = macros.template get<"density">(); auto& mpor_f = macros.template get<"porosity">(); diff --git a/modules/core/c++/macroscopic.hpp b/modules/core/c++/macroscopic.hpp index 19ab3e9..768fef0 100644 --- a/modules/core/c++/macroscopic.hpp +++ b/modules/core/c++/macroscopic.hpp @@ -28,6 +28,7 @@ void compute_rho_u ( } } + // TODO - Change to momentum, but that needs checks in every impacted function for(size_t i = 0; i < Desc::D; ++i){ vel().at({{i}}) = vel().at({{i}}) / rho().at({}); } diff --git a/modules/core/c++/particle.hpp b/modules/core/c++/particle.hpp index 2c60884..efcc3fb 100644 --- a/modules/core/c++/particle.hpp +++ b/modules/core/c++/particle.hpp @@ -25,7 +25,7 @@ public: // Compute forces - // Update particle velocity + // Update particle momentum verlet_step_lambda_old<T,Descriptor::D>(p,{1.0}); // Update porosity over lattice nodes diff --git a/modules/core/c++/psm.hpp b/modules/core/c++/psm.hpp index 15ad0f8..2daff5e 100644 --- a/modules/core/c++/psm.hpp +++ b/modules/core/c++/psm.hpp @@ -55,7 +55,7 @@ public: auto& porous_f = macros.template get<"porosity">(); auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); + auto& vel_f = macros.template get<"momentum">(); saw::data<sch::Scalar<T>>& rho = rho_f.at(index); saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index); @@ -105,116 +105,418 @@ public: component() = default; 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 { - - using dfi = df_info<T,Descriptor>; - bool is_even = ((time_step.get() % 2) == 0); - - saw::data<sch::Scalar<T>> one; - one.at({}) = 1.0; +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 { - auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + using dfi = df_info<T,Descriptor>; - auto& por_f = macros.template get<"porosity">(); - auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); - auto& force_f = macros.template get<"force">(); - - { - auto parts = pg.template get<"particles">(); - auto parts_size = parts.meta().at({0u}); - - auto& p_coll = pg.template get<"collision">().at({}); - auto& p_rad = p_coll.template get<"radius">(); - - auto& pi = parts.at(index); - auto& pirb = pi.template get<"rigid_body">(); - saw::data<sch::Vector<T,Descriptor::D>>& pirb_pos = pirb.template get<"position">(); - auto& pirb_pos_old = pirb.template get<"position_old">(); - - saw::data<sch::Scalar<T>> ts; - ts.at({}) = one.at({}) / sub_steps.template cast_to<T>(); - - saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> start; - saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> stop; - - auto eo_aabb = particle_aabb<typename ParticleSchema::ValueType>::calculate(pg,{{0u}},vel_f.meta()); - if(eo_aabb.is_error()){ - return; - } - auto& aabb = eo_aabb.get_value(); + bool is_even = ((time_step.get() % 2u) == 0u); - /// Ok, I iterate over the space which covers our particle? So lower bounds to upper bounds - start = aabb.template get<"a">(); - stop = aabb.template get<"b">(); + /* + * ------------------------------------------------------------ + * Fields + * ------------------------------------------------------------ + */ - saw::data<sch::Vector<T,Descriptor::D>> force_p; - for(uint64_t i{0u}; i < Descriptor::D; ++i){ - force_p.at({{i}}) = 0.0; - } - auto vel_p_old = (pirb_pos-pirb_pos_old) / ts; - - iterator<Descriptor::D>::apply([&](const auto& index_f) -> void{ - // ask for the d_k value here. - // For every value im iterating over I need sth - // std::cout<<"Pos: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl; - - auto& dfs = dfs_old_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; - saw::data<sch::Scalar<T>> eps; - eps.at({}) = 1.5f; - - auto& por = por_f.at(index_f); - por = particle_porosity<T,Descriptor::D,1u,por::ParticleSpheroid<T>>::calculate(rel_dist,p_rad,eps); - - if(por.at({}).get() >= 1.0f){ - return; - } - - saw::data<sch::Vector<T,Descriptor::D>> momentum; - for(uint64_t i{0u}; i < Descriptor::D; ++i){ - momentum.at({{i}}) = 0.0; - } - - for(uint64_t i{0u}; i < Descriptor::Q; ++i){ - saw::data<sch::Vector<T,Descriptor::D>> e_i; - saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> n_ind_i; - for(uint64_t k{0u}; k < Descriptor::D; ++k){ - e_i.at({{k}}) = (dfi::directions[i])[k]; - n_ind_i.at({k}) = index_f.at({k}) + (dfi::directions[i])[k]; - } - - uint64_t i_opp = dfi::opposite_index[i]; - - auto u_p_e = saw::math::dot(e_i,vel_p_old); - - saw::data<T> dfs_added = dfs.at({i})*(saw::data<T>{1}-u_p_e.at({})) + dfs_old_f.at(n_ind_i).at({i_opp})*(saw::data<T>{1}+u_p_e.at({})); - saw::data<sch::Scalar<T>> dfs_added_v; - dfs_added_v.at({}) = dfs_added; - auto ei_dfs = e_i * dfs_added_v; - - momentum = momentum + ei_dfs; - } - // technically needs to adjust for rotation as well - - auto& force = force_f.at(index_f); - auto& rho = rho_f.at(index_f); - // To Fluid - - auto flip_por = one - por; - force = momentum * flip_por; - // To Particle - force_p = force_p - force; - },start,stop); - - auto& pirb_acc = pirb.template get<"acceleration">(); - pirb_acc = force_p / pg.template get<"total_mass">().at({}); - - for(saw::data<sch::UInt64> i{0u}; i < sub_steps; ++i){ - verlet_step_lambda<T,Descriptor::D>(pi,ts); - } - } - } + auto& dfs_old_f = + (is_even) + ? field.template get<"dfs_old">() + : field.template get<"dfs">(); + + auto& por_f = macros.template get<"porosity">(); + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"momentum">(); + auto& force_f = macros.template get<"force">(); + + saw::data<sch::Scalar<T>> one; + one.at({}) = 1.0; + + saw::data<sch::Scalar<T>> half; + half.at({}) = 0.5; + + /* + * ------------------------------------------------------------ + * Particle + * ------------------------------------------------------------ + */ + + { + auto parts = pg.template get<"particles">(); + + auto& p_coll = + pg.template get<"collision">().at({}); + + auto& p_rad = + p_coll.template get<"radius">(); + + auto& pi = parts.at(index); + + auto& pirb = + pi.template get<"rigid_body">(); + + saw::data<sch::Vector<T,Descriptor::D>>& pirb_pos = + pirb.template get<"position">(); + + auto& pirb_pos_old = + pirb.template get<"position_old">(); + + + /* + * -------------------------------------------------------- + * Particle timestep + * -------------------------------------------------------- + */ + + saw::data<sch::Scalar<T>> ts; + + ts.at({}) = + one.at({}) + / sub_steps.template cast_to<T>(); + + + /* + * -------------------------------------------------------- + * Particle momentum + * -------------------------------------------------------- + */ + + auto vel_p_old = + (pirb_pos - pirb_pos_old) / ts; + + + /* + * -------------------------------------------------------- + * Find particle AABB + * -------------------------------------------------------- + */ + + auto eo_aabb = + particle_aabb< + typename ParticleSchema::ValueType + >::calculate( + pg, + {{0u}}, + vel_f.meta() + ); + + if(eo_aabb.is_error()){ + return; + } + + auto& aabb = eo_aabb.get_value(); + + saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> start; + saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> stop; + + start = aabb.template get<"a">(); + stop = aabb.template get<"b">(); + + + /* + * -------------------------------------------------------- + * Total force acting on particle + * -------------------------------------------------------- + */ + + saw::data<sch::Vector<T,Descriptor::D>> force_p; + + for(uint64_t k = 0u; k < Descriptor::D; ++k){ + force_p.at({{k}}) = 0.0; + } + + + /* + * -------------------------------------------------------- + * Iterate over particle AABB + * -------------------------------------------------------- + */ + + iterator<Descriptor::D>::apply( + [&](const auto& index_f) -> void { + + /* + * ------------------------------------------------ + * Distribution function at current cell + * ------------------------------------------------ + */ + + auto& dfs = + dfs_old_f.at(index_f); + + + /* + * ------------------------------------------------ + * Particle relative position + * ------------------------------------------------ + */ + + saw::data<sch::Vector<T,Descriptor::D>> rel_dist = + saw::math::vectorize_data(index_f) + .template cast_to<T>() + - pirb_pos; + + + /* + * ------------------------------------------------ + * PSM particle porosity + * ------------------------------------------------ + */ + + saw::data<sch::Scalar<T>> eps; + eps.at({}) = 1.5; + + auto& por = + por_f.at(index_f); + + por = + particle_porosity< + T, + Descriptor::D, + 1u, + por::ParticleSpheroid<T> + >::calculate( + rel_dist, + p_rad, + eps + ); + + + /* + * Outside particle. + */ + + if(por.at({}).get() >= 1.0){ + return; + } + + + /* + * ------------------------------------------------ + * Momentum exchange for this cell + * ------------------------------------------------ + */ + + saw::data<sch::Vector<T,Descriptor::D>> momentum; + + for(uint64_t k = 0u; k < Descriptor::D; ++k){ + momentum.at({{k}}) = 0.0; + } + + + /* + * ------------------------------------------------ + * D3Q27 / general Q momentum exchange + * + * We exclude i=0 because the rest population has + * no momentum. + * + * The original code summed both directions of + * every opposite pair. Therefore we apply 1/2 + * after the sum. + * ------------------------------------------------ + */ + + for(uint64_t i = 1u; i < Descriptor::Q; ++i){ + + uint64_t i_opp = + dfi::opposite_index[i]; + + + /* + * Direction vector e_i. + */ + + saw::data<sch::Vector<T,Descriptor::D>> e_i; + + saw::data< + sch::FixedArray< + sch::UInt64, + Descriptor::D + > + > n_ind_i; + + + for(uint64_t k = 0u; + k < Descriptor::D; + ++k){ + + e_i.at({{k}}) = + dfi::directions[i][k]; + + /* + * Neighbor in direction e_i. + */ + + n_ind_i.at({k}) = + index_f.at({k}) + + dfi::directions[i][k]; + } + + + /* + * ------------------------------------------------ + * Particle momentum projected onto e_i + * ------------------------------------------------ + */ + + auto u_p_e = + saw::math::dot( + e_i, + vel_p_old + ); + + + /* + * ------------------------------------------------ + * Momentum-exchange population + * + * f_i(x) + * + * + f_-i(x + e_i) + * + * with moving-wall correction. + * ------------------------------------------------ + */ + + saw::data<T> dfs_added = + dfs.at({i}) + * ( + saw::data<T>{1.0} + - u_p_e.at({}) + ) + + + dfs_old_f + .at(n_ind_i) + .at({i_opp}) + * ( + saw::data<T>{1.0} + + u_p_e.at({}) + ); + + + /* + * Convert scalar into vector multiplication. + */ + + saw::data<sch::Scalar<T>> dfs_added_v; + + dfs_added_v.at({}) = + dfs_added; + + + /* + * ------------------------------------------------ + * Delta p_i + * + * ------------------------------------------------ + */ + + momentum = + momentum + + e_i * dfs_added_v; + } + + + /* + * ------------------------------------------------ + * Since both i and i_opp were included above, + * compensate for double counting. + * + * This is the important change from your original + * implementation. + * ------------------------------------------------ + */ + + + + momentum = + momentum * half; + + + /* + * ------------------------------------------------ + * PSM solid fraction + * + * por = 1 -> fluid + * por = 0 -> particle + * + * Therefore: + * + * flip_por = 1 - por + * ------------------------------------------------ + */ + + auto flip_por = + one - por; + + + /* + * ------------------------------------------------ + * Force transferred TO FLUID + * ------------------------------------------------ + */ + + auto& force = + force_f.at(index_f); + + force = + momentum * flip_por; + + + /* + * ------------------------------------------------ + * Force transferred TO PARTICLE + * + * Newton's third law. + * ------------------------------------------------ + */ + + force_p = + force_p - force; + + }, + start, + stop + ); + + + /* + * ------------------------------------------------------------ + * Particle acceleration + * ------------------------------------------------------------ + */ + + auto& pirb_acc = + pirb.template get<"acceleration">(); + + pirb_acc = + force_p + / pg.template get<"total_mass">().at({}); + + + /* + * ------------------------------------------------------------ + * Integrate particle + * ------------------------------------------------------------ + */ + + for(saw::data<sch::UInt64> i{0u}; + i < sub_steps; + ++i){ + + verlet_step_lambda<T,Descriptor::D>( + pi, + ts + ); + } + } +} }; } diff --git a/modules/core/c++/write_csv.hpp b/modules/core/c++/write_csv.hpp index a60e208..aa28613 100644 --- a/modules/core/c++/write_csv.hpp +++ b/modules/core/c++/write_csv.hpp @@ -91,7 +91,7 @@ struct lbm_csv_writer<sch::Chunk<T,Ghost,D...>> { template<typename T, uint64_t D> struct lbm_csv_writer<sch::Vector<T,D>> { static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::Vector<T,D>>& field){ - static_assert(D > 0, "Non-dimensionality is bad for velocity."); + static_assert(D > 0, "Non-dimensionality is bad for momentum."); // csv_file<<"VECTORS "<<name<<" float\n"; for(uint64_t i = 0u; i < D; ++i){ diff --git a/modules/core/c++/write_vtk.hpp b/modules/core/c++/write_vtk.hpp index e852172..67e8a2f 100644 --- a/modules/core/c++/write_vtk.hpp +++ b/modules/core/c++/write_vtk.hpp @@ -115,7 +115,7 @@ struct lbm_vtk_writer<sch::Vector<T,D>> { return saw::make_void(); } static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::Vector<T,D>>& field){ - static_assert(D > 0, "Non-dimensionality is bad for velocity."); + static_assert(D > 0, "Non-dimensionality is bad for momentum."); static_assert(D <= 3, "4th dimension as well. Mostly due to vtk."); // vtk_file<<"VECTORS "<<name<<" float\n"; diff --git a/modules/core/tests/equilibrium.cpp b/modules/core/tests/equilibrium.cpp index 20a1f08..7d36600 100644 --- a/modules/core/tests/equilibrium.cpp +++ b/modules/core/tests/equilibrium.cpp @@ -20,7 +20,7 @@ void check_equilibrium(){ auto eq = lbm::equilibrium<lbm::sch::Float64,Descriptor>(rho,vel); for(saw::data<lbm::sch::UInt64> i{0u}; i.get() < Descriptor::Q; ++i){ - SAW_EXPECT(eq.at(i).get() == dfi::weights[i.get()], std::string{"No velocity and normalized rho should be exactly the weights: "} + std::to_string(eq.at(i).get()) + std::string{" "} + std::to_string(dfi::weights[i.get()])); + SAW_EXPECT(eq.at(i).get() == dfi::weights[i.get()], std::string{"No momentum and normalized rho should be exactly the weights: "} + std::to_string(eq.at(i).get()) + std::string{" "} + std::to_string(dfi::weights[i.get()])); } } diff --git a/modules/core/tests/vtk_write.cpp b/modules/core/tests/vtk_write.cpp index c99d752..1de25bf 100644 --- a/modules/core/tests/vtk_write.cpp +++ b/modules/core/tests/vtk_write.cpp @@ -31,7 +31,7 @@ using CellStruct = Struct< template<typename T, uint64_t D> using MacroStruct = Struct< - Member<FixedArray<T,D>, "velocity">, + Member<FixedArray<T,D>, "momentum">, Member<T, "pressure"> >; @@ -47,7 +47,7 @@ SAW_TEST("VTK Write test example"){ saw::data<sch::Array<sch::MacroStruct<sch::T,2>, 2>> cells{{{2u,2u}}}; auto& cell_0 = cells.at({{{0,0}}}); - cell_0.template get<"velocity">()= {{0.5,-0.1}}; + cell_0.template get<"momentum">()= {{0.5,-0.1}}; cell_0.template get<"pressure">().set(1.1); auto eov = lbm::impl::lbm_vtk_writer<sch::Array<sch::MacroStruct<sch::T,2>, 2>>::apply(sstream, cells); @@ -56,7 +56,7 @@ SAW_TEST("VTK Write test example"){ // I want to print it to see it for myself. For now I have no tooling to more easily view associated and potentially generated files std::cout<<sstream.str()<<std::endl; - static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS velocity float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n"; + static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS momentum float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n"; SAW_EXPECT(sstream.str() == comparison_str, "Expected different encoding"); // using Type = typename parameter_pack_type<i,T...>::type; @@ -74,7 +74,7 @@ SAW_TEST("VTK Write raw test example"){ auto meta = cells.dims(); auto& cell_0 = cells.at({{{0,0}}}); - cell_0.template get<"velocity">()= {{0.5,-0.1}}; + cell_0.template get<"momentum">()= {{0.5,-0.1}}; cell_0.template get<"pressure">().set(1.1); auto eov = lbm::impl::lbm_vtk_writer_raw<sch::MacroStruct<sch::T,2u>,2u>::apply(sstream, &cell_0, meta); @@ -83,7 +83,7 @@ SAW_TEST("VTK Write raw test example"){ // I want to print it to see it for myself. For now I have no tooling to more easily view associated and potentially generated files std::cout<<sstream.str()<<std::endl; - static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS velocity float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n"; + static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS momentum float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n"; SAW_EXPECT(sstream.str() == comparison_str, "Expected different encoding"); // using Type = typename parameter_pack_type<i,T...>::type; diff --git a/poc/triangle_stream/common.hpp b/poc/triangle_stream/common.hpp index 927f8f8..848d337 100644 --- a/poc/triangle_stream/common.hpp +++ b/poc/triangle_stream/common.hpp @@ -37,7 +37,7 @@ using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; template<typename T, typename Desc> using MacroStruct = Struct< - Member<VelChunk<T,Desc>, "velocity">, + Member<VelChunk<T,Desc>, "momentum">, Member<RhoChunk<T>, "density"> >; } |