summaryrefslogtreecommitdiff
path: root/modules
diff options
context:
space:
mode:
Diffstat (limited to 'modules')
-rw-r--r--modules/core/c++/boundary.hpp545
-rw-r--r--modules/core/c++/boundary/common.hpp8
-rw-r--r--modules/core/c++/boundary/zou_he_pressure.hpp9
-rw-r--r--modules/core/c++/collision.hpp4
-rw-r--r--modules/core/c++/component/psm.hpp2
-rw-r--r--modules/core/c++/fplbm.hpp10
-rw-r--r--modules/core/c++/hlbm.hpp4
-rw-r--r--modules/core/c++/macroscopic.hpp1
-rw-r--r--modules/core/c++/particle.hpp2
-rw-r--r--modules/core/c++/psm.hpp516
-rw-r--r--modules/core/c++/write_csv.hpp2
-rw-r--r--modules/core/c++/write_vtk.hpp2
-rw-r--r--modules/core/tests/equilibrium.cpp2
-rw-r--r--modules/core/tests/vtk_write.cpp10
14 files changed, 969 insertions, 148 deletions
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;
generated by cgit v1.2.3 (git 2.25.1) at 2026-08-30 11:22:35 +0000