diff --git a/applications/FluidDynamicsApplication/python_scripts/symbolic_generation/low_mach_navier_stokes/low_mach_navier_stokes_element_derivation.ipynb b/applications/FluidDynamicsApplication/python_scripts/symbolic_generation/low_mach_navier_stokes/low_mach_navier_stokes_element_derivation.ipynb
index 459d80347a5..b749d81eaa3 100644
--- a/applications/FluidDynamicsApplication/python_scripts/symbolic_generation/low_mach_navier_stokes/low_mach_navier_stokes_element_derivation.ipynb
+++ b/applications/FluidDynamicsApplication/python_scripts/symbolic_generation/low_mach_navier_stokes/low_mach_navier_stokes_element_derivation.ipynb
@@ -17,7 +17,7 @@
     "\\begin{equation*}\n",
     "\\begin{aligned}\n",
     "    \\frac{\\partial \\rho}{\\partial t} + \\nabla\\cdot\\left(\\rho\\bm{u}\\right) &= 0 \\\\\n",
-    "    \\rho\\frac{\\partial\\bm{u}}{\\partial t} + \\rho\\bm{u}\\cdot\\nabla\\bm{u} - \\nabla\\cdot\\boldsymbol{\\tau} + \\nabla p + \\sigma(\\bm{u} - \\bm{u}_{s}) + \\sigma\\left(\\bm{u} - \\bm{u}_{s}\\right) &= \\rho\\bm{g} \\\\\n",
+    "    \\rho\\frac{\\partial\\bm{u}}{\\partial t} + \\rho\\bm{u}\\cdot\\nabla\\bm{u} - \\nabla\\cdot\\boldsymbol{\\tau} + \\nabla p + \\sigma(\\bm{u} - \\bm{u}_{s}) &= \\rho\\bm{g} \\\\\n",
     "    \\rho c_{p}\\frac{\\partial T}{\\partial t} + \\rho c_{p}\\bm{u}\\cdot\\nabla T - \\nabla\\cdot\\left(\\kappa\\nabla T\\right) - \\alpha T\\frac{d p_{th}}{dt} &= Q\n",
     "\\end{aligned}\n",
     "\\; ,\n",
@@ -50,65 +50,78 @@
     "\\end{equation*}\n",
     "where the subscript $0$ indicates the initial magnitudes at $t=0$ and $\\Omega$ the problem domain.\n",
     "\n",
-    "<!-- ## Variational form\n",
-    "The variational form (i.e. the functional to be automatically differenctiated) is obtained by multiplying previous equations by the corresponding test functions $\\mathbf{v}$, $q$ and $w$. After integration by parts, this yields\n",
-    "\n",
-    "\\begin{align*}\n",
-    "\\left(w_{x},\\rho g_{x}\\right)_{\\Omega} - & \\left(w_{x},\\rho\\frac{\\partial u_{x}}{\\partial t}\\right)_{\\Omega} - \\left(w_{x},\\rho a_{y}\\frac{\\partial u_{x}}{\\partial y}\\right)_{\\Omega} - \\left(w_{x},\\rho a_{x}\\frac{\\partial u_{x}}{\\partial x}\\right)_{\\Omega} + \\left(\\frac{\\partial w_{x}}{\\partial x}, p\\right)_{\\Omega} \\\\\n",
-    "- & \\left(\\frac{\\partial w_{x}}{\\partial y}, \\mu\\frac{\\partial u_{x}}{\\partial y}\\right)_{\\Omega} - \\left(\\frac{\\partial w_{x}}{\\partial x}, \\mu\\frac{\\partial u_{x}}{\\partial x}\\right)_{\\Omega} + \\langle w_{x}, \\mu\\left(\\frac{\\partial u_{x}}{\\partial y}+\\frac{\\partial u_{x}}{\\partial x}\\right)n_{x} - p n_{x} \\rangle_{\\Gamma} \\\\\n",
-    "+ & \\left(w_{y},\\rho g_{y}\\right)_{\\Omega} - \\left(w_{y},\\rho\\frac{\\partial u_{y}}{\\partial t}\\right)_{\\Omega} - \\left(w_{y},\\rho a_{y}\\frac{\\partial u_{y}}{\\partial y}\\right)_{\\Omega} - \\left(w_{y},\\rho a_{x}\\frac{\\partial u_{y}}{\\partial x}\\right)_{\\Omega}\n",
-    "+ \\left(\\frac{w_{y}}{y} + \\frac{\\partial w_{y}}{\\partial y}, p\\right)_{\\Omega} \\\\\n",
-    "- & \\left(\\frac{\\partial w_{y}}{\\partial y}, \\mu\\frac{\\partial u_{y}}{\\partial y}\\right)_{\\Omega} - \\left(\\frac{\\partial w_{y}}{\\partial x}, \\mu\\frac{\\partial u_{y}}{\\partial x}\\right)_{\\Omega} + \\langle w_{y}, \\mu\\left(\\frac{\\partial u_{y}}{\\partial y}+\\frac{\\partial u_{y}}{\\partial x}\\right)n_{y} - p n_{y} \\rangle_{\\Gamma} \\\\\n",
-    "- & \\left(q,\\frac{u_{y}}{y}\\right)_{\\Omega} - \\left(q,\\frac{\\partial u_{y}}{\\partial y}\\right)_{\\Omega} - \\left(q, \\frac{\\partial u_{x}}{\\partial x}\\right)_{\\Omega} = 0\\; ,\n",
-    "\\end{align*}\n",
-    "being $\\mathbf{a}$ the convective velocity, which computed from the linearised previous iteration velocity $\\mathbf{u}^{k-1}$ and the mesh velocity $\\mathbf{u}_{M}$ as\n",
+    "## Variational form\n",
+    "The variational form (i.e. the functional to be automatically differenctiated) is obtained by multiplying previous equations by the corresponding pressure, velocity and temperature test functions $q$, $\\mathbf{v}$ and $w$. After integration by parts, this yields\n",
+    "\n",
     "\\begin{equation*}\n",
