summaryrefslogtreecommitdiff
path: root/modules/core/c++
diff options
context:
space:
mode:
Diffstat (limited to 'modules/core/c++')
-rw-r--r--modules/core/c++/fplbm.hpp100
-rw-r--r--modules/core/c++/hlbm.hpp38
-rw-r--r--modules/core/c++/lbm.hpp1
3 files changed, 9 insertions, 130 deletions
diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp
index 1f38f51..37074b0 100644
--- a/modules/core/c++/fplbm.hpp
+++ b/modules/core/c++/fplbm.hpp
@@ -80,20 +80,25 @@ public:
auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
auto& dfs = dfs_old_f.at(index);
+ saw::data<sch::Scalar<T>> half;
+ half.at({}) = 0.5f;
+
auto& rho_f = macros.template get<"density">();
saw::data<sch::Scalar<T>>& rho = rho_f.at(index);
auto& vel_f = macros.template get<"velocity">();
- saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);
+ saw::data<sch::Vector<T,Descriptor::D>> vel = vel_f.at(index);
+
+ auto& force_f = macros.template get<"force">();
+ auto& force = force_f.at(index);
+
+ vel = vel + force * half;
// compute_rho_u<T,Descriptor>(dfs,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& por_f = macros.template get<"porosity">();
auto& por = por_f.at(index);
@@ -214,93 +219,6 @@ public:
};
template<typename T, typename Descriptor, typename Encode>
-class component<T, Descriptor, cmpt::FpLbmOneParticleImplicit, Encode> final {
-public:
- using Component = cmpt::FpLbmOneParticle;
-private:
-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 {
- // 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">();
- auto& force_f = macros.template get<"force">();
-
- /**
- * The other methods all work with a flipped porosity.
- * For compat reasons this is also flipped
- */
- // TODO - Change to tuple later
- auto parts = pg.template get<"particles">();
-
- auto parts_size = parts.meta().at({0u});
- auto& pi = parts.at(index);
-
- auto& pirb = pi.template get<"rigid_body">();
- auto& pirb_pos = pirb.template get<"position">();
- auto& pirb_pos_old = pirb.template get<"position_old">();
-
-
- auto& p_coll = pg.template get<"collision">().at({});
- auto& p_rad = p_coll.template get<"radius">();
-
- 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::Scalar<T>> two;
- two.at({}).set(2.0f);
-
- saw::data<sch::Scalar<T>> one;
- one.at({}) = 1.0;
-
- saw::data<sch::Scalar<T>> eps;
- eps.at({}) = 1.5f;
-
- saw::data<sch::Scalar<T>> sss;
- sss.at({}) = one.at({}) / sub_steps.template cast_to<T>();
- auto vel_s = (pirb_pos-pirb_pos_old) / sss;
-
- saw::data<sch::Vector<T,Descriptor::D>> force_p{};
-
- iterator<Descriptor::D>::apply([&](const auto& index_f) -> void {
- auto& force = force_f.at(index_f);
- auto& vel = vel_f.at(index_f);
- auto& por = por_f.at(index_f);
- auto& rho = rho_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;
-
- por = particle_porosity<T,Descriptor::D,1u,por::ParticleSpheroid<T>>::calculate(rel_dist,p_rad,eps);
- auto flip_por = one - por;
-
- // vel_s is technically time the density of the particle?
-
- force = ( vel_s * rho - vel * rho ) * two * flip_por / (one + flip_por);
-
- force_p = force_p - force;
- }, aabb.template get<"a">(), aabb.template get<"b">());
-
- 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,sss);
- }
- }
-};
-
-template<typename T, typename Descriptor, typename Encode>
class component<T, Descriptor, cmpt::FpLbmOneParticleNoVelocity, Encode> final {
private:
saw::data<T> relaxation_;
diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp
index fc19728..cc21dcf 100644
--- a/modules/core/c++/hlbm.hpp
+++ b/modules/core/c++/hlbm.hpp
@@ -234,43 +234,5 @@ public:
}
}
};
-
-template<typename T, typename Desc, typename Encode>
-class component<T, Desc, cmpt::HlbmOneParticleMomentumExchange, Encode> final {
-public:
- 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,Desc::D>> index, saw::data<sch::UInt64> time_step) const {
- //
- using dfi = df_info<T,Desc>;
- 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);
-
- 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];
- }
-
- uint64_t i_opp = dfi::opposite_index[i];
-
- saw::data<T> dfs_added = dfs.at({i}) - dfs_old_f.at(n_ind_i).at({i_opp});
- 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;
- }
- auto& force_f = macros.template get<"force">();
-
- // Set Force
- force_f.at(index) = momentum * macros.template get<"porosity">().at(index);
- }
-};
}
}
diff --git a/modules/core/c++/lbm.hpp b/modules/core/c++/lbm.hpp
index 23d30b1..4e73537 100644
--- a/modules/core/c++/lbm.hpp
+++ b/modules/core/c++/lbm.hpp
@@ -1,6 +1,5 @@
#pragma once
-#include "args.hpp"
#include "schema.hpp"
#include "descriptor.hpp"
#include "boundary.hpp"
generated by cgit v1.2.3 (git 2.25.1) at 2026-08-11 13:10:05 +0000