Merge branch 'dev'

1 files changed, 523 insertions, 22 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;
 		}
 	}
 };