-    "\\mathbf{a} = \\mathbf{u}^{k-1} - \\mathbf{u}_{M}\\; .\n",
+    "\\begin{aligned}\n",
+    "    - \\left(q, \\frac{\\partial \\rho}{\\partial t}\\right)_{\\Omega} &- \\left(q, \\nabla\\cdot(\\rho\\mathbf{u})\\right)_{\\Omega} - \\left(\\mathbf{v}, \\rho\\frac{\\partial\\mathbf{u}}{\\partial t}\\right)_{\\Omega} - \\left(\\mathbf{v}, \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u}\\right)_{\\Omega} - \\left(\\nabla^{s}\\mathbf{v}, \\boldsymbol{\\tau}\\right)_{\\Omega} + \\left(\\nabla\\cdot\\mathbf{v}, p\\right)_{\\Omega} \\\\\n",
+    "    &- \\left(\\mathbf{v}, \\sigma(\\mathbf{u} - {\\mathbf{u}_{s}})\\right)_{\\Omega} + \\langle\\mathbf{v}, \\mathbf{t}\\rangle_{\\Gamma} + \\left(\\mathbf{v}, \\rho\\mathbf{g}\\right)_{\\Omega} + \\left(w, Q\\right)_{\\Omega} - \\left(w, \\rho c_{p} \\frac{\\partial T}{\\partial t}\\right)_{\\Omega} \\\\\n",
+    "    &- \\left(w, \\rho c_p \\mathbf{a}\\cdot\\nabla T\\right)_{\\Omega} - \\left(\\nabla w, \\kappa\\nabla T\\right)_{\\Omega} + \\langle w, \\kappa\\nabla T\\cdot \\mathbf{n}\\rangle_{\\Gamma} + \\left(w, \\alpha T \\frac{d p_{th}}{dt}\\right)_{\\Omega} = 0\\; ,\n",
+    "\\end{aligned}\n",
+    "\\end{equation*}\n",
+    "being $\\mathbf{a}$ the convective velocity, which computed from the linearised previous iteration velocity $\\mathbf{u}^{k-1}$ and the mesh velocity $\\mathbf{u}_{m}$ as\n",
+    "\\begin{equation*}\n",
+    "\\mathbf{a} = \\mathbf{u}^{k-1} - \\mathbf{u}_{m}\\; .\n",
     "\\end{equation*}\n",
-    "\n",
-    "Note that the cylindrical component of the divergence operator applied to the viscous stress term (see the strong form above) no longer appears in the variational form as we get rid of it by integrating by parts.\n",
-    "However, it appears in the test function divergence of the pressure term.\n",
     "\n",
     "## Variational MultiScales stabilization\n",
     "The scales separation reads\n",
     "\n",
