diff --git a/doc/notebooks/demo_central_moments.ipynb b/doc/notebooks/demo_central_moments.ipynb
index cc1a9490a504eb54217c4348bd35896ca77f44cf..30e5cf06aa3148ac506ead04b9eed369eba60734 100644
--- a/doc/notebooks/demo_central_moments.ipynb
+++ b/doc/notebooks/demo_central_moments.ipynb
@@ -255,7 +255,7 @@
" "
],
"text/plain": [
- "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f3c14ff4400>"
+ "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f61e2aea370>"
]
},
"execution_count": 3,
@@ -300,54 +300,54 @@
" <tr style=\"border:none\">\n",
" <td style=\"border:none\">$1$</td>\n",
" <td style=\"border:none\">$\\rho$</td>\n",
- " <td style=\"border:none\">$1.99$</td>\n",
+ " <td style=\"border:none\">$rr$</td>\n",
" </tr>\n",
"<tr style=\"border:none\">\n",
" <td style=\"border:none\">$x$</td>\n",
" <td style=\"border:none\">$\\rho u_{0}$</td>\n",
- " <td style=\"border:none\">$1.99$</td>\n",
+ " <td style=\"border:none\">$rr$</td>\n",
" </tr>\n",
"<tr style=\"border:none\">\n",
" <td style=\"border:none\">$y$</td>\n",
" <td style=\"border:none\">$\\rho u_{1}$</td>\n",
- " <td style=\"border:none\">$1.99$</td>\n",
+ " <td style=\"border:none\">$rr$</td>\n",
" </tr>\n",
"<tr style=\"border:none\">\n",
" <td style=\"border:none\">$x^{2}$</td>\n",
" <td style=\"border:none\">$\\rho u_{0}^{2} + \\frac{\\rho}{3}$</td>\n",
- " <td style=\"border:none\">$1.99$</td>\n",
+ " <td style=\"border:none\">$rr$</td>\n",
" </tr>\n",
"<tr style=\"border:none\">\n",
" <td style=\"border:none\">$y^{2}$</td>\n",
" <td style=\"border:none\">$\\rho u_{1}^{2} + \\frac{\\rho}{3}$</td>\n",
- " <td style=\"border:none\">$1.99$</td>\n",
+ " <td style=\"border:none\">$rr$</td>\n",
" </tr>\n",
"<tr style=\"border:none\">\n",
" <td style=\"border:none\">$x y$</td>\n",
" <td style=\"border:none\">$\\rho u_{0} u_{1}$</td>\n",
- " <td style=\"border:none\">$1.99$</td>\n",
+ " <td style=\"border:none\">$rr$</td>\n",
" </tr>\n",
"<tr style=\"border:none\">\n",
" <td style=\"border:none\">$x^{2} y$</td>\n",
" <td style=\"border:none\">$\\rho u_{0}^{2} u_{1} + \\frac{\\rho u_{1}}{3}$</td>\n",
- " <td style=\"border:none\">$1.99$</td>\n",
+ " <td style=\"border:none\">$rr$</td>\n",
" </tr>\n",
"<tr style=\"border:none\">\n",
" <td style=\"border:none\">$x y^{2}$</td>\n",
" <td style=\"border:none\">$\\rho u_{0} u_{1}^{2} + \\frac{\\rho u_{0}}{3}$</td>\n",
- " <td style=\"border:none\">$1.99$</td>\n",
+ " <td style=\"border:none\">$rr$</td>\n",
" </tr>\n",
"<tr style=\"border:none\">\n",
" <td style=\"border:none\">$x^{2} y^{2}$</td>\n",
" <td style=\"border:none\">$\\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}$</td>\n",
- " <td style=\"border:none\">$1.99$</td>\n",
+ " <td style=\"border:none\">$rr$</td>\n",
" </tr>\n",
"\n",
" </table>\n",
" "
],
"text/plain": [
- "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f3c14fbfdc0>"
+ "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f62081143a0>"
]
},
"execution_count": 4,
@@ -357,7 +357,7 @@
],
"source": [
"def modification_func(moment, eq, rr):\n",
- " rr = 1.99\n",
+ " rr = sp.Symbol('rr')\n",
" return moment, eq, rr\n",
"\n",
"method = create_lb_method_from_existing(method, modification_func)\n",
@@ -420,7 +420,7 @@
"S = sp.diag(*method.relaxation_rates)\n",
"c_eq = sp.Matrix(method.moment_equilibrium_values)\n",
"\n",
- "N = shift_matrix(raw_moments, stencil)\n",
+ "N = set_up_shift_matrix(raw_moments, stencil)\n",
"\n",
"N"
]
@@ -438,7 +438,7 @@
"metadata": {},
"outputs": [],
"source": [
- "collision_equation = f + C.inv() * N.inv() * S * N * (c_eq - C * f) + sp.Matrix(F(method))"
+ "collision_equation = f + C.inv() * N.inv() * S * N * (c_eq - C * f)# + sp.Matrix(F(method))"
]
},
{
@@ -450,7 +450,7 @@
},
{
"cell_type": "code",
- "execution_count": 7,
+ "execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
@@ -459,14 +459,80 @@
"@ps.kernel\n",
"def subexpressions():\n",
" ρ @= sum(f) \n",
- " u[0] @= sum(f.T * sp.Matrix(np.array(stencil)[:, 0]))/ρ + force[0]/(2*ρ)\n",
- " u[1] @= sum(f.T * sp.Matrix(np.array(stencil)[:, 1]))/ρ + force[1]/(2*ρ)\n",
+ " u[0] @= sum(f.T * sp.Matrix(np.array(stencil)[:, 0])) /ρ # + force[0]/(2*ρ)\n",
+ " u[1] @= sum(f.T * sp.Matrix(np.array(stencil)[:, 1])) /ρ # + force[1]/(2*ρ)\n",
"\n",
"collision_rule = LbmCollisionRule(method, main_assignments=collision_rule, subexpressions=[subexpressions[0],\n",
" subexpressions[1],\n",
" subexpressions[2]])"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<div>Subexpressions:</div><table style=\"border:none; width: 100%; \"><tr style=\"border:none\"> <td style=\"border:none\">$$\\rho \\leftarrow f_{0} + f_{1} + f_{2} + f_{3} + f_{4} + f_{5} + f_{6} + f_{7} + f_{8}$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$u_{0} \\leftarrow \\frac{- f_{3} + f_{4} - f_{5} + f_{6} - f_{7} + f_{8}}{\\rho}$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$u_{1} \\leftarrow \\frac{f_{1} - f_{2} + f_{5} + f_{6} - f_{7} - f_{8}}{\\rho}$$</td> </tr> </table><div>Main Assignments:</div><table style=\"border:none; width: 100%; \"><tr style=\"border:none\"> <td style=\"border:none\">$$d_{0} \\leftarrow f_{0} + rr \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right) + \\left(- rr u_{0}^{2} + rr \\left(u_{0}^{2} - 1\\right)\\right) \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right) + \\left(- rr u_{1}^{2} + rr \\left(u_{1}^{2} - 1\\right)\\right) \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right) + \\left(- 2 rr u_{0} \\left(u_{1}^{2} - 1\\right) + rr \\left(2 u_{0} u_{1}^{2} - 2 u_{0}\\right)\\right) \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right) + \\left(- 2 rr u_{1} \\left(u_{0}^{2} - 1\\right) + rr \\left(2 u_{0}^{2} u_{1} - 2 u_{1}\\right)\\right) \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right) + \\left(rr u_{0}^{2} u_{1}^{2} + rr u_{0}^{2} \\left(u_{1}^{2} - 1\\right) - rr u_{0} \\left(2 u_{0} u_{1}^{2} - 2 u_{0}\\right) + rr u_{1}^{2} \\left(u_{0}^{2} - 1\\right) - rr u_{1} \\left(2 u_{0}^{2} u_{1} - 2 u_{1}\\right) + rr \\left(u_{0}^{2} u_{1}^{2} - u_{0}^{2} - u_{1}^{2} + 1\\right)\\right) \\left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{1} \\leftarrow f_{1} - \\frac{rr \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\left(\\frac{rr u_{0}^{2}}{2} + rr \\left(\\frac{1}{2} - \\frac{u_{0}^{2}}{2}\\right)\\right) \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right) + \\left(rr u_{1} + rr \\left(- u_{1} - \\frac{1}{2}\\right)\\right) \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{0}^{2} u_{1} + \\frac{\\rho u_{1}}{3}\\right) + \\left(- 2 rr u_{0} \\left(- u_{1} - \\frac{1}{2}\\right) + rr \\left(- 2 u_{0} u_{1} - u_{0}\\right)\\right) \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right) + \\left(- \\frac{rr u_{1}^{2}}{2} - rr u_{1} \\left(- u_{1} - \\frac{1}{2}\\right) + rr \\left(- \\frac{u_{1}^{2}}{2} - \\frac{u_{1}}{2}\\right)\\right) \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right) + \\left(2 rr u_{0} u_{1} \\left(- u_{1} - \\frac{1}{2}\\right) - 2 rr u_{0} \\left(- \\frac{u_{1}^{2}}{2} - \\frac{u_{1}}{2}\\right) - rr u_{1} \\left(- 2 u_{0} u_{1} - u_{0}\\right) + rr \\left(- u_{0} u_{1}^{2} - u_{0} u_{1}\\right)\\right) \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right) + \\left(- rr u_{0}^{2} u_{1} + rr u_{0}^{2} \\left(- u_{1} - \\frac{1}{2}\\right) - rr u_{0} \\left(- 2 u_{0} u_{1} - u_{0}\\right) - 2 rr u_{1} \\left(\\frac{1}{2} - \\frac{u_{0}^{2}}{2}\\right) + rr \\left(- u_{0}^{2} u_{1} - \\frac{u_{0}^{2}}{2} + u_{1} + \\frac{1}{2}\\right)\\right) \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right) + \\left(\\frac{rr u_{0}^{2} u_{1}^{2}}{2} - rr u_{0}^{2} u_{1} \\left(- u_{1} - \\frac{1}{2}\\right) + rr u_{0}^{2} \\left(- \\frac{u_{1}^{2}}{2} - \\frac{u_{1}}{2}\\right) + rr u_{0} u_{1} \\left(- 2 u_{0} u_{1} - u_{0}\\right) - rr u_{0} \\left(- u_{0} u_{1}^{2} - u_{0} u_{1}\\right) + rr u_{1}^{2} \\left(\\frac{1}{2} - \\frac{u_{0}^{2}}{2}\\right) - rr u_{1} \\left(- u_{0}^{2} u_{1} - \\frac{u_{0}^{2}}{2} + u_{1} + \\frac{1}{2}\\right) + rr \\left(- \\frac{u_{0}^{2} u_{1}^{2}}{2} - \\frac{u_{0}^{2} u_{1}}{2} + \\frac{u_{1}^{2}}{2} + \\frac{u_{1}}{2}\\right)\\right) \\left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{2} \\leftarrow f_{2} - \\frac{rr \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\left(\\frac{rr u_{0}^{2}}{2} + rr \\left(\\frac{1}{2} - \\frac{u_{0}^{2}}{2}\\right)\\right) \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right) + \\left(rr u_{1} + rr \\left(\\frac{1}{2} - u_{1}\\right)\\right) \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{0}^{2} u_{1} + \\frac{\\rho u_{1}}{3}\\right) + \\left(- 2 rr u_{0} \\left(\\frac{1}{2} - u_{1}\\right) + rr \\left(- 2 u_{0} u_{1} + u_{0}\\right)\\right) \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right) + \\left(- \\frac{rr u_{1}^{2}}{2} - rr u_{1} \\left(\\frac{1}{2} - u_{1}\\right) + rr \\left(- \\frac{u_{1}^{2}}{2} + \\frac{u_{1}}{2}\\right)\\right) \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right) + \\left(2 rr u_{0} u_{1} \\left(\\frac{1}{2} - u_{1}\\right) - 2 rr u_{0} \\left(- \\frac{u_{1}^{2}}{2} + \\frac{u_{1}}{2}\\right) - rr u_{1} \\left(- 2 u_{0} u_{1} + u_{0}\\right) + rr \\left(- u_{0} u_{1}^{2} + u_{0} u_{1}\\right)\\right) \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right) + \\left(- rr u_{0}^{2} u_{1} + rr u_{0}^{2} \\left(\\frac{1}{2} - u_{1}\\right) - rr u_{0} \\left(- 2 u_{0} u_{1} + u_{0}\\right) - 2 rr u_{1} \\left(\\frac{1}{2} - \\frac{u_{0}^{2}}{2}\\right) + rr \\left(- u_{0}^{2} u_{1} + \\frac{u_{0}^{2}}{2} + u_{1} - \\frac{1}{2}\\right)\\right) \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right) + \\left(\\frac{rr u_{0}^{2} u_{1}^{2}}{2} - rr u_{0}^{2} u_{1} \\left(\\frac{1}{2} - u_{1}\\right) + rr u_{0}^{2} \\left(- \\frac{u_{1}^{2}}{2} + \\frac{u_{1}}{2}\\right) + rr u_{0} u_{1} \\left(- 2 u_{0} u_{1} + u_{0}\\right) - rr u_{0} \\left(- u_{0} u_{1}^{2} + u_{0} u_{1}\\right) + rr u_{1}^{2} \\left(\\frac{1}{2} - \\frac{u_{0}^{2}}{2}\\right) - rr u_{1} \\left(- u_{0}^{2} u_{1} + \\frac{u_{0}^{2}}{2} + u_{1} - \\frac{1}{2}\\right) + rr \\left(- \\frac{u_{0}^{2} u_{1}^{2}}{2} + \\frac{u_{0}^{2} u_{1}}{2} + \\frac{u_{1}^{2}}{2} - \\frac{u_{1}}{2}\\right)\\right) \\left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{3} \\leftarrow f_{3} - \\frac{rr \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\left(rr u_{0} + rr \\left(\\frac{1}{2} - u_{0}\\right)\\right) \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0} u_{1}^{2} + \\frac{\\rho u_{0}}{3}\\right) + \\left(\\frac{rr u_{1}^{2}}{2} + rr \\left(\\frac{1}{2} - \\frac{u_{1}^{2}}{2}\\right)\\right) \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right) + \\left(- 2 rr u_{1} \\left(\\frac{1}{2} - u_{0}\\right) + rr \\left(- 2 u_{0} u_{1} + u_{1}\\right)\\right) \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right) + \\left(- \\frac{rr u_{0}^{2}}{2} - rr u_{0} \\left(\\frac{1}{2} - u_{0}\\right) + rr \\left(- \\frac{u_{0}^{2}}{2} + \\frac{u_{0}}{2}\\right)\\right) \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right) + \\left(2 rr u_{0} u_{1} \\left(\\frac{1}{2} - u_{0}\\right) - rr u_{0} \\left(- 2 u_{0} u_{1} + u_{1}\\right) - 2 rr u_{1} \\left(- \\frac{u_{0}^{2}}{2} + \\frac{u_{0}}{2}\\right) + rr \\left(- u_{0}^{2} u_{1} + u_{0} u_{1}\\right)\\right) \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right) + \\left(- rr u_{0} u_{1}^{2} - 2 rr u_{0} \\left(\\frac{1}{2} - \\frac{u_{1}^{2}}{2}\\right) + rr u_{1}^{2} \\left(\\frac{1}{2} - u_{0}\\right) - rr u_{1} \\left(- 2 u_{0} u_{1} + u_{1}\\right) + rr \\left(- u_{0} u_{1}^{2} + u_{0} + \\frac{u_{1}^{2}}{2} - \\frac{1}{2}\\right)\\right) \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right) + \\left(\\frac{rr u_{0}^{2} u_{1}^{2}}{2} + rr u_{0}^{2} \\left(\\frac{1}{2} - \\frac{u_{1}^{2}}{2}\\right) - rr u_{0} u_{1}^{2} \\left(\\frac{1}{2} - u_{0}\\right) + rr u_{0} u_{1} \\left(- 2 u_{0} u_{1} + u_{1}\\right) - rr u_{0} \\left(- u_{0} u_{1}^{2} + u_{0} + \\frac{u_{1}^{2}}{2} - \\frac{1}{2}\\right) + rr u_{1}^{2} \\left(- \\frac{u_{0}^{2}}{2} + \\frac{u_{0}}{2}\\right) - rr u_{1} \\left(- u_{0}^{2} u_{1} + u_{0} u_{1}\\right) + rr \\left(- \\frac{u_{0}^{2} u_{1}^{2}}{2} + \\frac{u_{0}^{2}}{2} + \\frac{u_{0} u_{1}^{2}}{2} - \\frac{u_{0}}{2}\\right)\\right) \\left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{4} \\leftarrow f_{4} - \\frac{rr \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\left(rr u_{0} + rr \\left(- u_{0} - \\frac{1}{2}\\right)\\right) \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0} u_{1}^{2} + \\frac{\\rho u_{0}}{3}\\right) + \\left(\\frac{rr u_{1}^{2}}{2} + rr \\left(\\frac{1}{2} - \\frac{u_{1}^{2}}{2}\\right)\\right) \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right) + \\left(- 2 rr u_{1} \\left(- u_{0} - \\frac{1}{2}\\right) + rr \\left(- 2 u_{0} u_{1} - u_{1}\\right)\\right) \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right) + \\left(- \\frac{rr u_{0}^{2}}{2} - rr u_{0} \\left(- u_{0} - \\frac{1}{2}\\right) + rr \\left(- \\frac{u_{0}^{2}}{2} - \\frac{u_{0}}{2}\\right)\\right) \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right) + \\left(2 rr u_{0} u_{1} \\left(- u_{0} - \\frac{1}{2}\\right) - rr u_{0} \\left(- 2 u_{0} u_{1} - u_{1}\\right) - 2 rr u_{1} \\left(- \\frac{u_{0}^{2}}{2} - \\frac{u_{0}}{2}\\right) + rr \\left(- u_{0}^{2} u_{1} - u_{0} u_{1}\\right)\\right) \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right) + \\left(- rr u_{0} u_{1}^{2} - 2 rr u_{0} \\left(\\frac{1}{2} - \\frac{u_{1}^{2}}{2}\\right) + rr u_{1}^{2} \\left(- u_{0} - \\frac{1}{2}\\right) - rr u_{1} \\left(- 2 u_{0} u_{1} - u_{1}\\right) + rr \\left(- u_{0} u_{1}^{2} + u_{0} - \\frac{u_{1}^{2}}{2} + \\frac{1}{2}\\right)\\right) \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right) + \\left(\\frac{rr u_{0}^{2} u_{1}^{2}}{2} + rr u_{0}^{2} \\left(\\frac{1}{2} - \\frac{u_{1}^{2}}{2}\\right) - rr u_{0} u_{1}^{2} \\left(- u_{0} - \\frac{1}{2}\\right) + rr u_{0} u_{1} \\left(- 2 u_{0} u_{1} - u_{1}\\right) - rr u_{0} \\left(- u_{0} u_{1}^{2} + u_{0} - \\frac{u_{1}^{2}}{2} + \\frac{1}{2}\\right) + rr u_{1}^{2} \\left(- \\frac{u_{0}^{2}}{2} - \\frac{u_{0}}{2}\\right) - rr u_{1} \\left(- u_{0}^{2} u_{1} - u_{0} u_{1}\\right) + rr \\left(- \\frac{u_{0}^{2} u_{1}^{2}}{2} + \\frac{u_{0}^{2}}{2} - \\frac{u_{0} u_{1}^{2}}{2} + \\frac{u_{0}}{2}\\right)\\right) \\left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{5} \\leftarrow f_{5} + \\frac{rr \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4} + \\left(- \\frac{rr u_{0}}{2} + rr \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right)\\right) \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0} u_{1}^{2} + \\frac{\\rho u_{0}}{3}\\right) + \\left(- \\frac{rr u_{1}}{2} + rr \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right)\\right) \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{0}^{2} u_{1} + \\frac{\\rho u_{1}}{3}\\right) + \\left(\\frac{rr u_{0}^{2}}{4} - rr u_{0} \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right) + rr \\left(\\frac{u_{0}^{2}}{4} - \\frac{u_{0}}{4}\\right)\\right) \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right) + \\left(\\frac{rr u_{1}^{2}}{4} - rr u_{1} \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right) + rr \\left(\\frac{u_{1}^{2}}{4} + \\frac{u_{1}}{4}\\right)\\right) \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right) + \\left(rr u_{0} u_{1} - 2 rr u_{0} \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right) - 2 rr u_{1} \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right) + rr \\left(u_{0} u_{1} + \\frac{u_{0}}{2} - \\frac{u_{1}}{2} - \\frac{1}{4}\\right)\\right) \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right) + \\left(- \\frac{rr u_{0} u_{1}^{2}}{2} + 2 rr u_{0} u_{1} \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right) - 2 rr u_{0} \\left(\\frac{u_{1}^{2}}{4} + \\frac{u_{1}}{4}\\right) + rr u_{1}^{2} \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right) - rr u_{1} \\left(u_{0} u_{1} + \\frac{u_{0}}{2} - \\frac{u_{1}}{2} - \\frac{1}{4}\\right) + rr \\left(\\frac{u_{0} u_{1}^{2}}{2} + \\frac{u_{0} u_{1}}{2} - \\frac{u_{1}^{2}}{4} - \\frac{u_{1}}{4}\\right)\\right) \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right) + \\left(- \\frac{rr u_{0}^{2} u_{1}}{2} + rr u_{0}^{2} \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right) + 2 rr u_{0} u_{1} \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right) - rr u_{0} \\left(u_{0} u_{1} + \\frac{u_{0}}{2} - \\frac{u_{1}}{2} - \\frac{1}{4}\\right) - 2 rr u_{1} \\left(\\frac{u_{0}^{2}}{4} - \\frac{u_{0}}{4}\\right) + rr \\left(\\frac{u_{0}^{2} u_{1}}{2} + \\frac{u_{0}^{2}}{4} - \\frac{u_{0} u_{1}}{2} - \\frac{u_{0}}{4}\\right)\\right) \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right) + \\left(\\frac{rr u_{0}^{2} u_{1}^{2}}{4} - rr u_{0}^{2} u_{1} \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right) + rr u_{0}^{2} \\left(\\frac{u_{1}^{2}}{4} + \\frac{u_{1}}{4}\\right) - rr u_{0} u_{1}^{2} \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right) + rr u_{0} u_{1} \\left(u_{0} u_{1} + \\frac{u_{0}}{2} - \\frac{u_{1}}{2} - \\frac{1}{4}\\right) - rr u_{0} \\left(\\frac{u_{0} u_{1}^{2}}{2} + \\frac{u_{0} u_{1}}{2} - \\frac{u_{1}^{2}}{4} - \\frac{u_{1}}{4}\\right) + rr u_{1}^{2} \\left(\\frac{u_{0}^{2}}{4} - \\frac{u_{0}}{4}\\right) - rr u_{1} \\left(\\frac{u_{0}^{2} u_{1}}{2} + \\frac{u_{0}^{2}}{4} - \\frac{u_{0} u_{1}}{2} - \\frac{u_{0}}{4}\\right) + rr \\left(\\frac{u_{0}^{2} u_{1}^{2}}{4} + \\frac{u_{0}^{2} u_{1}}{4} - \\frac{u_{0} u_{1}^{2}}{4} - \\frac{u_{0} u_{1}}{4}\\right)\\right) \\left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{6} \\leftarrow f_{6} + \\frac{rr \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4} + \\left(- \\frac{rr u_{0}}{2} + rr \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right)\\right) \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0} u_{1}^{2} + \\frac{\\rho u_{0}}{3}\\right) + \\left(- \\frac{rr u_{1}}{2} + rr \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right)\\right) \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{0}^{2} u_{1} + \\frac{\\rho u_{1}}{3}\\right) + \\left(\\frac{rr u_{0}^{2}}{4} - rr u_{0} \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right) + rr \\left(\\frac{u_{0}^{2}}{4} + \\frac{u_{0}}{4}\\right)\\right) \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right) + \\left(\\frac{rr u_{1}^{2}}{4} - rr u_{1} \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right) + rr \\left(\\frac{u_{1}^{2}}{4} + \\frac{u_{1}}{4}\\right)\\right) \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right) + \\left(rr u_{0} u_{1} - 2 rr u_{0} \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right) - 2 rr u_{1} \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right) + rr \\left(u_{0} u_{1} + \\frac{u_{0}}{2} + \\frac{u_{1}}{2} + \\frac{1}{4}\\right)\\right) \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right) + \\left(- \\frac{rr u_{0} u_{1}^{2}}{2} + 2 rr u_{0} u_{1} \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right) - 2 rr u_{0} \\left(\\frac{u_{1}^{2}}{4} + \\frac{u_{1}}{4}\\right) + rr u_{1}^{2} \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right) - rr u_{1} \\left(u_{0} u_{1} + \\frac{u_{0}}{2} + \\frac{u_{1}}{2} + \\frac{1}{4}\\right) + rr \\left(\\frac{u_{0} u_{1}^{2}}{2} + \\frac{u_{0} u_{1}}{2} + \\frac{u_{1}^{2}}{4} + \\frac{u_{1}}{4}\\right)\\right) \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right) + \\left(- \\frac{rr u_{0}^{2} u_{1}}{2} + rr u_{0}^{2} \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right) + 2 rr u_{0} u_{1} \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right) - rr u_{0} \\left(u_{0} u_{1} + \\frac{u_{0}}{2} + \\frac{u_{1}}{2} + \\frac{1}{4}\\right) - 2 rr u_{1} \\left(\\frac{u_{0}^{2}}{4} + \\frac{u_{0}}{4}\\right) + rr \\left(\\frac{u_{0}^{2} u_{1}}{2} + \\frac{u_{0}^{2}}{4} + \\frac{u_{0} u_{1}}{2} + \\frac{u_{0}}{4}\\right)\\right) \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right) + \\left(\\frac{rr u_{0}^{2} u_{1}^{2}}{4} - rr u_{0}^{2} u_{1} \\left(\\frac{u_{1}}{2} + \\frac{1}{4}\\right) + rr u_{0}^{2} \\left(\\frac{u_{1}^{2}}{4} + \\frac{u_{1}}{4}\\right) - rr u_{0} u_{1}^{2} \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right) + rr u_{0} u_{1} \\left(u_{0} u_{1} + \\frac{u_{0}}{2} + \\frac{u_{1}}{2} + \\frac{1}{4}\\right) - rr u_{0} \\left(\\frac{u_{0} u_{1}^{2}}{2} + \\frac{u_{0} u_{1}}{2} + \\frac{u_{1}^{2}}{4} + \\frac{u_{1}}{4}\\right) + rr u_{1}^{2} \\left(\\frac{u_{0}^{2}}{4} + \\frac{u_{0}}{4}\\right) - rr u_{1} \\left(\\frac{u_{0}^{2} u_{1}}{2} + \\frac{u_{0}^{2}}{4} + \\frac{u_{0} u_{1}}{2} + \\frac{u_{0}}{4}\\right) + rr \\left(\\frac{u_{0}^{2} u_{1}^{2}}{4} + \\frac{u_{0}^{2} u_{1}}{4} + \\frac{u_{0} u_{1}^{2}}{4} + \\frac{u_{0} u_{1}}{4}\\right)\\right) \\left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{7} \\leftarrow f_{7} + \\frac{rr \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4} + \\left(- \\frac{rr u_{0}}{2} + rr \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right)\\right) \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0} u_{1}^{2} + \\frac{\\rho u_{0}}{3}\\right) + \\left(- \\frac{rr u_{1}}{2} + rr \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right)\\right) \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{0}^{2} u_{1} + \\frac{\\rho u_{1}}{3}\\right) + \\left(\\frac{rr u_{0}^{2}}{4} - rr u_{0} \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right) + rr \\left(\\frac{u_{0}^{2}}{4} - \\frac{u_{0}}{4}\\right)\\right) \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right) + \\left(\\frac{rr u_{1}^{2}}{4} - rr u_{1} \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right) + rr \\left(\\frac{u_{1}^{2}}{4} - \\frac{u_{1}}{4}\\right)\\right) \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right) + \\left(rr u_{0} u_{1} - 2 rr u_{0} \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right) - 2 rr u_{1} \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right) + rr \\left(u_{0} u_{1} - \\frac{u_{0}}{2} - \\frac{u_{1}}{2} + \\frac{1}{4}\\right)\\right) \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right) + \\left(- \\frac{rr u_{0} u_{1}^{2}}{2} + 2 rr u_{0} u_{1} \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right) - 2 rr u_{0} \\left(\\frac{u_{1}^{2}}{4} - \\frac{u_{1}}{4}\\right) + rr u_{1}^{2} \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right) - rr u_{1} \\left(u_{0} u_{1} - \\frac{u_{0}}{2} - \\frac{u_{1}}{2} + \\frac{1}{4}\\right) + rr \\left(\\frac{u_{0} u_{1}^{2}}{2} - \\frac{u_{0} u_{1}}{2} - \\frac{u_{1}^{2}}{4} + \\frac{u_{1}}{4}\\right)\\right) \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right) + \\left(- \\frac{rr u_{0}^{2} u_{1}}{2} + rr u_{0}^{2} \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right) + 2 rr u_{0} u_{1} \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right) - rr u_{0} \\left(u_{0} u_{1} - \\frac{u_{0}}{2} - \\frac{u_{1}}{2} + \\frac{1}{4}\\right) - 2 rr u_{1} \\left(\\frac{u_{0}^{2}}{4} - \\frac{u_{0}}{4}\\right) + rr \\left(\\frac{u_{0}^{2} u_{1}}{2} - \\frac{u_{0}^{2}}{4} - \\frac{u_{0} u_{1}}{2} + \\frac{u_{0}}{4}\\right)\\right) \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right) + \\left(\\frac{rr u_{0}^{2} u_{1}^{2}}{4} - rr u_{0}^{2} u_{1} \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right) + rr u_{0}^{2} \\left(\\frac{u_{1}^{2}}{4} - \\frac{u_{1}}{4}\\right) - rr u_{0} u_{1}^{2} \\left(\\frac{u_{0}}{2} - \\frac{1}{4}\\right) + rr u_{0} u_{1} \\left(u_{0} u_{1} - \\frac{u_{0}}{2} - \\frac{u_{1}}{2} + \\frac{1}{4}\\right) - rr u_{0} \\left(\\frac{u_{0} u_{1}^{2}}{2} - \\frac{u_{0} u_{1}}{2} - \\frac{u_{1}^{2}}{4} + \\frac{u_{1}}{4}\\right) + rr u_{1}^{2} \\left(\\frac{u_{0}^{2}}{4} - \\frac{u_{0}}{4}\\right) - rr u_{1} \\left(\\frac{u_{0}^{2} u_{1}}{2} - \\frac{u_{0}^{2}}{4} - \\frac{u_{0} u_{1}}{2} + \\frac{u_{0}}{4}\\right) + rr \\left(\\frac{u_{0}^{2} u_{1}^{2}}{4} - \\frac{u_{0}^{2} u_{1}}{4} - \\frac{u_{0} u_{1}^{2}}{4} + \\frac{u_{0} u_{1}}{4}\\right)\\right) \\left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{8} \\leftarrow f_{8} + \\frac{rr \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4} + \\left(- \\frac{rr u_{0}}{2} + rr \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right)\\right) \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0} u_{1}^{2} + \\frac{\\rho u_{0}}{3}\\right) + \\left(- \\frac{rr u_{1}}{2} + rr \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right)\\right) \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{0}^{2} u_{1} + \\frac{\\rho u_{1}}{3}\\right) + \\left(\\frac{rr u_{0}^{2}}{4} - rr u_{0} \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right) + rr \\left(\\frac{u_{0}^{2}}{4} + \\frac{u_{0}}{4}\\right)\\right) \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right) + \\left(\\frac{rr u_{1}^{2}}{4} - rr u_{1} \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right) + rr \\left(\\frac{u_{1}^{2}}{4} - \\frac{u_{1}}{4}\\right)\\right) \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right) + \\left(rr u_{0} u_{1} - 2 rr u_{0} \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right) - 2 rr u_{1} \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right) + rr \\left(u_{0} u_{1} - \\frac{u_{0}}{2} + \\frac{u_{1}}{2} - \\frac{1}{4}\\right)\\right) \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right) + \\left(- \\frac{rr u_{0} u_{1}^{2}}{2} + 2 rr u_{0} u_{1} \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right) - 2 rr u_{0} \\left(\\frac{u_{1}^{2}}{4} - \\frac{u_{1}}{4}\\right) + rr u_{1}^{2} \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right) - rr u_{1} \\left(u_{0} u_{1} - \\frac{u_{0}}{2} + \\frac{u_{1}}{2} - \\frac{1}{4}\\right) + rr \\left(\\frac{u_{0} u_{1}^{2}}{2} - \\frac{u_{0} u_{1}}{2} + \\frac{u_{1}^{2}}{4} - \\frac{u_{1}}{4}\\right)\\right) \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right) + \\left(- \\frac{rr u_{0}^{2} u_{1}}{2} + rr u_{0}^{2} \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right) + 2 rr u_{0} u_{1} \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right) - rr u_{0} \\left(u_{0} u_{1} - \\frac{u_{0}}{2} + \\frac{u_{1}}{2} - \\frac{1}{4}\\right) - 2 rr u_{1} \\left(\\frac{u_{0}^{2}}{4} + \\frac{u_{0}}{4}\\right) + rr \\left(\\frac{u_{0}^{2} u_{1}}{2} - \\frac{u_{0}^{2}}{4} + \\frac{u_{0} u_{1}}{2} - \\frac{u_{0}}{4}\\right)\\right) \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right) + \\left(\\frac{rr u_{0}^{2} u_{1}^{2}}{4} - rr u_{0}^{2} u_{1} \\left(\\frac{u_{1}}{2} - \\frac{1}{4}\\right) + rr u_{0}^{2} \\left(\\frac{u_{1}^{2}}{4} - \\frac{u_{1}}{4}\\right) - rr u_{0} u_{1}^{2} \\left(\\frac{u_{0}}{2} + \\frac{1}{4}\\right) + rr u_{0} u_{1} \\left(u_{0} u_{1} - \\frac{u_{0}}{2} + \\frac{u_{1}}{2} - \\frac{1}{4}\\right) - rr u_{0} \\left(\\frac{u_{0} u_{1}^{2}}{2} - \\frac{u_{0} u_{1}}{2} + \\frac{u_{1}^{2}}{4} - \\frac{u_{1}}{4}\\right) + rr u_{1}^{2} \\left(\\frac{u_{0}^{2}}{4} + \\frac{u_{0}}{4}\\right) - rr u_{1} \\left(\\frac{u_{0}^{2} u_{1}}{2} - \\frac{u_{0}^{2}}{4} + \\frac{u_{0} u_{1}}{2} - \\frac{u_{0}}{4}\\right) + rr \\left(\\frac{u_{0}^{2} u_{1}^{2}}{4} - \\frac{u_{0}^{2} u_{1}}{4} + \\frac{u_{0} u_{1}^{2}}{4} - \\frac{u_{0} u_{1}}{4}\\right)\\right) \\left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho\\right)$$</td> </tr> </table>"
+ ],
+ "text/plain": [
+ "AssignmentCollection: d_0, d_5, d_6, d_3, d_4, d_8, d_2, d_7, d_1 <- f(f_0, f_3, f_4, f_5, f_2, f_8, f_7, rr, f_1, f_6)"
+ ]
+ },
+ "execution_count": 15,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "collision_rule"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "collision_rule.simplification_hints['density'] = ρ\n",
+ "collision_rule.simplification_hints['velocity'] = (u[0], u[1])\n",
+ "collision_rule.simplification_hints['relaxation_rates'] = [sp.Symbol('rr')]*9"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "simplification = create_simplification_strategy(method, split_inner_loop=False)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<div>Subexpressions:</div><table style=\"border:none; width: 100%; \"><tr style=\"border:none\"> <td style=\"border:none\">$$\\rho \\leftarrow f_{0} + f_{1} + f_{2} + f_{3} + f_{4} + f_{5} + f_{6} + f_{7} + f_{8}$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$\\xi_{0} \\leftarrow \\frac{1}{\\rho}$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$u_{0} \\leftarrow \\xi_{0} \\left(- f_{3} + f_{4} - f_{5} + f_{6} - f_{7} + f_{8}\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$u_{1} \\leftarrow \\xi_{0} \\left(f_{1} - f_{2} + f_{5} + f_{6} - f_{7} - f_{8}\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$u0Mu1 \\leftarrow u_{0} - u_{1}$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$u0Pu1 \\leftarrow u_{0} + u_{1}$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$f_{eq common} \\leftarrow \\frac{9 \\rho u_{0}^{2} u_{1}^{2}}{4} - \\frac{3 \\rho u_{0}^{2}}{2} - \\frac{3 \\rho u_{1}^{2}}{2} + \\rho$$</td> </tr> </table><div>Main Assignments:</div><table style=\"border:none; width: 100%; \"><tr style=\"border:none\"> <td style=\"border:none\">$$d_{0} \\leftarrow f_{0} + rr \\left(- f_{0} + \\frac{4 f_{eq common}}{9}\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{1} \\leftarrow f_{1} + rr \\left(- f_{1} + \\frac{f_{eq common}}{9} + \\rho \\left(- \\frac{3 u_{0}^{2} u_{1}^{2}}{4} - \\frac{u_{0}^{2} u_{1}}{2} + \\frac{u_{1}^{2}}{2} + \\frac{u_{1}}{3}\\right)\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{2} \\leftarrow f_{2} + rr \\left(- f_{2} + \\frac{f_{eq common}}{9} + \\rho \\left(- \\frac{3 u_{0}^{2} u_{1}^{2}}{4} + \\frac{u_{0}^{2} u_{1}}{2} + \\frac{u_{1}^{2}}{2} - \\frac{u_{1}}{3}\\right)\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{3} \\leftarrow f_{3} + rr \\left(- f_{3} + \\frac{f_{eq common}}{9} + \\rho \\left(- \\frac{3 u_{0}^{2} u_{1}^{2}}{4} + \\frac{u_{0}^{2}}{2} + \\frac{u_{0} u_{1}^{2}}{2} - \\frac{u_{0}}{3}\\right)\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{4} \\leftarrow f_{4} + rr \\left(- f_{4} + \\frac{f_{eq common}}{9} + \\rho \\left(- \\frac{3 u_{0}^{2} u_{1}^{2}}{4} + \\frac{u_{0}^{2}}{2} - \\frac{u_{0} u_{1}^{2}}{2} + \\frac{u_{0}}{3}\\right)\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{5} \\leftarrow f_{5} + rr \\left(- f_{5} + \\frac{f_{eq common}}{36} + \\rho \\left(\\frac{u0Mu1^{2}}{8} - \\frac{u0Mu1}{12} + \\frac{3 u_{0}^{2} u_{1}^{2}}{16} + \\frac{u_{0}^{2} u_{1}}{4} - \\frac{u_{0} u_{1}^{2}}{4}\\right)\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{6} \\leftarrow f_{6} + rr \\left(- f_{6} + \\frac{f_{eq common}}{36} + \\rho \\left(\\frac{u0Pu1^{2}}{8} + \\frac{u0Pu1}{12} + \\frac{3 u_{0}^{2} u_{1}^{2}}{16} + \\frac{u_{0}^{2} u_{1}}{4} + \\frac{u_{0} u_{1}^{2}}{4}\\right)\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{7} \\leftarrow f_{7} + rr \\left(- f_{7} + \\frac{f_{eq common}}{36} + \\rho \\left(\\frac{u0Pu1^{2}}{8} - \\frac{u0Pu1}{12} + \\frac{3 u_{0}^{2} u_{1}^{2}}{16} - \\frac{u_{0}^{2} u_{1}}{4} - \\frac{u_{0} u_{1}^{2}}{4}\\right)\\right)$$</td> </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$d_{8} \\leftarrow f_{8} + rr \\left(- f_{8} + \\frac{f_{eq common}}{36} + \\rho \\left(\\frac{u0Mu1^{2}}{8} + \\frac{u0Mu1}{12} + \\frac{3 u_{0}^{2} u_{1}^{2}}{16} - \\frac{u_{0}^{2} u_{1}}{4} + \\frac{u_{0} u_{1}^{2}}{4}\\right)\\right)$$</td> </tr> </table>"
+ ],
+ "text/plain": [
+ "AssignmentCollection: d_0, d_5, d_6, d_3, d_4, d_8, d_2, d_7, d_1 <- f(f_0, f_3, f_4, f_5, f_2, f_8, f_7, rr, f_1, f_6)"
+ ]
+ },
+ "execution_count": 18,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "simplification(collision_rule)"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index eef6d3a901941b2077c55c1906c97ba00110ab3e..5427a3ab174f2ecf7b9301c08c9026a577bacb07 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -179,6 +179,7 @@ For example, to modify the AST one can run::
"""
from copy import copy
+from warnings import warn
import sympy as sp
@@ -522,6 +523,7 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
'force_model': 'none',
'force': (0, 0, 0),
'maxwellian_moments': True,
+ 'central_moments': False,
'cumulant': False,
'initial_velocity': None,
@@ -613,6 +615,16 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
raise ValueError("Unknown optimization parameter(s): " + ",".join(unknown_opt_params))
params_result = copy(default_method_description)
+
+ # TODO: the equlibrium order should be set automatically for all methods
+ if 'central_moments' in params and 'equilibrium_order' in params:
+ if params['central_moments'] and params['equilibrium_order'] < 4:
+ warn("The central moments should be used with an equilibrium order depending on the highest order moment")
+
+ if 'central_moments' in params and 'equilibrium_order' not in params:
+ if params['central_moments']:
+ params['equilibrium_order'] = 4
+
params_result.update(params)
opt_params_result = copy(default_optimization_description)
opt_params_result.update(opt_params)
diff --git a/lbmpy/methods/momentbased.py b/lbmpy/methods/momentbased.py
index 53b901a25dfee07d82bfa6b75228607026714515..ae4074ae6d53fd620a41a69bf8589dd2f5c66a52 100644
--- a/lbmpy/methods/momentbased.py
+++ b/lbmpy/methods/momentbased.py
@@ -5,7 +5,7 @@ import sympy as sp
from lbmpy.maxwellian_equilibrium import get_weights
from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
-from lbmpy.moments import MOMENT_SYMBOLS, moment_matrix, shift_matrix
+from lbmpy.moments import MOMENT_SYMBOLS, moment_matrix, set_up_shift_matrix
from pystencils import Assignment
from pystencils.sympyextensions import subs_additive
@@ -17,6 +17,10 @@ class MomentBasedLbMethod(AbstractLbMethod):
These methods work by transforming the pdfs into moment space using a linear transformation. In the moment
space each component (moment) is relaxed to an equilibrium moment by a certain relaxation rate. These
equilibrium moments can e.g. be determined by taking the equilibrium moments of the continuous Maxwellian.
+ The collision can also be applied in the central moment space. For this Moment based LBM needs to set up a
+ shift matrix which is used to shift from moment space to central moment space. The shift matrix can be set up by
+ calling set_shift_matrix which either sets up the matrix based on the moments and the stencil of
+ Moment based LBM or takes an individual shift matrix as argument which is used for further calculations.
Args:
stencil: see :func:`lbmpy.stencils.get_stencil`
@@ -35,6 +39,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
self._momentToRelaxationInfoDict = OrderedDict(moment_to_relaxation_info_dict.items())
self._conservedQuantityComputation = conserved_quantity_computation
self._weights = None
+ self._shift_matrix = None
@property
def force_model(self):
@@ -74,13 +79,23 @@ class MomentBasedLbMethod(AbstractLbMethod):
self._weights = self._compute_weights()
return self._weights
+ def set_shift_matrix(self, shift_matrix=None):
+ if shift_matrix:
+ self._shift_matrix = shift_matrix
+ else:
+ self._shift_matrix = set_up_shift_matrix(self.moments, self.stencil)
+
def override_weights(self, weights):
assert len(weights) == len(self.stencil)
self._weights = weights
def get_equilibrium(self, conserved_quantity_equations=None, include_force_terms=False):
relaxation_matrix = sp.eye(len(self.relaxation_rates))
- return self._collision_rule_with_relaxation_matrix(relaxation_matrix,
+ if self.shift_matrix:
+ sm = sp.eye(len(self.relaxation_rates))
+ else:
+ sm = None
+ return self._collision_rule_with_relaxation_matrix(relaxation_matrix, sm,
conserved_quantity_equations=conserved_quantity_equations,
include_force_terms=include_force_terms)
@@ -90,8 +105,9 @@ class MomentBasedLbMethod(AbstractLbMethod):
def get_collision_rule(self, conserved_quantity_equations=None, keep_rrs_symbolic=True):
d = self.relaxation_matrix
+ sm = self.shift_matrix
relaxation_rate_sub_expressions, d = self._generate_relaxation_matrix(d, keep_rrs_symbolic)
- ac = self._collision_rule_with_relaxation_matrix(d, relaxation_rate_sub_expressions,
+ ac = self._collision_rule_with_relaxation_matrix(d, sm, relaxation_rate_sub_expressions,
True, conserved_quantity_equations)
return ac
@@ -120,13 +136,21 @@ class MomentBasedLbMethod(AbstractLbMethod):
def collision_matrix(self):
pdfs_to_moments = self.moment_matrix
d = self.relaxation_matrix
- return pdfs_to_moments.inv() * d * pdfs_to_moments
+ sm = self._shift_matrix
+ if sm:
+ return pdfs_to_moments.inv() * sm.inv() * d * sm * pdfs_to_moments
+ else:
+ return pdfs_to_moments.inv() * d * pdfs_to_moments
@property
def inverse_collision_matrix(self):
pdfs_to_moments = self.moment_matrix
inverse_relaxation_matrix = self.relaxation_matrix.inv()
- return pdfs_to_moments.inv() * inverse_relaxation_matrix * pdfs_to_moments
+ sm = self._shift_matrix
+ if sm:
+ return pdfs_to_moments.inv() * sm.inv() * inverse_relaxation_matrix * sm * pdfs_to_moments
+ else:
+ return pdfs_to_moments.inv() * inverse_relaxation_matrix * pdfs_to_moments
@property
def moment_matrix(self):
@@ -134,7 +158,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
@property
def shift_matrix(self):
- return shift_matrix(self.moments, self.stencil)
+ return self._shift_matrix
@property
def relaxation_matrix(self):
@@ -202,13 +226,18 @@ class MomentBasedLbMethod(AbstractLbMethod):
weights.append(value)
return weights
- def _collision_rule_with_relaxation_matrix(self, d, additional_subexpressions=(), include_force_terms=True,
+ def _collision_rule_with_relaxation_matrix(self, d, sm=None, additional_subexpressions=(), include_force_terms=True,
conserved_quantity_equations=None):
f = sp.Matrix(self.pre_collision_pdf_symbols)
pdf_to_moment = self.moment_matrix
m_eq = sp.Matrix(self.moment_equilibrium_values)
- collision_rule = f + pdf_to_moment.inv() * d * (m_eq - pdf_to_moment * f)
+ if sm:
+ print("Using collision with shift matrix")
+ collision_rule = f + pdf_to_moment.inv() * sm.inv() * d * sm * (m_eq - pdf_to_moment * f)
+ else:
+ print("Using collision without shift matrix")
+ collision_rule = f + pdf_to_moment.inv() * d * (m_eq - pdf_to_moment * f)
collision_eqs = [Assignment(lhs, rhs) for lhs, rhs in zip(self.post_collision_pdf_symbols, collision_rule)]
if conserved_quantity_equations is None:
diff --git a/lbmpy/moments.py b/lbmpy/moments.py
index 04a060c1ddbde0a026750e122c64162ff2138c7e..a78cc9a305abd2443004a5ca09aa689b59560901 100644
--- a/lbmpy/moments.py
+++ b/lbmpy/moments.py
@@ -369,9 +369,9 @@ def moment_matrix(moments, stencil, shift_velocity=None):
return sp.Matrix(len(moments), len(stencil), generator)
-def shift_matrix(moments, stencil):
+def set_up_shift_matrix(moments, stencil):
"""
- Returns shift matrix to shift raw moments to central moment space
+ Sets up a shift matrix to shift raw moments to central moment space
"""
x, y, z = MOMENT_SYMBOLS
dim = len(stencil[0])
diff --git a/lbmpy/simplificationfactory.py b/lbmpy/simplificationfactory.py
index af3a0d5854a19f82bc2d6967a5e082c7149bea96..43348e730c1cd8c373f013a44d1aae801b058003 100644
--- a/lbmpy/simplificationfactory.py
+++ b/lbmpy/simplificationfactory.py
@@ -16,22 +16,27 @@ def create_simplification_strategy(lb_method, split_inner_loop=False):
expand = apply_to_all_assignments(sp.expand)
if isinstance(lb_method, MomentBasedLbMethod):
- if len(set(lb_method.relaxation_rates)) <= 2:
- s.add(expand)
- s.add(replace_second_order_velocity_products)
- s.add(expand)
- s.add(factor_relaxation_rates)
- s.add(replace_density_and_velocity)
- s.add(replace_common_quadratic_and_constant_term)
- s.add(factor_density_after_factoring_relaxation_times)
- s.add(subexpression_substitution_in_main_assignments)
- if split_inner_loop:
- s.add(create_lbm_split_groups)
- s.add(add_subexpressions_for_divisions)
- else:
- s.add(subexpression_substitution_in_main_assignments)
- if split_inner_loop:
- s.add(create_lbm_split_groups)
+ if not lb_method.shift_matrix:
+ if len(set(lb_method.relaxation_rates)) <= 2:
+ s.add(expand)
+ s.add(replace_second_order_velocity_products)
+ s.add(expand)
+ s.add(factor_relaxation_rates)
+ s.add(replace_density_and_velocity)
+ s.add(replace_common_quadratic_and_constant_term)
+ s.add(factor_density_after_factoring_relaxation_times)
+ s.add(subexpression_substitution_in_main_assignments)
+ if split_inner_loop:
+ s.add(create_lbm_split_groups)
+ s.add(add_subexpressions_for_divisions)
+ else:
+ s.add(subexpression_substitution_in_main_assignments)
+ if split_inner_loop:
+ s.add(create_lbm_split_groups)
+ # else:
+ # s.add(expand)
+ # s.add(factor_relaxation_rates)
+ # s.add(add_subexpressions_for_divisions)
elif isinstance(lb_method, CumulantBasedLbMethod):
s.add(expand)
s.add(factor_relaxation_rates)
diff --git a/lbmpy_tests/test_moments.py b/lbmpy_tests/test_moments.py
index 565f436d2d104aaaaf08b97dd2a122e0b0224a9e..e23005ff9b8301b82df88607205d2955deba7df9 100644
--- a/lbmpy_tests/test_moments.py
+++ b/lbmpy_tests/test_moments.py
@@ -111,7 +111,7 @@ def test_shift_matrix():
method = create_lb_method(stencil=stencil, method='mrt_raw', compressible=True, equilibrium_order=6)
moments = method.moments
- sm = shift_matrix(moments, stencil)
+ sm = set_up_shift_matrix(moments, stencil)
m_eq = sp.Matrix(method.moment_equilibrium_values)
m_eq_rel = sp.simplify(sm * m_eq)
@@ -119,6 +119,6 @@ def test_shift_matrix():
m_eq_sol = sp.Matrix(method.moment_equilibrium_values)
for i in range(0, len(stencil)):
- m_eq_sol[i] = m_eq_sol[i].subs({sp.Symbol("u_0"):0, sp.Symbol("u_1"):0, sp.Symbol("u_2"):0})
+ m_eq_sol[i] = m_eq_sol[i].subs({sp.Symbol("u_0"): 0, sp.Symbol("u_1"): 0, sp.Symbol("u_2"): 0})
assert m_eq_rel == m_eq_sol