diff --git a/doc/notebooks/05_tutorial_nondimensionalization_and_scaling.ipynb b/doc/notebooks/05_tutorial_nondimensionalization_and_scaling.ipynb
index d40de1bd7cb3f16764b1cd23a9b880726d0f23cf..07066795ad8dc8311e9ce9d3dba438b25d3734bc 100644
--- a/doc/notebooks/05_tutorial_nondimensionalization_and_scaling.ipynb
+++ b/doc/notebooks/05_tutorial_nondimensionalization_and_scaling.ipynb
@@ -52,7 +52,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
- "model_id": "8013073325cc4f7da83800d88fc944bf",
+ "model_id": "267d2750d0114490bce98ee465856b8b",
"version_major": 2,
"version_minor": 0
},
@@ -385,7 +385,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Python 3",
+ "display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
@@ -399,7 +399,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.8.2"
+ "version": "3.9.9"
}
},
"nbformat": 4,
diff --git a/lbmpy/methods/momentbased/recursiveregularisedmethod.py b/lbmpy/methods/momentbased/recursiveregularisedmethod.py
index b48715b0ca8c4c0aea6f9b914d7fe5186c6cd1c1..08d56557ed663f90a69c8525359e2208e3aa9d50 100644
--- a/lbmpy/methods/momentbased/recursiveregularisedmethod.py
+++ b/lbmpy/methods/momentbased/recursiveregularisedmethod.py
@@ -335,6 +335,8 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
if additional_subexpressions is None:
additional_subexpressions = list()
+ half = sp.Rational(1, 2)
+
omega1 = d[4, 4]
omega2 = d[7, 7]
omega3 = d[10, 10]
@@ -356,10 +358,25 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
cs4 = c_s_sq * c_s_sq
u = sp.symbols("u_:3")
- taup = rho * c_s_sq / omega1
+ a0xx, a0xy, a0yy, a0xz, a0zz, a0yz, a0xxy, a0xyy, a0xxz, a0xzz, a0yyz, a0yzz = sp.symbols("a0xx a0xy a0yy a0xz a0zz a0yz a0xxy a0xyy a0xxz a0xzz a0yyz a0yzz")
+ a1xx, a1xy, a1yy, a1xz, a1zz, a1yz, a1xxy, a1xyy, a1xxz, a1xzz, a1yyz, a1yzz = sp.symbols("a1xx a1xy a1yy a1xz a1zz a1yz a1xxy a1xyy a1xxz a1xzz a1yyz a1yzz")
+
+ tauS = 1 / omega1
+
+ a12FD = [sp.Symbol(f"a12FD_{i}") for i in range(6)]
+ p1ij = [sp.Symbol(f"p1ij_{i}") for i in range(6)]
+ P = [sp.Symbol(f"P{i}") for i in range(6)]
+ inv36 = sp.Rational(1, 36)
+ inv72 = sp.Rational(1, 72)
+ inv12 = sp.Rational(1, 12)
+ inv24 = sp.Rational(1, 24)
+
+
+
RR = [sp.Symbol(f"M_{index.name[1:]}") for index in Ordering]
- RReq = [sp.Symbol(f"eq_{index.name[1:]}") for index in Ordering]
+ fEq = [sp.Symbol(f"fEq_{index.name[1:]}") for index in Ordering]
+ fNEq = [sp.Symbol(f"fNEq_{index.name[1:]}") for index in Ordering]
# RMeq = [sp.Symbol(f"RMeq_{index.name[1:]}") for index in Ordering]
# eq = sp.symbols(f"eq_:{self.stencil.Q}")
RRneq = [sp.Symbol(f"neq_{index.name[1:]}") for index in Ordering]
@@ -368,25 +385,25 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
stencil = self.stencil
F000 = stencil.index((0, 0, 0))
- FM00 = stencil.index((-1, 0, 0))
- F0M0 = stencil.index((0, -1, 0))
- F00M = stencil.index((0, 0, -1))
- FMM0 = stencil.index((-1, -1, 0))
- FMP0 = stencil.index((-1, 1, 0))
- FM0M = stencil.index((-1, 0, -1))
- FM0P = stencil.index((-1, 0, 1))
- F0MM = stencil.index((0, -1, -1))
+ FM00 = stencil.index((0, 0, 1))
+ F0M0 = stencil.index((0, 0, -1))
+ F00M = stencil.index((0, 1, 0))
+ FMM0 = stencil.index((0, -1, 0))
+ FMP0 = stencil.index((1, 0, 0))
+ FM0M = stencil.index((-1, 0, 0))
+ FM0P = stencil.index((0, 1, 1))
+ F0MM = stencil.index((0, 1, -1))
F0MP = stencil.index((0, -1, 1))
- FP00 = stencil.index((1, 0, 0))
- F0P0 = stencil.index((0, 1, 0))
- F00P = stencil.index((0, 0, 1))
- FPP0 = stencil.index((1, 1, 0))
- FPM0 = stencil.index((1, -1, 0))
- FP0P = stencil.index((1, 0, 1))
- FP0M = stencil.index((1, 0, -1))
- F0PP = stencil.index((0, 1, 1))
- F0PM = stencil.index((0, 1, -1))
+ FP00 = stencil.index((0, -1, -1))
+ F0P0 = stencil.index((1, 0, 1))
+ F00P = stencil.index((1, 0, -1))
+ FPP0 = stencil.index((-1, 0, 1))
+ FPM0 = stencil.index((-1, 0, -1))
+ FP0P = stencil.index((1, 1, 0))
+ FP0M = stencil.index((1, -1, 0))
+ F0PP = stencil.index((-1, 1, 0))
+ F0PM = stencil.index((-1, -1, 0))
M000 = Ordering.M000
@@ -415,6 +432,10 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
gradient = []
sr = []
+ ux = u[0]
+ uy = u[1]
+ uz = u[2]
+
if self._velocity_field:
gradient = [Assignment(sp.Symbol("du_11"), self._velocity_field[1, 0, 0](0) - self._velocity_field[-1, 0, 0](0)),
Assignment(sp.Symbol("du_12"), self._velocity_field[0, 1, 0](0) - self._velocity_field[0, -1, 0](0)),
@@ -426,12 +447,12 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
Assignment(sp.Symbol("du_32"), self._velocity_field[0, 1, 0](2) - self._velocity_field[0, -1, 0](2)),
Assignment(sp.Symbol("du_33"), self._velocity_field[0, 0, 1](2) - self._velocity_field[0, 0, -1](2))]
- sr = [Assignment(sp.Symbol("SR1"), 2.0 * sp.Symbol("du_11")),
- Assignment(sp.Symbol("SR2"), 2.0 * sp.Symbol("du_22")),
- Assignment(sp.Symbol("SR3"), 2.0 * sp.Symbol("du_33")),
- Assignment(sp.Symbol("SR4"), sp.Symbol("du_12") + sp.Symbol("du_21")),
- Assignment(sp.Symbol("SR5"), sp.Symbol("du_13") + sp.Symbol("du_31")),
- Assignment(sp.Symbol("SR6"), sp.Symbol("du_23") + sp.Symbol("du_32"))]
+ sr = [Assignment(sp.Symbol("SR1"), sp.Symbol("du_11")),
+ Assignment(sp.Symbol("SR2"), sp.Symbol("du_22")),
+ Assignment(sp.Symbol("SR3"), sp.Symbol("du_33")),
+ Assignment(sp.Symbol("SR4"), half * (sp.Symbol("du_12") + sp.Symbol("du_21"))),
+ Assignment(sp.Symbol("SR5"), half * (sp.Symbol("du_13") + sp.Symbol("du_31"))),
+ Assignment(sp.Symbol("SR6"), half * (sp.Symbol("du_23") + sp.Symbol("du_32")))]
moments = [Assignment(X_M1, f[FM00] + f[FMM0] + f[FMP0] + f[FM0M] + f[FM0P]),
Assignment(X_P1, f[FP00] + f[FPP0] + f[FPM0] + f[FP0P] + f[FP0M]),
@@ -445,162 +466,190 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
Assignment(inv_rho, 1.0 / rho),
Assignment(RR[M000], 1.0),
- Assignment(RR[M100], inv_rho * (X_P1 - X_M1)),
- Assignment(RR[M010], inv_rho * (Y_P1 - Y_M1)),
- Assignment(RR[M001], inv_rho * (Z_P1 - Z_M1)),
- Assignment(RR[M200], inv_rho * (X_M1 + X_P1) - c_s_sq),
- Assignment(RR[M020], inv_rho * (Y_M1 + Y_P1) - c_s_sq),
- Assignment(RR[M002], inv_rho * (Z_M1 + Z_P1) - c_s_sq),
- Assignment(RR[M110], inv_rho * (f[FMM0] - f[FMP0] + f[FPP0] - f[FPM0])),
- Assignment(RR[M101], inv_rho * (f[FM0M] - f[FM0P] + f[FP0P] - f[FP0M])),
- Assignment(RR[M011], inv_rho * (f[F0MM] - f[F0MP] + f[F0PP] - f[F0PM])),
+ Assignment(p1ij[0], inv_rho * (X_M1 + X_P1) - c_s_sq - u[0] * u[0]),
+ Assignment(p1ij[1], inv_rho * (Y_M1 + Y_P1) - c_s_sq - u[1] * u[1]),
+ Assignment(p1ij[2], inv_rho * (Z_M1 + Z_P1) - c_s_sq - u[2] * u[2]),
+ Assignment(p1ij[3], inv_rho * (f[FMM0] - f[FMP0] + f[FPP0] - f[FPM0]) - u[0] * u[1]),
+ Assignment(p1ij[4], inv_rho * (f[FM0M] - f[FM0P] + f[FP0P] - f[FP0M]) - u[0] * u[2]),
+ Assignment(p1ij[5], inv_rho * (f[F0MM] - f[F0MP] + f[F0PP] - f[F0PM]) - u[1] * u[2]),
+
+ Assignment(a12FD[0], -2.0 * tauS * rho * c_s_sq * sp.Symbol("SR1")),
+ Assignment(a12FD[1], -2.0 * tauS * rho * c_s_sq * sp.Symbol("SR2")),
+ Assignment(a12FD[2], -2.0 * tauS * rho * c_s_sq * sp.Symbol("SR3")),
+ Assignment(a12FD[3], -2.0 * tauS * rho * c_s_sq * sp.Symbol("SR4")),
+ Assignment(a12FD[4], -2.0 * tauS * rho * c_s_sq * sp.Symbol("SR5")),
+ Assignment(a12FD[5], -2.0 * tauS * rho * c_s_sq * sp.Symbol("SR6")),
]
- eq_moments = [Assignment(RReq[M000], 1.0),
- Assignment(RReq[M100], u[0]),
- Assignment(RReq[M010], u[1]),
- Assignment(RReq[M001], u[2]),
- Assignment(RReq[M200], u[0] * u[0]),
- Assignment(RReq[M020], u[1] * u[1]),
- Assignment(RReq[M002], u[2] * u[2]),
- Assignment(RReq[M110], u[0] * u[1]),
- Assignment(RReq[M101], u[0] * u[2]),
- Assignment(RReq[M011], u[1] * u[2]),
- Assignment(RReq[M210], RReq[M200] * u[1]),
- Assignment(RReq[M201], RReq[M200] * u[2]),
- Assignment(RReq[M021], RReq[M020] * u[2]),
- Assignment(RReq[M120], RReq[M020] * u[0]),
- Assignment(RReq[M102], RReq[M002] * u[0]),
- Assignment(RReq[M012], RReq[M002] * u[1]),
- Assignment(RReq[M220], RReq[M200] * RReq[M020]),
- Assignment(RReq[M202], RReq[M200] * RReq[M002]),
- Assignment(RReq[M022], RReq[M020] * RReq[M002])]
-
neq = [
- Assignment(RRneq[M200], (RR[M200] - RReq[M200]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR1"))),
- Assignment(RRneq[M020], (RR[M020] - RReq[M020]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR2"))),
- Assignment(RRneq[M002], (RR[M002] - RReq[M002]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR3"))),
- Assignment(RRneq[M110], (RR[M110] - RReq[M110]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR4"))),
- Assignment(RRneq[M101], (RR[M101] - RReq[M101]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR5"))),
- Assignment(RRneq[M011], (RR[M011] - RReq[M011]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR6"))),
-
- Assignment(RRneq[M210], u[1] * RRneq[M200] + 2.0 * u[0] * RRneq[M110]),
- Assignment(RRneq[M201], u[2] * RRneq[M200] + 2.0 * u[0] * RRneq[M101]),
- Assignment(RRneq[M021], u[2] * RRneq[M020] + 2.0 * u[1] * RRneq[M011]),
- Assignment(RRneq[M120], u[0] * RRneq[M020] + 2.0 * u[1] * RRneq[M110]),
- Assignment(RRneq[M102], u[0] * RRneq[M002] + 2.0 * u[2] * RRneq[M101]),
- Assignment(RRneq[M012], u[1] * RRneq[M002] + 2.0 * u[2] * RRneq[M011]),
-
- Assignment(RRneq[M220],
- u[1] * u[1] * RRneq[M200] + u[0] * u[0] * RRneq[M020] + 4.0 * u[0] * u[1] *
- RRneq[M110]),
- Assignment(RRneq[M202],
- u[2] * u[2] * RRneq[M200] + u[0] * u[0] * RRneq[M002] + 4.0 * u[0] * u[2] *
- RRneq[M101]),
- Assignment(RRneq[M022],
- u[2] * u[2] * RRneq[M020] + u[1] * u[1] * RRneq[M002] + 4.0 * u[1] * u[2] *
- RRneq[M011])
- ]
+ Assignment(P[0], p1ij[0] * sigma + (sp.Integer(1) - sigma) * a12FD[0]),
+ Assignment(P[1], p1ij[1] * sigma + (sp.Integer(1) - sigma) * a12FD[1]),
+ Assignment(P[2], p1ij[2] * sigma + (sp.Integer(1) - sigma) * a12FD[2]),
+ Assignment(P[3], p1ij[3] * sigma + (sp.Integer(1) - sigma) * a12FD[3]),
+ Assignment(P[4], p1ij[4] * sigma + (sp.Integer(1) - sigma) * a12FD[4]),
+ Assignment(P[5], p1ij[5] * sigma + (sp.Integer(1) - sigma) * a12FD[5]),
+
+ Assignment(a0xx, ux * ux),
+ Assignment(a0xy, ux * uy),
+ Assignment(a0yy, uy * uy),
+ Assignment(a0xz, ux * uz),
+ Assignment(a0zz, uz * uz),
+ Assignment(a0yz, uy * uz),
+
+ Assignment(a0xxy, uy * a0xx),
+ Assignment(a0xyy, ux * a0yy),
+ Assignment(a0xxz, uz * a0xx),
+ Assignment(a0xzz, ux * a0zz),
+ Assignment(a0yyz, uz * a0yy),
+ Assignment(a0yzz, uy * a0zz),
+
+ Assignment(fEq[0],
+ -((sp.Float(3.0) * a0zz + sp.Float(3.0) * a0yy + sp.Float(3.0) * a0xx - 2) * rho) * sp.Float(
+ 1 / 6)),
+ Assignment(fEq[1], (rho * (
+ sp.Float(6.0) * uz + sp.Float(6.0) * a0zz - sp.Float(9.0) * a0yyz - sp.Float(3.0) * a0yy - sp.Float(
+ 9.0) * a0xxz - sp.Float(3.0) * a0xx + 2)) * (inv36)),
+ Assignment(fEq[2], -(rho * (
+ sp.Float(6.0) * uz - sp.Float(6.0) * a0zz - sp.Float(9.0) * a0yyz + sp.Float(3.0) * a0yy - sp.Float(
+ 9.0) * a0xxz + sp.Float(3.0) * a0xx - 2)) * (inv36)),
+ Assignment(fEq[3], (rho * (
+ sp.Float(6.0) * uy - sp.Float(3.0) * a0zz - sp.Float(9.0) * a0yzz + sp.Float(6.0) * a0yy - sp.Float(
+ 9.0) * a0xxy - sp.Float(3.0) * a0xx + 2)) * (inv36)),
+ Assignment(fEq[4], -(rho * (
+ sp.Float(6.0) * uy + sp.Float(3.0) * a0zz - sp.Float(9.0) * a0yzz - sp.Float(6.0) * a0yy - sp.Float(
+ 9.0) * a0xxy + sp.Float(3.0) * a0xx - 2)) * (inv36)),
+ Assignment(fEq[5], (rho * (
+ sp.Float(6.0) * ux - sp.Float(3.0) * a0zz - sp.Float(3.0) * a0yy - sp.Float(9.0) * a0xzz - sp.Float(
+ 9.0) * a0xyy + sp.Float(6.0) * a0xx + 2)) * (inv36)),
+ Assignment(fEq[6], -(rho * (
+ sp.Float(6.0) * ux + sp.Float(3.0) * a0zz + sp.Float(3.0) * a0yy - sp.Float(9.0) * a0xzz - sp.Float(
+ 9.0) * a0xyy - sp.Float(6.0) * a0xx - 2)) * (inv36)),
+ Assignment(fEq[7], (rho * (sp.Float(6.0) * uz + sp.Float(6.0) * uy + sp.Float(
+ 6.0) * a0zz + 18 * a0yzz + 18 * a0yz + 18 * a0yyz + sp.Float(6.0) * a0yy - sp.Float(
+ 9.0) * a0xxz - sp.Float(
+ 9.0) * a0xxy - sp.Float(3.0) * a0xx + 2)) * (inv72)),
+ Assignment(fEq[8], -(rho * (sp.Float(6.0) * uz - sp.Float(6.0) * uy - sp.Float(
+ 6.0) * a0zz - 18 * a0yzz + 18 * a0yz + 18 * a0yyz - sp.Float(6.0) * a0yy - sp.Float(
+ 9.0) * a0xxz + sp.Float(
+ 9.0) * a0xxy + sp.Float(3.0) * a0xx - 2)) * (inv72)),
+ Assignment(fEq[9], (rho * (sp.Float(6.0) * uz - sp.Float(6.0) * uy + sp.Float(
+ 6.0) * a0zz - 18 * a0yzz - 18 * a0yz + 18 * a0yyz + sp.Float(6.0) * a0yy - sp.Float(
+ 9.0) * a0xxz + sp.Float(
+ 9.0) * a0xxy - sp.Float(3.0) * a0xx + 2)) * (inv72)),
+ Assignment(fEq[10], -(rho * (sp.Float(6.0) * uz + sp.Float(6.0) * uy - sp.Float(
+ 6.0) * a0zz + 18 * a0yzz - 18 * a0yz + 18 * a0yyz - sp.Float(6.0) * a0yy - sp.Float(
+ 9.0) * a0xxz - sp.Float(
+ 9.0) * a0xxy + sp.Float(3.0) * a0xx - 2)) * (inv72)),
+ Assignment(fEq[11], (rho * (
+ sp.Float(6.0) * uz + sp.Float(6.0) * ux + sp.Float(6.0) * a0zz - sp.Float(9.0) * a0yyz - sp.Float(
+ 3.0) * a0yy + 18 * a0xzz + 18 * a0xz - sp.Float(9.0) * a0xyy + 18 * a0xxz + sp.Float(
+ 6.0) * a0xx + 2)) * (
+ inv72)),
+ Assignment(fEq[12], -(rho * (
+ sp.Float(6.0) * uz - sp.Float(6.0) * ux - sp.Float(6.0) * a0zz - sp.Float(9.0) * a0yyz + sp.Float(
+ 3.0) * a0yy - 18 * a0xzz + 18 * a0xz + sp.Float(9.0) * a0xyy + 18 * a0xxz - sp.Float(
+ 6.0) * a0xx - 2)) * (
+ inv72)),
+ Assignment(fEq[13], (rho * (sp.Float(6.0) * uz - sp.Float(6.0) * ux + sp.Float(6.0) * a0zz - sp.Float(9.0) * a0yyz - sp.Float(3.0) * a0yy - 18 * a0xzz - 18 * a0xz + sp.Float(9.0) * a0xyy + 18 * a0xxz + sp.Float(6.0) * a0xx + 2)) * (inv72)),
+ Assignment(fEq[14], -(rho * (
+ sp.Float(6.0) * uz + sp.Float(6.0) * ux - sp.Float(6.0) * a0zz - sp.Float(9.0) * a0yyz + sp.Float(
+ 3.0) * a0yy + 18 * a0xzz - 18 * a0xz - sp.Float(9.0) * a0xyy + 18 * a0xxz - sp.Float(
+ 6.0) * a0xx - 2)) * (
+ inv72)),
+ Assignment(fEq[15], (rho * (
+ sp.Float(6.0) * uy + sp.Float(6.0) * ux - sp.Float(3.0) * a0zz - sp.Float(9.0) * a0yzz + sp.Float(
+ 6.0) * a0yy - sp.Float(9.0) * a0xzz + 18 * a0xyy + 18 * a0xy + 18 * a0xxy + sp.Float(
+ 6.0) * a0xx + 2)) * (
+ inv72)),
+ Assignment(fEq[16], -(rho * (
+ sp.Float(6.0) * uy - sp.Float(6.0) * ux + sp.Float(3.0) * a0zz - sp.Float(9.0) * a0yzz - sp.Float(
+ 6.0) * a0yy + sp.Float(9.0) * a0xzz - 18 * a0xyy + 18 * a0xy + 18 * a0xxy - sp.Float(
+ 6.0) * a0xx - 2)) * (
+ inv72)),
+ Assignment(fEq[17], (rho * (
+ sp.Float(6.0) * uy - sp.Float(6.0) * ux - sp.Float(3.0) * a0zz - sp.Float(9.0) * a0yzz + sp.Float(
+ 6.0) * a0yy + sp.Float(9.0) * a0xzz - 18 * a0xyy - 18 * a0xy + 18 * a0xxy + sp.Float(
+ 6.0) * a0xx + 2)) * (
+ inv72)),
+ Assignment(fEq[18], -(rho * (
+ sp.Float(6.0) * uy + sp.Float(6.0) * ux + sp.Float(3.0) * a0zz - sp.Float(9.0) * a0yzz - sp.Float(
+ 6.0) * a0yy - sp.Float(9.0) * a0xzz + 18 * a0xyy - 18 * a0xy + 18 * a0xxy - sp.Float(
+ 6.0) * a0xx - 2)) * (
+ inv72))
- col = [
- Assignment(RRcoll[M200], RR[M200] - omegaPlus * RRneq[M200] - omegaMinus * RRneq[
- M020] - omegaMinus * RRneq[M002]),
- Assignment(RRcoll[M020], RR[M020] - omegaMinus * RRneq[M200] - omegaPlus * RRneq[
- M020] - omegaMinus * RRneq[M002]),
- Assignment(RRcoll[M002],
- RR[M002] - omegaMinus * RRneq[M200] - omegaMinus * RRneq[
- M020] - omegaPlus * RRneq[M002]),
-
- Assignment(RRcoll[M110], (1.0 - omega2) * RRneq[M110] + RReq[M110]),
- Assignment(RRcoll[M101], (1.0 - omega2) * RRneq[M101] + RReq[M101]),
- Assignment(RRcoll[M011], (1.0 - omega2) * RRneq[M011] + RReq[M011]),
-
- Assignment(RRcoll[M210], (1.0 - omega3) * RRneq[M210] + RReq[M210]),
- Assignment(RRcoll[M201], (1.0 - omega3) * RRneq[M201] + RReq[M201]),
- Assignment(RRcoll[M021], (1.0 - omega3) * RRneq[M021] + RReq[M021]),
- Assignment(RRcoll[M120], (1.0 - omega3) * RRneq[M120] + RReq[M120]),
- Assignment(RRcoll[M102], (1.0 - omega3) * RRneq[M102] + RReq[M102]),
- Assignment(RRcoll[M012], (1.0 - omega3) * RRneq[M012] + RReq[M012]),
-
- Assignment(RRcoll[M220], (1.0 - omega4) * RRneq[M220] + RReq[M220]),
- Assignment(RRcoll[M202], (1.0 - omega4) * RRneq[M202] + RReq[M202]),
- Assignment(RRcoll[M022], (1.0 - omega4) * RRneq[M022] + RReq[M022]),
-
- Assignment(RMcoll[M200], RRcoll[M200] + c_s_sq),
- Assignment(RMcoll[M020], RRcoll[M020] + c_s_sq),
- Assignment(RMcoll[M002], RRcoll[M002] + c_s_sq),
-
- Assignment(RMcoll[M110], RRcoll[M110]),
- Assignment(RMcoll[M101], RRcoll[M101]),
- Assignment(RMcoll[M011], RRcoll[M011]),
-
- Assignment(RMcoll[M210], RRcoll[M210] + c_s_sq * u[1]),
- Assignment(RMcoll[M201], RRcoll[M201] + c_s_sq * u[2]),
- Assignment(RMcoll[M021], RRcoll[M021] + c_s_sq * u[2]),
- Assignment(RMcoll[M120], RRcoll[M120] + c_s_sq * u[0]),
- Assignment(RMcoll[M102], RRcoll[M102] + c_s_sq * u[0]),
- Assignment(RMcoll[M012], RRcoll[M012] + c_s_sq * u[1]),
-
- Assignment(RMcoll[M220],
- RRcoll[M220] + c_s_sq * (RRcoll[M200] + RRcoll[M020]) + cs4),
- Assignment(RMcoll[M202],
- RRcoll[M202] + c_s_sq * (RRcoll[M200] + RRcoll[M002]) + cs4),
- Assignment(RMcoll[M022],
- RRcoll[M022] + c_s_sq * (RRcoll[M020] + RRcoll[M002]) + cs4),
]
- tmp1 = sp.Rational(1, 2) * rho * (
- u[0] + RMcoll[M200] - RMcoll[M120] - RMcoll[M102] - RMcoll[
- M220] - RMcoll[M202])
- tmp2 = sp.Rational(1, 2) * rho * (
- u[1] + RMcoll[M020] - RMcoll[M210] - RMcoll[M012] - RMcoll[
- M220] - RMcoll[M022])
- tmp3 = sp.Rational(1, 2) * rho * (
- u[2] + RMcoll[M002] - RMcoll[M201] - RMcoll[M021] - RMcoll[
- M202] - RMcoll[M022])
- tmp4 = sp.Rational(1, 4) * rho * (
- RMcoll[M110] + RMcoll[M210] + RMcoll[M120] + RMcoll[M220])
-
- tmp5 = sp.Rational(1, 4) * rho * (
- RMcoll[M101] + RMcoll[M201] + RMcoll[M102] + RMcoll[M202])
-
- tmp6 = sp.Rational(1, 4) * rho * (
- RMcoll[M011] + RMcoll[M021] + RMcoll[M012] + RMcoll[M022])
-
- main = [
- Assignment(d[F000], rho * (
- 1.0 - RMcoll[M200] - RMcoll[M020] - RMcoll[M002] + RMcoll[
- M220] + RMcoll[M202] + RMcoll[M022])),
-
- Assignment(d[FP00], tmp1),
- Assignment(d[FM00], rho * (-u[0] + RMcoll[M120] + RMcoll[M102]) + tmp1),
-
- Assignment(d[F0P0], tmp2),
- Assignment(d[F0M0], rho * (-u[1] + RMcoll[M210] + RMcoll[M012]) + tmp2),
-
- Assignment(d[F00P], tmp3),
- Assignment(d[F00M], rho * (-u[2] + RMcoll[M201] + RMcoll[M021]) + tmp3),
-
- Assignment(d[FPP0], tmp4),
- Assignment(d[FMP0], sp.Rational(1, 2) * rho * (-RMcoll[M110] - RMcoll[M120]) + tmp4),
- Assignment(d[FPM0], sp.Rational(1, 2) * rho * (-RMcoll[M110] - RMcoll[M210]) + tmp4),
- Assignment(d[FMM0], sp.Rational(1, 2) * rho * (- RMcoll[M210] - RMcoll[M120]) + tmp4),
-
- Assignment(d[FP0P], tmp5),
- Assignment(d[FM0P], sp.Rational(1, 2) * rho * (-RMcoll[M101] - RMcoll[M102]) + tmp5),
- Assignment(d[FP0M], sp.Rational(1, 2) * rho * (-RMcoll[M101] - RMcoll[M201]) + tmp5),
- Assignment(d[FM0M], sp.Rational(1, 2) * rho * (- RMcoll[M201] - RMcoll[M102]) + tmp5),
+ col = [
- Assignment(d[F0PP], tmp6),
- Assignment(d[F0MP], sp.Rational(1, 2) * rho * (-RMcoll[M011] - RMcoll[M012]) + tmp6),
- Assignment(d[F0PM], sp.Rational(1, 2) * rho * (-RMcoll[M011] - RMcoll[M021]) + tmp6),
- Assignment(d[F0MM], sp.Rational(1, 2) * rho * (- RMcoll[M021] - RMcoll[M012]) + tmp6),
+ Assignment(a1xx, p1ij[0]),
+ Assignment(a1xy, p1ij[1]),
+ Assignment(a1xz, p1ij[2]),
+ Assignment(a1yy, p1ij[3]),
+ Assignment(a1yz, p1ij[4]),
+ Assignment(a1zz, p1ij[5]),
+
+ Assignment(a1xxy, sp.Float(2.0) * ux * a1xy + uy * a1xx),
+ Assignment(a1xxz, sp.Float(2.0) * ux * a1xz + uz * a1xx),
+ Assignment(a1xyy, sp.Float(2.0) * uy * a1xy + ux * a1yy),
+ Assignment(a1xzz, sp.Float(2.0) * uz * a1xz + ux * a1zz),
+ Assignment(a1yyz, sp.Float(2.0) * uy * a1yz + uz * a1yy),
+ Assignment(a1yzz, sp.Float(2.0) * uz * a1yz + uy * a1zz),
+
+ Assignment(fNEq[0], -(a1zz + a1yy + a1xx) * sp.Float(0.5)),
+ Assignment(fNEq[1],
+ (sp.Float(2.0) * a1zz - sp.Float(3.0) * a1yyz - a1yy - sp.Float(3.0) * a1xxz - a1xx) * (inv12)),
+ Assignment(fNEq[2],
+ (sp.Float(2.0) * a1zz + sp.Float(3.0) * a1yyz - a1yy + sp.Float(3.0) * a1xxz - a1xx) * (inv12)),
+ Assignment(fNEq[3],
+ -(a1zz + sp.Float(3.0) * a1yzz - sp.Float(2.0) * a1yy + sp.Float(3.0) * a1xxy + a1xx) * (inv12)),
+ Assignment(fNEq[4],
+ -(a1zz - sp.Float(3.0) * a1yzz - sp.Float(2.0) * a1yy - sp.Float(3.0) * a1xxy + a1xx) * (inv12)),
+ Assignment(fNEq[5],
+ -(a1zz + a1yy + sp.Float(3.0) * a1xzz + sp.Float(3.0) * a1xyy - sp.Float(2.0) * a1xx) * (inv12)),
+ Assignment(fNEq[6],
+ -(a1zz + a1yy - sp.Float(3.0) * a1xzz - sp.Float(3.0) * a1xyy - sp.Float(2.0) * a1xx) * (inv12)),
+ Assignment(fNEq[7], (sp.Float(2.0) * a1zz + sp.Float(6.0) * a1yzz + sp.Float(6.0) * a1yz + sp.Float(
+ 6.0) * a1yyz + sp.Float(2.0) * a1yy - sp.Float(3.0) * a1xxz - sp.Float(3.0) * a1xxy - a1xx) * (inv24)),
+ Assignment(fNEq[8], (sp.Float(2.0) * a1zz + sp.Float(6.0) * a1yzz - sp.Float(6.0) * a1yz - sp.Float(
+ 6.0) * a1yyz + sp.Float(2.0) * a1yy + sp.Float(3.0) * a1xxz - sp.Float(3.0) * a1xxy - a1xx) * (inv24)),
+ Assignment(fNEq[9], (sp.Float(2.0) * a1zz - sp.Float(6.0) * a1yzz - sp.Float(6.0) * a1yz + sp.Float(
+ 6.0) * a1yyz + sp.Float(2.0) * a1yy - sp.Float(3.0) * a1xxz + sp.Float(3.0) * a1xxy - a1xx) * (inv24)),
+ Assignment(fNEq[10], (sp.Float(2.0) * a1zz - sp.Float(6.0) * a1yzz + sp.Float(6.0) * a1yz - sp.Float(
+ 6.0) * a1yyz + sp.Float(2.0) * a1yy + sp.Float(3.0) * a1xxz + sp.Float(3.0) * a1xxy - a1xx) * (inv24)),
+ Assignment(fNEq[11], (
+ sp.Float(2.0) * a1zz - sp.Float(3.0) * a1yyz - a1yy + sp.Float(6.0) * a1xzz + sp.Float(
+ 6.0) * a1xz - sp.Float(3.0) * a1xyy + sp.Float(6.0) * a1xxz + sp.Float(2.0) * a1xx) * (inv24)),
+ Assignment(fNEq[12], (
+ sp.Float(2.0) * a1zz + sp.Float(3.0) * a1yyz - a1yy + sp.Float(6.0) * a1xzz - sp.Float(
+ 6.0) * a1xz - sp.Float(3.0) * a1xyy - sp.Float(6.0) * a1xxz + sp.Float(2.0) * a1xx) * (inv24)),
+ Assignment(fNEq[13], (
+ sp.Float(2.0) * a1zz - sp.Float(3.0) * a1yyz - a1yy - sp.Float(6.0) * a1xzz - sp.Float(
+ 6.0) * a1xz + sp.Float(3.0) * a1xyy + sp.Float(6.0) * a1xxz + sp.Float(2.0) * a1xx) * (inv24)),
+ Assignment(fNEq[14], (
+ sp.Float(2.0) * a1zz + sp.Float(3.0) * a1yyz - a1yy - sp.Float(6.0) * a1xzz + sp.Float(
+ 6.0) * a1xz + sp.Float(3.0) * a1xyy - sp.Float(6.0) * a1xxz + sp.Float(2.0) * a1xx) * (inv24)),
+ Assignment(fNEq[15], -(
+ a1zz + sp.Float(3.0) * a1yzz - sp.Float(2.0) * a1yy + sp.Float(3.0) * a1xzz - sp.Float(
+ 6.0) * a1xyy - sp.Float(6.0) * a1xy - sp.Float(6.0) * a1xxy - sp.Float(2.0) * a1xx) * (inv24)),
+ Assignment(fNEq[16], -(
+ a1zz - sp.Float(3.0) * a1yzz - sp.Float(2.0) * a1yy + sp.Float(3.0) * a1xzz - sp.Float(
+ 6.0) * a1xyy + sp.Float(6.0) * a1xy + sp.Float(6.0) * a1xxy - sp.Float(2.0) * a1xx) * (inv24)),
+ Assignment(fNEq[17], -(
+ a1zz + sp.Float(3.0) * a1yzz - sp.Float(2.0) * a1yy - sp.Float(3.0) * a1xzz + sp.Float(
+ 6.0) * a1xyy + sp.Float(6.0) * a1xy - sp.Float(6.0) * a1xxy - sp.Float(2.0) * a1xx) * (inv24)),
+ Assignment(fNEq[18], -(
+ a1zz - sp.Float(3.0) * a1yzz - sp.Float(2.0) * a1yy - sp.Float(3.0) * a1xzz + sp.Float(
+ 6.0) * a1xyy - sp.Float(6.0) * a1xy + sp.Float(6.0) * a1xxy - sp.Float(2.0) * a1xx) * (inv24))
]
- subexpressions = list(additional_subexpressions) + gradient + sr + cqe.all_assignments + moments + eq_moments + neq + col
+ main = []
+ for i in range(19):
+ main.append(Assignment(d[i], fEq[i] + (1 - omega1) * fNEq[i]))
+
+ subexpressions = list(additional_subexpressions) + gradient + sr + cqe.all_assignments + moments + neq + col
ac = AssignmentCollection(subexpressions=subexpressions, main_assignments=main)
ac = ac.new_without_unused_subexpressions()
diff --git a/lbmpy/stencils.py b/lbmpy/stencils.py
index be935bef3f7ba35f2ad9d02208e32449dd8477eb..d122acc449b4fa9826ef932db5a1b417367cadd5 100644
--- a/lbmpy/stencils.py
+++ b/lbmpy/stencils.py
@@ -229,7 +229,41 @@ def _predefined_stencils(stencil: str, ordering: str):
(-1, 0, 1), (0, -1, -1), (0, -1, 1),
(1, 0, 0), (0, 1, 0), (0, 0, 1),
(1, 1, 0), (1, -1, 0), (1, 0, 1),
- (1, 0, -1), (0, 1, 1), (0, 1, -1))
+ (1, 0, -1), (0, 1, 1), (0, 1, -1)),
+
+ 'prolb': ((0, 0, 0), (0, 0, 1),(0, 0,-1),
+
+ (0, 1, 0),
+
+ (0,-1, 0),
+
+ (1, 0, 0),
+
+ (-1, 0, 0),
+
+ (0, 1, 1),
+
+ (0, 1,-1),
+
+ (0,-1, 1),
+
+ (0,-1,-1),
+
+ (1, 0, 1),
+
+ (1, 0,-1),
+
+ (-1, 0, 1),
+
+ (-1, 0,-1),
+
+ (1, 1, 0),
+
+ (1,-1, 0),
+
+ (-1, 1, 0),
+
+ (-1,-1, 0),)
},
'D3Q27': {
'walberla': ((0, 0, 0),