-    "\\begin{align*}\n",
-    "\\mathbf{u} & = \\mathbf{u_{h}}+\\mathbf{u_{s}}\\\\\n",
-    "p & = p_{h}+p_{s}\n",
-    "\\end{align*}\n",
+    "\\begin{equation*}\n",
+    "\\begin{aligned}\n",
+    "    p &= p_{h} + \\tilde{p} \\\\\n",
+    "    \\mathbf{u} &= \\mathbf{u}_{h} + \\tilde{\\mathbf{u}} \\\\\n",
+    "    T &= T_{h} + \\tilde{T}\n",
+    "\\end{aligned}\\; ,\n",
+    "\\end{equation*}\n",
+    "\n",
+    "being $p_{h}$, $\\mathbf{u}_{h}$ and $T_{h}$ the finite element scales and $\\tilde{p}$, $\\tilde{\\mathbf{u}}$ and $\\tilde{T}$ the subgrid scales, which are modelled from the finite element residuals as\n",
+    "\n",
+    "\\begin{equation*}\n",
+    "\\begin{aligned}\n",
+    "\\tilde{p} & := \\tau_{c} R_{c}(\\mathbf{u}_{h},T_{h}) \\\\\n",
+    "\\tilde{\\mathbf{u}} & := \\tau_{m} \\mathbf{R}_{m}(p_{h}, \\mathbf{u}_{h}, T_{h}) \\\\\n",
+    "\\tilde{T} & := \\tau_{e} R_{e}(\\mathbf{u}_{h}, T_{h}) \\\\\n",
+    "\\end{aligned}\\; ,\n",
+    "\\end{equation*}\n",
+    "\n",
+    "being $\\tau_{c}$, $\\tau_{m}$ and $\\tau_{e}$ the stabilization parameters which are equal to\n",
     "\n",
-    "being $\\mathbf{u}_{h}$ and $p_{h}$ the finite element scales and $\\mathbf{u}_{s}$ and $p_{s}$ the subgrid scales, which are modelled from the finite element residuals as\n",
+    "\\begin{equation*}\n",
+    "\\begin{aligned}\n",
+    "    \\tau_{c} &:= \\frac{\\mu}{\\rho} + \\frac{c_{2}}{c_{1}}h\\lVert u\\rVert +\\frac{c_{3}}{c_{1}}\\frac{\\sigma h}{\\rho} \\\\\n",
+    "    \\tau_{u} &:= \\left(c_{1}\\frac{\\mu}{h^2} + c_{2}\\frac{\\rho\\lVert u\\rVert}{h} + c_{3} \\frac{\\sigma}{h}\\right)^{-1} \\\\\n",
+    "    \\tau_{e} &:= \\left(c_{1}\\frac{\\kappa}{h^2} + c_{2}\\frac{\\rho c_{p}\\lVert u\\rVert}{h}\\right)^{-1}\n",
+    "\\end{aligned} \\; .\n",
+    "\\end{equation*}\n",
     "\n",
-    "\\begin{align*}\n",
-    "u_{s,x} & = \\tau_{1} R^{M}_{x}(u_{h,x},p_{h})\\; , \\\\\n",
-    "u_{s,y} & = \\tau_{1} R^{M}_{y}(u_{h,y},p_{h})\\; , \\\\\n",
-    "p_{s} & = \\tau_{2} R^{C}(\\mathbf{u}_{h},p_{h})\\; ,\n",
-    "\\end{align*}\n",
+    "$c_{1}$, $c_{2}$, and $c_{3}$ are algorithmic constants defaulted in the header file to $4$, $2$, and $2$, respectively. $h$ the characteristic element size.\n",
+    "\n",
+    "The finite element residuals $R^{c}(\\mathbf{u}_{h}, T_{h})$, $\\mathbf{R}^{m}(p_{h}, \\mathbf{u}_{h}, T_{h})$ and $R^{e}(p_{h}, \\mathbf{u}_{h}, T_{h})$ are defined as\n",
     "\n",
-    "being $\\tau_{1}$ and $\\tau_{2}$ the stabilization parameters which in this case equal\n",
     "\\begin{equation*}\n",
