summaryrefslogtreecommitdiff
path: root/modules
diff options
context:
space:
mode:
Diffstat (limited to 'modules')
-rw-r--r--modules/core/c++/descriptor/common.hpp9
-rw-r--r--modules/core/c++/descriptor/d1q3.hpp9
-rw-r--r--modules/core/c++/hlbm.hpp38
-rw-r--r--modules/core/c++/particle.hpp2
-rw-r--r--modules/core/c++/particle/particle.hpp43
-rw-r--r--modules/core/tests/schema.cpp19
6 files changed, 68 insertions, 52 deletions
diff --git a/modules/core/c++/descriptor/common.hpp b/modules/core/c++/descriptor/common.hpp
new file mode 100644
index 0000000..f821510
--- /dev/null
+++ b/modules/core/c++/descriptor/common.hpp
@@ -0,0 +1,9 @@
+#pragma once
+
+#include "../common.hpp"
+
+namespace kel {
+namespace lbm {
+
+}
+}
diff --git a/modules/core/c++/descriptor/d1q3.hpp b/modules/core/c++/descriptor/d1q3.hpp
new file mode 100644
index 0000000..7fc8591
--- /dev/null
+++ b/modules/core/c++/descriptor/d1q3.hpp
@@ -0,0 +1,9 @@
+#pragma once
+
+#include "common.hpp"
+
+namespace kel {
+namespace lbm {
+
+}
+}
diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp
index 25b6bf0..799d2b5 100644
--- a/modules/core/c++/hlbm.hpp
+++ b/modules/core/c++/hlbm.hpp
@@ -79,6 +79,7 @@ public:
auto& N = particle_N_f.at(index);
auto& D = particle_D_f.at(index);
+
// Convex combination of velocities
vel = vel * porosity + [&]() -> saw::data<sch::Vector<T,Desc::D>> {
return (D.at({}).get() > 0.0 ? N * flip_porosity / D : N);
@@ -113,13 +114,20 @@ public:
using dfi = df_info<T,Desc>;
bool is_even = ((time_step.get() % 2) == 0);
+ saw::data<sch::Scalar<T>> one;
+ one.at({}) = 1.0;
+
auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
auto& part_spheroid_group = part_group;
- auto& mvel = macros.template get<"velocity">();
+ auto& mvel_f = macros.template get<"velocity">();
auto& mpor_f = macros.template get<"porosity">();
+
+ auto& particle_N_f = field.template get<"particle_N">();
+ auto& particle_D_f = field.template get<"particle_D">();
+
{
- auto& parts = part_spheroid_group.template get<"particles">();
+ auto parts = part_spheroid_group.template get<"particles">();
auto parts_size = parts.meta().at({0u});
auto& p_coll = part_spheroid_group.template get<"collision">().at({});
@@ -132,7 +140,7 @@ public:
saw::data<sch::FixedArray<sch::UInt64,Desc::D>> start;
saw::data<sch::FixedArray<sch::UInt64,Desc::D>> stop;
- auto eo_aabb = particle_aabb<typename ParticleSchema::ValueType>::calculate(part_spheroid_group,{{0u}},mvel.meta());
+ auto eo_aabb = particle_aabb<typename ParticleSchema::ValueType>::calculate(part_spheroid_group,{{0u}},mvel_f.meta());
if(eo_aabb.is_error()){
return;
}
@@ -143,21 +151,22 @@ public:
stop = aabb.template get<"b">();
saw::data<sch::Vector<T,Desc::D>> force_p;
-
- iterator<Desc::D>::apply([&](const auto& index){
+
+ iterator<Desc::D>::apply([&](const auto& index_f){
// 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);
+ auto& dfs = dfs_old_f.at(index_f);
saw::data<sch::Vector<T,Desc::D>> momentum;
+
for(uint64_t i = 0u; i < Desc::Q; ++i){
saw::data<sch::Vector<T,Desc::D>> e_i;
saw::data<sch::FixedArray<sch::UInt64,Desc::D>> n_ind_i;
for(uint64_t k{0u}; k < Desc::D; ++k){
e_i.at({{k}}) = (dfi::directions[i])[k];
- n_ind_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];
@@ -170,19 +179,26 @@ public:
momentum = momentum + ei_dfs;
}
- auto& mpor = mpor_f.at(index);
- auto rel_dist = saw::math::vectorize_data(index).template cast_to<T>() - pirb_pos;
+ auto& mpor = mpor_f.at(index_f);
+ auto rel_dist = saw::math::vectorize_data(index_f).template cast_to<T>() - pirb_pos;
saw::data<sch::Scalar<T>> eps;
eps.at({}) = 1.5f;
mpor = particle_porosity<T,Desc::D,1u,por::ParticleSpheroid<T>>::calculate(rel_dist,p_rad,eps);
force_p = force_p + momentum * mpor;
+ auto& N = particle_N_f.at(index_f);
+ auto& D = particle_D_f.at(index_f);
+ auto flip_porosity = one - mpor;
+ D = D + flip_porosity;
+ N = N + mvel_f.at(index_f) * flip_porosity;
},start,stop);
-
+
auto& pirb_acc = pirb.template get<"acceleration">();
pirb_acc = force_p;
- verlet_step_lambda(part_spheroid_group,{{0u}},{{1u}});
+ saw::data<sch::Scalar<T>> ts;
+ ts.at({}) = 1u;
+ verlet_step_lambda<T,Desc::D>(pi,ts);
// Check
}
diff --git a/modules/core/c++/particle.hpp b/modules/core/c++/particle.hpp
index 691a74b..2c60884 100644
--- a/modules/core/c++/particle.hpp
+++ b/modules/core/c++/particle.hpp
@@ -26,7 +26,7 @@ public:
// Compute forces
// Update particle velocity
- verlet_step_lambda<T,Descriptor::D>(p,{1.0});
+ verlet_step_lambda_old<T,Descriptor::D>(p,{1.0});
// Update porosity over lattice nodes
diff --git a/modules/core/c++/particle/particle.hpp b/modules/core/c++/particle/particle.hpp
index 6d57930..f759c8d 100644
--- a/modules/core/c++/particle/particle.hpp
+++ b/modules/core/c++/particle/particle.hpp
@@ -137,7 +137,7 @@ saw::data<sch::Particle<T,D, sch::ParticleCollisionSpheroid<T>>> create_spheroid
*/
template<typename T,uint64_t D>
-constexpr auto verlet_step_lambda_old = [](saw::data<sch::Particle<T,D>>& particle, saw::data<sch::Scalar<T>> time_step_delta){
+constexpr auto verlet_step_lambda = [](saw::data<sch::Particle<T,D>>& particle, saw::data<sch::Scalar<T>> time_step_delta){
auto& body = particle.template get<"rigid_body">();
auto& pos = body.template get<"position">();
@@ -158,51 +158,14 @@ constexpr auto verlet_step_lambda_old = [](saw::data<sch::Particle<T,D>>& partic
two.at({}).set(2.0);
pos_new = pos * two - pos_old + pos_acc * tsd_squared;
- // Angular
- saw::data<typename sch::impl::rotation_type_helper<T,D>::Schema> rot_new;
- rot_new = rot * two - rot_old + rot_acc * tsd_squared;
-
- // Swap - Could be std::swap?
+ // Swap - Could be std::swap? No, this is device code
pos_old = pos;
pos = pos_new;
- rot_old = rot;
- rot = rot_new;
-};
-
-template<typename T,uint64_t D,uint64_t PC, typename Coll>
-constexpr auto verlet_step_lambda = [](saw::data<sch::ParticleGroup<T,D,PC,Coll>>& pgrp, saw::data<sch::FixedArray<T,1u>> index, saw::data<sch::Scalar<T>> time_step_delta){
- auto& parts = pgrp.template get<"particles">();
- auto& particle = parts.at(index);
-
- auto& body = particle.template get<"rigid_body">();
-
- auto& pos = body.template get<"position">();
- auto& pos_old = body.template get<"position_old">();
-
- auto& pos_acc = body.template get<"acceleration">();
-
- auto& rot = body.template get<"rotation">();
- auto& rot_old = body.template get<"rotation_old">();
-
- auto& rot_acc = body.template get<"angular_acceleration">();
-
- auto tsd_squared = time_step_delta * time_step_delta;
-
- saw::data<sch::Vector<T,D>> pos_new;
- // Actual step
- saw::data<sch::Scalar<T>> two;
- two.at({}).set(2.0);
- pos_new = pos * two - pos_old + pos_acc * tsd_squared;
-
// Angular
saw::data<typename sch::impl::rotation_type_helper<T,D>::Schema> rot_new;
rot_new = rot * two - rot_old + rot_acc * tsd_squared;
-
- // Swap - Could be std::swap?
- pos_old = pos;
- pos = pos_new;
-
+
rot_old = rot;
rot = rot_new;
};
diff --git a/modules/core/tests/schema.cpp b/modules/core/tests/schema.cpp
new file mode 100644
index 0000000..b705af0
--- /dev/null
+++ b/modules/core/tests/schema.cpp
@@ -0,0 +1,19 @@
+#include <forstio/test/suite.hpp>
+
+#include "../c++/abstract/schema.hpp"
+
+namespace {
+namespace ts {
+using namespace kel::lbm::sch;
+
+using TPtr = Ptr<UInt16>;
+using TVal = UInt16;
+}
+
+SAW_TEST("Abstract/Schema/Ptr"){
+ using namespace kel::lbm;
+
+ // static_assert(std::is_same_v<typename schema_decay<ts::TPtr>::Schema, ts::TVal>, "schema decay failed");
+}
+
+}
generated by cgit v1.2.3 (git 2.25.1) at 2026-07-21 22:46:00 +0000