summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--c++/collision.hpp12
-rw-r--r--c++/component.hpp5
-rw-r--r--c++/descriptor.hpp7
-rw-r--r--c++/equilibrium.hpp16
-rw-r--r--c++/examples/cavity_2d.cpp68
-rw-r--r--tests/equilibrium.cpp20
6 files changed, 72 insertions, 56 deletions
diff --git a/c++/collision.hpp b/c++/collision.hpp
index 4e16832..87187c2 100644
--- a/c++/collision.hpp
+++ b/c++/collision.hpp
@@ -9,16 +9,22 @@ namespace cmpt {
struct BGK {};
}
-template<typename T>
-class component<T, cmpt::BGK> {
+template<typename T, typename Descriptor>
+class component<T, Descriptor, cmpt::BGK> {
public:
+ using Component = cmpt::BGK;
+
using access = cmpt::access_tuple<
- cmpt::access<"dfs", 0, true, true>
+ cmpt::access<"dfs", 1, true, true, true>
>;
static constexpr saw::string_literal name = "collision";
static constexpr saw::string_literal after = "";
static constexpr saw::string_literal before = "streaming";
+
+ void apply(saw::data<sch::CellField>& field){
+
+ }
};
}
}
diff --git a/c++/component.hpp b/c++/component.hpp
index 14f908a..c67387b 100644
--- a/c++/component.hpp
+++ b/c++/component.hpp
@@ -6,13 +6,14 @@ namespace kel {
namespace lbm {
namespace cmpt {
-template<saw::string_literal Name, uint64_t Dist, bool Read, bool Write>
+template<saw::string_literal Name, uint64_t Dist, bool Read, bool Write, bool SkipSync = false>
class access {
public:
static constexpr saw::string_literal name = Name;
static constexpr uint64_t distance = Dist;
static constexpr bool read = Read;
static constexpr bool write = Write;
+ static constexpr bool skip_sync = SkipSync;
};
template<typename... Acc>
@@ -24,7 +25,7 @@ public:
/**
* Compponent class which forms the basis of the
*/
-template<typename T>
+template<typename T, typename Descriptor, typename Cmpt>
class component {};
}
}
diff --git a/c++/descriptor.hpp b/c++/descriptor.hpp
index 1228e10..014327c 100644
--- a/c++/descriptor.hpp
+++ b/c++/descriptor.hpp
@@ -196,12 +196,11 @@ private:
public:
data() = default;
data(const data<MetaSchema,Encode>& inner_meta__):
- inner_{inner_meta__},
- inner_meta_{inner_meta__}
+ inner_{inner_meta__}
{}
- const data<MetaSchema, Encode>& meta() const {
- return inner_.dims();
+ const data<MetaSchema, Encode> meta() const {
+ return inner_.get_dims();
}
template<uint64_t i>
diff --git a/c++/equilibrium.hpp b/c++/equilibrium.hpp
index 326fc9e..ac36dbc 100644
--- a/c++/equilibrium.hpp
+++ b/c++/equilibrium.hpp
@@ -5,7 +5,7 @@
namespace kel {
namespace lbm {
template<typename T, typename Descriptor>
-saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<sch::T> rho, saw::data<sch::FixedArray<T,Descriptor::D>> vel){
+saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<T> rho, saw::data<sch::FixedArray<T,Descriptor::D>> vel){
using dfi = df_info<T, Descriptor>;
saw::data<sch::FixedArray<T,Descriptor::Q>> eq;
@@ -18,26 +18,26 @@ saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<sch::T> rho,
vel_vel = vel_vel + vel.at(j) * vel.at(j);
}
- for(uint64_t i = 0u; i < eq.get_dim_size<0u>(); ++i){
+ for(uint64_t i = 0u; i < eq.template get_dim_size<0u>(); ++i){
saw::data<T> vel_c{};
for(uint64_t j = 0u; j < Descriptor::D; ++j){
- vel_c = vel_c + (vel.at(j) * {dfi::directions[i][j]});
+ vel_c = vel_c + (vel.at(j) * saw::data<T>{static_cast<saw::native_data_type<T>::type>(dfi::directions[i][j])});
}
- auto vel_c_cs2 = vel_c * {dfi::inv_cs2};
+ auto vel_c_cs2 = vel_c * saw::data<T>{dfi::inv_cs2};
eq.at(i).set(
dfi::weights[i] * rho.get() *
(
1.0
- + vel_c_cs2
+ + vel_c_cs2.get()
- dfi::inv_cs2 * 0.5 * vel_vel.get()
- + vel_c_cs2 * vel_c_cs2 * 0.5
+ + vel_c_cs2.get() * vel_c_cs2.get() * 0.5
)
);
-
- return eq;
}
+
+ return eq;
}
}
}
diff --git a/c++/examples/cavity_2d.cpp b/c++/examples/cavity_2d.cpp
index 49f47cd..88a7702 100644
--- a/c++/examples/cavity_2d.cpp
+++ b/c++/examples/cavity_2d.cpp
@@ -38,11 +38,12 @@ using CellInfo = Cell<UInt8, D2Q9, 1, 0, 0>;
template<typename Desc>
using CellStruct = Struct<
Member<DfCell<Desc>, "dfs">,
+ Member<DfCell<Desc>, "dfs_old">,
Member<CellInfo<Desc>, "info">
>;
-using CavityFieldD2Q9 = CellFields<D2Q9, CellStruct<D2Q9>>;
+using CavityFieldD2Q9 = CellField<D2Q9, CellStruct<D2Q9>>;
}
/**
@@ -237,8 +238,8 @@ public:
}
constexpr size_t dim_size = 2;
-constexpr size_t dim_x = 256;
-constexpr size_t dim_y = 256;
+constexpr size_t dim_x = 128;
+constexpr size_t dim_y = 128;
template<typename Func>
@@ -247,7 +248,7 @@ void apply_for_cells(Func&& func, saw::data<kel::lbm::sch::CavityFieldD2Q9>& dat
for(std::size_t j = 0; j < dat.meta().at({0u}).get(); ++j){
saw::data<saw::schema::UInt64> di{i};
saw::data<saw::schema::UInt64> dj{j};
- auto& cell_v = dat.template get<"dfs">().at({{dj,di}});
+ auto& cell_v = dat({{dj,di}});
func(cell_v, j, i);
}
}
@@ -285,9 +286,11 @@ void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
(void) i;
(void) j;
auto& dfs = cell.template get<"dfs">();
+ auto& dfs_old = cell.template get<"dfs_old">();
auto info = cell.template get<"info">()(0u).get();
for(uint64_t k = 0; k < sch::D2Q9::Q; ++k){
dfs(k).set(eq[k]);
+ dfs_old(k).set(eq[k]);
}
}, latt);
}
@@ -300,10 +303,12 @@ void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
(void) i;
(void) j;
auto& dfs = cell.template get<"dfs">();
+ auto& dfs_old = cell.template get<"dfs_old">();
auto info = cell.template get<"info">()(0u).get();
if(info == 2u){
for(uint64_t k = 0; k < sch::D2Q9::Q; ++k){
dfs(k).set(eq[k]);
+ dfs_old(k).set(eq[k]);
}
}
}, latt);
@@ -311,16 +316,16 @@ void set_initial_conditions(saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt){
}
void lbm_step(
- saw::data<kel::lbm::sch::CavityFieldD2Q9>& old_latt,
- saw::data<kel::lbm::sch::CavityFieldD2Q9>& new_latt
+ saw::data<kel::lbm::sch::CavityFieldD2Q9>& latt,
+ bool even_step
){
using namespace kel::lbm;
using dfi = df_info<sch::T,sch::D2Q9>;
component<cmpt::BGK,sch::D2Q9> coll;
- coll.relaxation_ = 0.653846;
+ coll.relaxation_ = 0.5384;
component<cmpt::ConstRhoBGK,sch::D2Q9> rho_coll;
- rho_coll.relaxation_ = 0.501538;
+ rho_coll.relaxation_ = 0.5384;
rho_coll.rho_ = 1.0;
component<cmpt::BounceBack,sch::D2Q9> bb;
@@ -329,13 +334,13 @@ void lbm_step(
// Collide
apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){
- auto& df = cell.template get<"dfs">();
+ auto& df = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
auto& info = cell.template get<"info">();
auto info_val = info({0u}).get();
switch(info_val){
case 1u:
- rho_coll.apply(df);
+ coll.apply(df);
break;
case 2u:
// bb.apply(df);
@@ -345,23 +350,25 @@ void lbm_step(
bb.apply(df);
break;
}
- }, old_latt);
+ }, latt);
// Stream
- for(uint64_t i = 1; (i+1) < old_latt.template get_dim_size<0>().get(); ++i){
- for(uint64_t j = 1; (j+1) < old_latt.template get_dim_size<1>().get(); ++j){
- auto& df_new = new_latt.template get<"dfs">().at({{i,j}});
- auto& info_new = new_latt.template get<"info">().at({{i,j}});
+ for(uint64_t i = 1; (i+1) < latt.template get_dim_size<0>().get(); ++i){
+ for(uint64_t j = 1; (j+1) < latt.template get_dim_size<1>().get(); ++j){
+ auto& cell = latt({{i,j}});
+ auto& df_new = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">();
+ auto& info_new = cell.template get<"info">();
if(info_new({0u}).get() > 0u && info_new({0u}).get() != 2u){
for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){
auto dir = dfi::directions[dfi::opposite_index[k]];
+ auto& cell_dir_old = latt({{i+dir[0],j+dir[1]}});
- auto& df_old = old_latt.template get<"dfs">().at({{i+dir[0],j+dir[1]}});
- auto& info_old = old_latt.template get<"info">().at({{i+dir[0],j+dir[1]}});
+ auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">();
+ auto& info_old = cell_dir_old.template get<"info">();
if( info_old({0}).get() == 2u ){
- auto& df_old_loc = old_latt.template get<"dfs">().at({{i,j}});
+ auto& df_old_loc = even_step ? latt({{i,j}}).template get<"dfs_old">() : latt({{i,j}}).template get<"dfs">();
df_new({k}) = df_old_loc({dfi::opposite_index.at(k)}) - 2.0 * dfi::inv_cs2 * dfi::weights.at(k) * 1.0 * ( bb_lid.lid_vel[0] * dir[0] + bb_lid.lid_vel[1] * dir[1]);
// dfs({dfi::opposite_index.at(i)}) = dfs_cpy({i}) - 2.0 * dfi::weights[i] * 1.0 * ( lid_vel[0] * dfi::directions[i][0] + lid_vel[1] * dfi::directions[i][1]) * dfi::inv_cs2;
} else {
@@ -378,8 +385,7 @@ int main(){
saw::data<sch::FixedArray<sch::UInt64,sch::D2Q9::D>> dim{{dim_x, dim_y}};
- saw::data<sch::CavityFieldD2Q9, saw::encode::Native> old_lattice{dim};
- saw::data<sch::CavityFieldD2Q9, saw::encode::Native> new_lattice{dim};
+ saw::data<sch::CavityFieldD2Q9, saw::encode::Native> lattice{dim};
// auto& df_field = lattices.at(0).template get<"dfs">();
//for(uint64_t i = 0; i < df_field.get_dim_size<0u>(); ++i){
@@ -389,13 +395,11 @@ int main(){
/**
* Set meta information describing what this cell is
*/
- set_geometry(old_lattice);
- set_geometry(new_lattice);
+ set_geometry(lattice);
/**
*
*/
- set_initial_conditions(old_lattice);
- set_initial_conditions(new_lattice);
+ set_initial_conditions(lattice);
/**
* Timeloop
@@ -413,19 +417,15 @@ int main(){
}else{
std::cout<<"\n";
}
- }, old_lattice);
+ }, lattice);
- uint64_t lattice_steps = 74000u;
+ uint64_t lattice_steps = 512000u;
bool even_step = true;
uint64_t print_every = 256u;
uint64_t file_no = 0u;
for(uint64_t step = 0; step < lattice_steps; ++step){
- auto& old_lat = even_step ? old_lattice : new_lattice;
- auto& new_lat = even_step ? new_lattice : old_lattice;
-
-
if(step % print_every == 0u){
std::cout<<"\n";
@@ -479,7 +479,7 @@ int main(){
}else{
std::cout<<"\n";
}
- }, old_lat);
+ }, lattice);
std::cout<<"Average rho: "<<(sum / ((dim_x-4)*(dim_y-4)))<<std::endl;
@@ -496,20 +496,20 @@ int main(){
vtk_file <<"VECTORS Velocity float\n";
apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){
- auto dfs = cell.template get<"dfs">();
+ auto dfs = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">();
typename saw::native_data_type<sch::T>::type rho;
std::array<typename saw::native_data_type<sch::T>::type, sch::D2Q9::D> vel;
compute_rho_u<sch::D2Q9>(dfs,rho,vel);
vtk_file << static_cast<float>(vel[0u]) << " " << static_cast<float>(vel[1u])<<" 0.0\n";
- }, old_lat);
+ }, lattice);
++file_no;
}
- lbm_step(old_lat, new_lat);
+ lbm_step(lattice, even_step);
even_step = !even_step;
}
diff --git a/tests/equilibrium.cpp b/tests/equilibrium.cpp
index bb26dd0..9201e55 100644
--- a/tests/equilibrium.cpp
+++ b/tests/equilibrium.cpp
@@ -9,17 +9,17 @@ template<typename Descriptor>
void check_equilibrium(){
using namespace kel;
- using dfi = df_info<lbm::sch::Float64,Descriptor>;
+ using dfi = lbm::df_info<lbm::sch::Float64,Descriptor>;
saw::data<lbm::sch::Float64> rho{1.0};
- saw::data<lbm::sch::FixedArray<lbm::sch::Float64,>> vel;
- for(saw::data<lbm::sch::UInt64> i{0u}; i < {Descriptor::D}; ++i){
+ saw::data<lbm::sch::FixedArray<lbm::sch::Float64,Descriptor::D>> vel;
+ for(saw::data<lbm::sch::UInt64> i{0u}; i.get() < Descriptor::D; ++i){
vel.at(i) = {0.0};
}
auto eq = lbm::equilibrium<lbm::sch::Float64,Descriptor>(rho,vel);
- for(saw::data<lbm::sch::UInt64> i{0u}; i < {Descriptor::Q}; ++i){
- SAW_CHECK(eq(i) == dfi::weights[i.get()], "No velocity and normalized rho should be exactly the weights");
+ 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()]));
}
}
@@ -27,4 +27,14 @@ SAW_TEST("Equilibrium at rest D1Q3"){
using namespace kel;
check_equilibrium<lbm::sch::Descriptor<1,3>>();
}
+
+SAW_TEST("Equilibrium at rest D2Q5"){
+ using namespace kel;
+ check_equilibrium<lbm::sch::Descriptor<2,5>>();
+}
+
+SAW_TEST("Equilibrium at rest D2Q9"){
+ using namespace kel;
+ check_equilibrium<lbm::sch::Descriptor<2,9>>();
+}
}
generated by cgit v1.2.3 (git 2.25.1) at 2025-04-19 23:04:43 +0000