summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--default.nix5
-rw-r--r--examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp9
-rw-r--r--lib/core/c++/collision.hpp9
-rw-r--r--lib/core/c++/fplbm.hpp177
-rw-r--r--lib/core/c++/lbm.hpp5
5 files changed, 198 insertions, 7 deletions
diff --git a/default.nix b/default.nix
index ef68cab..8b9c2ef 100644
--- a/default.nix
+++ b/default.nix
@@ -167,6 +167,11 @@ in rec {
inherit kel;
};
+ poiseulle_particles_2d_fplbm_gpu = pkgs.callPackage ./examples/poiseulle_particles_2d_fplbm_gpu/.nix/derivation.nix {
+ inherit pname version stdenv forstio adaptive-cpp;
+ inherit kel;
+ };
+
poiseulle_3d = pkgs.callPackage ./examples/poiseulle_3d/.nix/derivation.nix {
inherit pname version stdenv forstio adaptive-cpp;
inherit kel;
diff --git a/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp
index 378154c..93d4b46 100644
--- a/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_fplbm_gpu/sim.cpp
@@ -45,6 +45,9 @@ template<typename T>
using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>;
template<typename T, typename Desc>
+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<RhoChunk<T>, "density">,
@@ -112,7 +115,7 @@ saw::error_or<void> setup_initial_conditions(
[&](auto& index){
auto& df = df_f.at(index);
auto& rho = rho_f.at(index);
- por_f.at(index).at({}) = {1};
+ por_f.at(index).at({}) = {0};
rho.at({}) = {1};
auto& vel = vel_f.at(index);
auto eq = equilibrium<T,Desc>(rho,vel);
@@ -153,7 +156,7 @@ saw::error_or<void> setup_initial_conditions(
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;
+ porous_f.at(index).at({}) = 1.0;
}
},
{},// 0-index
@@ -176,7 +179,7 @@ saw::error_or<void> step(
// auto coll_ev =
q.submit([&](acpp::sycl::handler& h){
- component<T,Desc,cmpt::HLBM,encode::Sycl<saw::encode::Native>> collision{0.65};
+ component<T,Desc,cmpt::FpLbmOneParticleNoVelocity,encode::Sycl<saw::encode::Native>> collision{0.65};
component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb;
component<T,Desc,cmpt::AntiBounceBack<0u>,encode::Sycl<saw::encode::Native>> abb;
diff --git a/lib/core/c++/collision.hpp b/lib/core/c++/collision.hpp
index cd941bd..9c76c1a 100644
--- a/lib/core/c++/collision.hpp
+++ b/lib/core/c++/collision.hpp
@@ -143,6 +143,7 @@ public:
auto& force = force_f.at(index);
auto total_force = force + global_force_;
+
saw::data<sch::Scalar<T>> half;
half.at({}).set(0.5);
saw::data<sch::Vector<T,Descriptor::D>> vel = vel_f.at(index) + total_force * ( half / rho );
@@ -168,7 +169,7 @@ public:
auto ci_dot_u = saw::math::dot(ci,vel);
// saw::data<sch::Vector<T,Descriptor::D>> F_i;
- // F_i = f * (c_i - u * ics2 + <c_i,u> * c_i * ics2 * ics2) * w_i;
+ // F_i = f * ((c_i - u) * ics2 + <c_i,u> * c_i * ics2 * ics2) * w_i;
saw::data<sch::Scalar<T>> w;
w.at({}).set(dfi::weights[i]);
@@ -187,7 +188,11 @@ public:
auto F_i = w * force_projection;
- dfs.at({i}) = dfs.at({i}) + frequency_ * (eq.at(i) - dfs.at({i}) ) + F_i.at({}) * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_);
+ dfs.at({i}) = dfs.at({i})
+ + frequency_
+ * (eq.at(i) - dfs.at({i}) )
+ + F_i.at({})
+ * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_);
}
}
};
diff --git a/lib/core/c++/fplbm.hpp b/lib/core/c++/fplbm.hpp
new file mode 100644
index 0000000..61def6b
--- /dev/null
+++ b/lib/core/c++/fplbm.hpp
@@ -0,0 +1,177 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+namespace cmpt {
+struct FpLbm {};
+struct FpLbmOneParticleNoVelocity{};
+}
+
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::FpLbm, Encode> final {
+private:
+ saw::data<T> relaxation_;
+ saw::data<T> frequency_;
+public:
+ component(typename saw::native_data_type<T>::type relaxation__):
+ relaxation_{{relaxation__}},
+ frequency_{saw::data<T>{1} / relaxation_}
+ {}
+
+ component(const saw::data<T>& relaxation__):
+ relaxation_{relaxation__},
+ frequency_{saw::data<T>{1} / relaxation_}
+ {}
+
+ using Component = cmpt::FpLbm;
+
+
+ template<typename CellFieldSchema, typename MacroFieldSchema>
+ void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+ // void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, saw::data<sch::UInt64> time_step){
+ bool is_even = ((time_step.get() % 2) == 0);
+
+ auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+ auto& dfs = dfs_old_f.at(index);
+
+ auto& rho_f = macros.template get<"density">();
+ auto& vel_f = macros.template get<"velocity">();
+
+ saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
+ 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);
+
+ using dfi = df_info<T,Descriptor>;
+
+ auto& force_f = macros.template get<"force">();
+ auto& force = force_f.at(index);
+
+ auto& porosity_f = macros.template get<"porosity">();
+ auto& porosity = porosity_f.at(index);
+
+ saw::data<sch::Scalar<T>> dfi_inv_cs2;
+ dfi_inv_cs2.at({}).set(dfi::inv_cs2);
+
+ for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+ // saw::data<T> ci_min_u{0};
+ saw::data<sch::Vector<T,Descriptor::D>> ci;
+ for(uint64_t d = 0u; d < Descriptor::D; ++d){
+ ci.at({{d}}).set(static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d]));
+ }
+ auto ci_dot_u = saw::math::dot(ci,vel);
+
+ // saw::data<sch::Vector<T,Descriptor::D>> F_i;
+ // F_i = f * (c_i - u * ics2 + <c_i,u> * c_i * ics2 * ics2) * w_i;
+ saw::data<sch::Scalar<T>> w;
+ w.at({}).set(dfi::weights[i]);
+
+ auto F_i_d = saw::math::dot(force * w, (ci - vel * dfi_inv_cs2 + ci * ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2 ));
+ /*
+ saw::data<sch::Scalar<T>> F_i_sum;
+ for(uint64_t d = 0u; d < Descriptor::D; ++d){
+ saw::data<sch::Scalar<T>> F_i_d;
+ F_i_d.at({}) = F_i.at({{d}});
+ F_i_sum = F_i_sum + F_i_d;
+ }
+ */
+
+ dfs.at({i}) = dfs.at({i}) + frequency_ * (eq.at(i) - dfs.at({i}) ) + F_i_d.at({}) * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_);
+ }
+ }
+};
+
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::FpLbmOneParticleNoVelocity, Encode> final {
+private:
+ saw::data<T> relaxation_;
+ saw::data<T> frequency_;
+public:
+ component(typename saw::native_data_type<T>::type relaxation__):
+ relaxation_{{relaxation__}},
+ frequency_{saw::data<T>{1} / relaxation_}
+ {}
+
+ component(const saw::data<T>& relaxation__):
+ relaxation_{relaxation__},
+ frequency_{saw::data<T>{1} / relaxation_}
+ {}
+
+ using Component = cmpt::FpLbmOneParticleNoVelocity;
+
+ template<typename CellFieldSchema, typename MacroFieldSchema>
+ void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
+ // void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, saw::data<sch::UInt64> time_step){
+ bool is_even = ((time_step.get() % 2) == 0);
+
+ auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
+ auto& dfs = dfs_old_f.at(index);
+
+ auto& rho_f = macros.template get<"density">();
+ auto& vel_f = macros.template get<"velocity">();
+ auto& por_f = macros.template get<"porosity">();
+
+ saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
+
+ saw::data<sch::Scalar<T>> half;
+ half.at({}).set(0.5);
+ saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);// + total_force * ( half / rho );
+
+ compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel);
+ auto eq = equilibrium<T,Descriptor>(rho,vel);
+
+ using dfi = df_info<T,Descriptor>;
+
+ saw::data<sch::Scalar<T>> min_two;
+ min_two.at({}).set(-2);
+
+ auto force = vel * rho * min_two * por_f.at(index);
+
+ saw::data<sch::Scalar<T>> dfi_inv_cs2;
+ dfi_inv_cs2.at({}).set(dfi::inv_cs2);
+
+
+ // auto vel = vel_f.at(index);
+
+ for(uint64_t i = 0u; i < Descriptor::Q; ++i){
+ // saw::data<T> ci_min_u{0};
+ saw::data<sch::Vector<T,Descriptor::D>> ci;
+ for(uint64_t d = 0u; d < Descriptor::D; ++d){
+ ci.at({{d}}).set(static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d]));
+ }
+ auto ci_dot_u = saw::math::dot(ci,vel);
+
+ // saw::data<sch::Vector<T,Descriptor::D>> F_i;
+ // F_i = f * ((c_i - u) * ics2 + <c_i,u> * c_i * ics2 * ics2) * w_i;
+ saw::data<sch::Scalar<T>> w;
+ w.at({}).set(dfi::weights[i]);
+
+ /*
+ saw::data<sch::Scalar<T>> F_i_sum;
+ for(uint64_t d = 0u; d < Descriptor::D; ++d){
+ saw::data<sch::Scalar<T>> F_i_d;
+ F_i_d.at({}) = F_i.at({{d}});
+ F_i_sum = F_i_sum + F_i_d;
+ }
+ */
+ auto term1 = (ci-vel) * dfi_inv_cs2;
+ auto term2 = ci * (ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2);
+
+ auto force_projection = saw::math::dot(term1 + term2, force);
+
+ auto F_i = w * force_projection;
+
+ dfs.at({i}) = dfs.at({i})
+ + frequency_
+ * (eq.at(i) - dfs.at({i}) )
+ + F_i.at({})
+ * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_);
+ }
+ }
+
+};
+}
+}
diff --git a/lib/core/c++/lbm.hpp b/lib/core/c++/lbm.hpp
index 29de83e..fbad908 100644
--- a/lib/core/c++/lbm.hpp
+++ b/lib/core/c++/lbm.hpp
@@ -2,8 +2,6 @@
#include "args.hpp"
#include "schema.hpp"
-#include "flatten.hpp"
-#include "chunk.hpp"
#include "descriptor.hpp"
#include "boundary.hpp"
#include "converter.hpp"
@@ -12,6 +10,9 @@
#include "component.hpp"
#include "environment.hpp"
#include "equilibrium.hpp"
+#include "flatten.hpp"
+#include "fplbm.hpp"
+#include "chunk.hpp"
#include "iterator.hpp"
#include "hlbm.hpp"
#include "macroscopic.hpp"
generated by cgit v1.2.3 (git 2.25.1) at 2026-05-11 12:10:47 +0000