-    "    \\tau_{1} = \\left(\\frac{\\rho\\tau_{d}}{\\Delta t} + \\frac{c_{1}\\mu}{h^{2}} + \\frac{c_{2}\\rho\\lVert\\mathbf{a}\\rVert}{h}\\right)^{-1}\n",
+    "\\begin{aligned}\n",
+    "    &R_{c}(\\mathbf{u}_{h}, T_{h}) = -\\frac{\\partial\\rho}{\\partial t} + \\alpha\\nabla T_{h}\\cdot\\mathbf{u}_{h} - \\rho\\nabla\\cdot\\mathbf{u}_{h}\\\\\n",
+    "    &\\mathbf{R}_{m}(p_{h}, \\mathbf{u}_{h}, T_{h}) = \\rho\\mathbf{g} - \\rho\\frac{\\partial\\mathbf{u}_{h}}{\\partial t} - \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u}_{h} - \\nabla p - \\sigma(\\mathbf{u}_{h} - \\mathbf{u}_{s})\\\\\n",
+    "    &R_{e}(p_{h}, \\mathbf{u}_{h}, T_{h}) = Q - \\rho c_p \\frac{\\partial T}{\\partial t} - \\rho c_p \\mathbf{u}_{h}\\nabla T_{h} + \\alpha T_{h} \\frac{d p_{th}}{dt}\n",
+    "\\end{aligned}\n",
     "\\end{equation*}\n",
-    "and\n",
+    "\n",
+    "Once the subscales are defined, one can derive the stabilization terms to be added to the Galerkin functional. Hence, after doing the required integration by parts to remove the subscale derivatives and dropping the high order terms (the use of linear elements is assumed), one obtains the final functional to which the automatic differentiation is to be applied\n",
+    "\n",
     "\\begin{equation*}\n",
-    "    \\tau_{2} = \\mu+\\frac{c_{2}\\rho h\\lVert\\mathbf{a}\\rVert}{c_{1}} \\; ,\n",
+    "\\begin{aligned}\n",
+    "    - \\left(q, \\frac{\\partial \\rho}{\\partial t}\\right)_{\\Omega} &- \\left(q, \\nabla\\rho\\cdot\\mathbf{u}_{h}\\right)_{\\Omega} - \\left(q, \\rho\\nabla\\cdot\\mathbf{u}_{h}\\right)_{\\Omega} + \\left(\\nabla q, \\rho\\tilde{\\mathbf{u}}\\right)_{\\Omega} - \\left(\\mathbf{v}, \\rho\\frac{\\partial\\mathbf{u}_{h}}{\\partial t}\\right)_{\\Omega} - \\left(\\mathbf{v}, \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u}_{h}\\right)_{\\Omega} \\\\\n",
+    "    &+ \\left((\\nabla\\rho\\cdot\\mathbf{a})\\mathbf{v}, \\tilde{\\mathbf{u}}\\right)_{\\Omega} + \\left((\\rho\\nabla\\cdot\\mathbf{a})\\mathbf{v}, \\tilde{\\mathbf{u}} \\right)_{\\Omega} + \\left(\\rho\\mathbf{a}\\cdot\\nabla\\mathbf{v}, \\tilde{\\mathbf{u}}\\right)_{\\Omega} - \\left(\\nabla^{s}\\mathbf{v}, \\boldsymbol{\\tau}\\right)_{\\Omega} \\\\\n",
+    "    &+ \\left(\\nabla\\cdot\\mathbf{v}, p_{h}\\right)_{\\Omega} + \\left(\\nabla\\cdot\\mathbf{v}, \\tilde{p}\\right)_{\\Omega} - \\left(\\mathbf{v}, \\sigma(\\mathbf{u}_{h} + \\tilde{\\mathbf{u}} - {\\mathbf{u}_{s}})\\right)_{\\Omega} + \\left(\\mathbf{v}, \\rho\\mathbf{g}\\right)_{\\Omega} \\\\\n",
+    "    &+ \\left(w, Q\\right)_{\\Omega} - \\left(w, \\rho c_{p} \\frac{\\partial T_{h}}{\\partial t}\\right)_{\\Omega} - \\left(w, \\rho c_p \\mathbf{a}\\cdot\\nabla T_{h}\\right)_{\\Omega} + \\left((\\nabla w)\\rho c_p \\mathbf{a}, \\tilde{T}\\right)_{\\Omega} \\\\\n",
+    "    &+ \\left(wc_p (\\nabla\\rho\\cdot\\mathbf{a}), \\tilde{T} \\right)_{\\Omega} + \\left(w\\rho c_p (\\nabla\\cdot\\mathbf{a}), \\tilde{T}\\right)_{\\Omega} - \\left(\\nabla w, \\kappa\\nabla T_h\\right)_{\\Omega} + \\left(w, \\alpha (T_h + \\tilde{T}) \\frac{d p_{th}}{dt}\\right)_{\\Omega} = 0\n",
+    "\\end{aligned}\n",
     "\\end{equation*}\n",
