Corrections on the N-S code

chapter2/navierstokes.md

@@ -75,24 +75,24 @@ We now move on to the second step in  our splitting scheme for the incompressibl
 We may now use the computed tentative velocity to compute the new pressure $p^n$:
 ```{math}
 :label: ipcs-two
-    \langle \nabla p^{n+1}, \nabla q \rangle = \langle \nabla p^n, \nabla q\rangle - \Delta t^{-1}\langle \nabla \cdot u^*, q\rangle.
+    \langle \nabla p^{n+1}, \nabla q \rangle = \langle \nabla p^n, \nabla q\rangle - \frac{\rho}{\Delta t}\langle \nabla \cdot u^*, q\rangle.
 ```
 Note here that $q$ is a scalar-valued test function from the pressure space, whereas the test function $v$ in [](ipcs-one) is a vector-valued test function from the velocity space.
 
 One way to think about this step is to subtract the Navier-Stokes momentum equation [](navier-stokes) expressed in terms of the tentative velocity $u^*$ and the pressure $p^n$ from the momentum equation expressed in terms of the velocity $u^{n+1}$ and pressure $p^{n+1}$. This results in the equation
 ```{math}
-\frac{u^{n+1}-u^*}{\Delta t}+\nabla p^{n+1}- \nabla p^n = 0.
+\frac{\rho (u^{n+1}-u^*)}{\Delta t}+\nabla p^{n+1}- \nabla p^n = 0.
 ```
 Taking the divergence and requiring that $\nabla \cdot u^{n+1}=0$ by the Navier-Stokes continuity equation, we obtain the equation
 ```{math}
 :label: ipcs-tmp
- - \frac{\nabla\cdot  u^*}{\Delta t}+ \nabla^2p^{n+1}-\nabla^2p^n=0,
+ - \frac{\rho \nabla\cdot  u^*}{\Delta t}+ \nabla^2p^{n+1}-\nabla^2p^n=0,
 ```
 which is the Poisson problem for the pressure $p^{n+1} resulting in the variational formulation [](ipcs-two).
 
 Finally, we compute the corrected velocity $u^{n+1}$ from the equation [](ipcs-tmp). Multiplying this equation by a test function $v$, we obtain
 ```{math}
-    \langle u^{n+1}, v\rangle=\langle u^*, v\rangle -\Delta t\langle \nabla(p^{n+1}-p^n), v\rangle
+    \rho \langle (u^{n+1} - u^*), v\rangle= -\Delta t\langle \nabla(p^{n+1}-p^n), v\rangle
 ```
 
 In summary, we may thus solve the incompressible Navier-Stokes equations efficiently by solving a sequence of three linear variational problems in each step.

chapter2/ns_code1.ipynb

@@ -282,15 +282,15 @@
     "# Define variational problem for step 2\n",
     "u_ = Function(V)\n",
     "a2 = form(dot(nabla_grad(p), nabla_grad(q))*dx)\n",
-    "L2 = form(dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx)\n",
+    "L2 = form(dot(nabla_grad(p_n), nabla_grad(q))*dx - (rho/k)*div(u_)*q*dx)\n",
     "A2 = assemble_matrix(a2, bcs=bcp)\n",
     "A2.assemble()\n",
     "b2 = create_vector(L2)\n",
     "\n",
     "# Define variational problem for step 3\n",
     "p_ = Function(Q)\n",
-    "a3 = form(dot(u, v)*dx)\n",
-    "L3 = form(dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx)\n",
+    "a3 = form(rho*dot(u, v)*dx)\n",
+    "L3 = form(rho*dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx)\n",
     "A3 = assemble_matrix(a3)\n",
     "A3.assemble()\n",
     "b3 = create_vector(L3)"

chapter2/ns_code1.py

@@ -180,15 +180,15 @@ def sigma(u, p):
 # Define variational problem for step 2
 u_ = Function(V)
 a2 = form(dot(nabla_grad(p), nabla_grad(q))*dx)
-L2 = form(dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx)
+L2 = form(dot(nabla_grad(p_n), nabla_grad(q))*dx - (rho/k)*div(u_)*q*dx)
 A2 = assemble_matrix(a2, bcs=bcp)
 A2.assemble()
 b2 = create_vector(L2)
 
 # Define variational problem for step 3
 p_ = Function(Q)
-a3 = form(dot(u, v)*dx)
-L3 = form(dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx)
+a3 = form(rho*dot(u, v)*dx)
+L3 = form(rho*dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx)
 A3 = assemble_matrix(a3)
 A3.assemble()
 b3 = create_vector(L3)

chapter2/ns_code2.ipynb