-    "being $c_{1}$ and $c_{2}$ the user-definable stabilization constants, $\\tau_{d}$ a dynamic coefficient ranging between 0 and 1 and $h$ the characteristic element size.\n",
-    "\n",
-    "The finite element residuals $\\mathbf{R}^{M}(\\mathbf{u}_{h},p_{h}$ and $R^{C}(\\mathbf{u}_{h},p_{h})$ are defined as\n",
-    "\\begin{align*}\n",
-    "R^{M}_{x}(u_{h,x},p_{h}) = & \\rho g_{x} - \\rho\\frac{\\partial u_{h,x}}{\\partial t} - \\rho a_{h,y}\\frac{\\partial u_{h,x}}{\\partial y} - \\rho a_{h,x}\\frac{\\partial u_{h,x}}{\\partial x} - \\frac{\\partial p_{h}}{\\partial x} + \\mu\\left(\\frac{1}{y}\\frac{\\partial u_{h,x}}{\\partial y} + \\frac{\\partial^{2} u_{h,x}}{\\partial y^{2}} + \\frac{\\partial^{2} u_{h,x}}{\\partial x^{2}}\\right)\\; , \\\\\n",
-    "R^{M}_{y}(u_{h,y},p_{h}) = & \\rho g_{y} - \\rho\\frac{\\partial u_{h,y}}{\\partial t} - \\rho a_{h,y}\\frac{\\partial u_{h,y}}{\\partial y} - \\rho a_{h,x}\\frac{\\partial u_{h,y}}{\\partial x} - \\frac{\\partial p_{h}}{\\partial y} + \\mu\\left(- \\frac{u_{h,y}}{y^{2}} + \\frac{1}{y}\\frac{\\partial u_{h,y}}{\\partial y} + \\frac{\\partial^{2} u_{h,y}}{\\partial y^{2}} + \\frac{\\partial^{2} u_{h,y}}{\\partial x^{2}}\\right)\\; , \\\\\n",
-    "R^{C}(\\mathbf{u}_{h},p_{h}) = & -\\frac{1}{y}u_{h,y} - \\frac{\\partial u_{h,y}}{\\partial y} - \\frac{\\partial u_{h,x}}{\\partial x} \\; .\n",
-    "\\end{align*}\n",
-    "\n",
-    "Once the subscales are defined, one can derive the stabilization terms to be added to the functional above. Hence, after doing the required integration by parts to remove the subscale derivatives and dropping the high order terms (the use of linear elements is assumed), one obtains the stabilization functional\n",
-    "\\begin{align*}\n",
-    "\\left(\\frac{\\partial w_{x}}{\\partial y}\\rho a_{y} + w_{x}\\rho\\frac{\\partial a_{y}}{\\partial y}, u_{s,x}\\right)_{\\Omega} & + \\left(\\frac{\\partial w_{x}}{\\partial x}\\rho a_{x} + w_{x}\\rho\\frac{\\partial a_{x}}{\\partial x}, u_{s,x}\\right)_{\\Omega} + \\left(\\frac{\\partial w_{x}}{\\partial x},p_{s}\\right)_{\\Omega} \\\\\n",
-    "& + \\left(\\frac{\\partial w_{y}}{\\partial y}\\rho a_{y} + w_{y}\\rho\\frac{\\partial a_{y}}{\\partial y}, u_{s,y}\\right)_{\\Omega} + \\left(\\frac{\\partial w_{y}}{\\partial x}\\rho a_{x} + w_{y}\\rho\\frac{\\partial a_{x}}{\\partial x}, u_{s,y}\\right)_{\\Omega} + \\left(\\frac{w_{y}}{y} + \\frac{\\partial w_{y}}{\\partial y},p_{s}\\right)_{\\Omega} \\\\\n",
-    "& - \\left(q,\\frac{1}{y}u_{s,y}\\right)_{\\Omega} + \\left(\\frac{\\partial q}{\\partial y},u_{s,y}\\right)_{\\Omega} + \\left(\\frac{\\partial q}{\\partial x},u_{s,x}\\right)_{\\Omega}\n",
-    "\\end{align*}\n",
-    "in which the quasi-static subscales assumption has been considered. -->\n",
+    "\n",
+    "Note that the quasi-static subscales assumption is taken so all the subscales time derivatives are dropped.\n",
+    "\n",
     "<!-- \n",
     "## IMPLEMENTATION\n",
     "The remaining implementation is similar to that of the weakly compressible Navier-Stokes element but with the added simplification of assuming a Newtonian constitutive model within the element.  --> -->\n"