@@ -372,7 +372,7 @@
     "$$\n",
     "where we have used the two previous time steps in the temporal derivative for the velocity, and compute the pressure staggered in time, at the time between the previous and current solution. The second step becomes\n",
     "$$\n",
-    "-\\nabla \\phi = -\\frac{1}{\\delta t} \\nabla \\cdot u^* \\qquad\\text{in } \\Omega,\n",
+    "\\nabla \\phi = -\\frac{\\rho}{\\delta t} \\nabla \\cdot u^* \\qquad\\text{in } \\Omega,\n",
     "$$\n",
     "\n",
     "$$\n",
@@ -387,7 +387,7 @@
     "Finally, the third step is\n",
     "\n",
     "$$\n",
-    "u^{n+1} = u^{*}-\\delta t \\phi. \n",
+    "\\rho (u^{n+1}-u^{*}) = -\\delta t \\phi. \n",
     "$$\n",
     "\n",
     "We start by defining all the variables used in the variational formulations."
@@ -451,7 +451,7 @@
    "outputs": [],
    "source": [
     "a2 = form(dot(grad(p), grad(q))*dx)\n",
-    "L2 = form(-1/k * dot(div(u_s), q) * dx)\n",
+    "L2 = form(-rho / k * dot(div(u_s), q) * dx)\n",
     "A2 = assemble_matrix(a2, bcs=bcp)\n",
     "A2.assemble()\n",
     "b2 = create_vector(L2)"
@@ -470,8 +470,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "a3 = form(dot(u, v)*dx)\n",
-    "L3 = form(dot(u_s, v)*dx - k * dot(nabla_grad(phi), v)*dx)\n",
+    "a3 = form(rho * dot(u, v)*dx)\n",
+    "L3 = form(rho * dot(u_s, v)*dx - k * dot(nabla_grad(phi), v)*dx)\n",
     "A3 = assemble_matrix(a3)\n",
     "A3.assemble()\n",
     "b3 = create_vector(L3)"
@@ -805,7 +805,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.4"
+   "version": "3.10.6"
   }
  },
  "nbformat": 4,

chapter2/ns_code2.py

@@ -6,7 +6,7 @@
 #       extension: .py
 #       format_name: light
 #       format_version: '1.5'
-#       jupytext_version: 1.14.1
+#       jupytext_version: 1.13.8
 #   kernelspec:
 #     display_name: Python 3 (ipykernel)
 #     language: python
@@ -228,7 +228,7 @@ def __call__(self, x):
 # $$
 # where we have used the two previous time steps in the temporal derivative for the velocity, and compute the pressure staggered in time, at the time between the previous and current solution. The second step becomes
 # $$
-# -\nabla \phi = -\frac{1}{\delta t} \nabla \cdot u^* \qquad\text{in } \Omega,
+# \nabla \phi = -\frac{\rho}{\delta t} \nabla \cdot u^* \qquad\text{in } \Omega,
 # $$
 #
 # $$
@@ -243,7 +243,7 @@ def __call__(self, x):
 # Finally, the third step is
 #
 # $$
-# u^{n+1} = u^{*}-\delta t \phi. 
+# \rho (u^{n+1}-u^{*}) = -\delta t \phi. 
 # $$
 #
 # We start by defining all the variables used in the variational formulations.
@@ -265,9 +265,9 @@ def __call__(self, x):
 
 f = Constant(mesh, PETSc.ScalarType((0,0)))
 F1 = rho / k * dot(u - u_n, v) * dx 
-F1 += inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
+F1 += rho * inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
 F1 += 0.5 * mu * inner(grad(u + u_n), grad(v))*dx - dot(p_, div(v))*dx
-F1 += dot(f, v) * dx
+F1 -= dot(f, v) * dx
 a1 = form(lhs(F1))
 L1 = form(rhs(F1))
 A1 = create_matrix(a1)
@@ -276,15 +276,15 @@ def __call__(self, x):
 # Next we define the second step
 
 a2 = form(dot(grad(p), grad(q))*dx)
-L2 = form(-1/k * dot(div(u_s), q) * dx)
+L2 = form(-rho / k * dot(div(u_s), q) * dx)
 A2 = assemble_matrix(a2, bcs=bcp)
 A2.assemble()
 b2 = create_vector(L2)
 
 # We finally create the last step
 
-a3 = form(dot(u, v)*dx)
-L3 = form(dot(u_s, v)*dx - k * dot(nabla_grad(phi), v)*dx)
+a3 = form(rho * dot(u, v)*dx)
+L3 = form(rho * dot(u_s, v)*dx - k * dot(nabla_grad(phi), v)*dx)
 A3 = assemble_matrix(a3)
 A3.assemble()
 b3 = create_vector(L3)