WIP( Work in progress): add monolithic CR–P0 Navier–Stokes solver #291
Conversation
|
Hello @jorgensd, I hope you are doing well. I wanted to kindly follow up on this PR. When you have a moment, I would greatly appreciate any initial feedback or workflow approval, so that I can continue developing and refining the contribution. Thank you very much for your time and consideration. |
Hi @akarve1507 |
|
Hi @jorgensd, I hope your week is going well. I wanted to gently check in regarding this PR. Any feedback whenever you have time would be very helpful. Thank you for your time. |
I think you should embed the markdown part within the tutorial. I have some other comments as well. |
jorgensd
left a comment
jorgensd
left a comment
There was a problem hiding this comment.
Good start, but a few things to consider along the way.
|
|
||
| --- | ||
|
|
||
| We solve the steady incompressible Navier–Stokes equations in a 2D channel using a **monolithic** formulation with **Crouzeix–Raviart (CR)** velocity and **DG-0** pressure. Dirichlet velocities on the left, top, and bottom are imposed **weakly via Nitsche’s method**, suitable for the nonconforming CR element. |
There was a problem hiding this comment.
| We solve the steady incompressible Navier–Stokes equations in a 2D channel using a **monolithic** formulation with **Crouzeix–Raviart (CR)** velocity and **DG-0** pressure. Dirichlet velocities on the left, top, and bottom are imposed **weakly via Nitsche’s method**, suitable for the nonconforming CR element. | |
| We solve the steady incompressible Navier–Stokes equations in a 2D channel using | |
| a **monolithic** formulation with [Crouzeix–Raviart](https://defelement.org/elements/crouzeix-raviart.html) | |
| (CR) for the velocity and **DG-0** pressure. | |
| We will impose Dirichlet conditions for the velocity field on the top, bottom and left part of the channel weakly using | |
| **Nitsche's method** (see [Weak imposition of Dirichlet conditions for the Poisson equation](../chapter1/nitsche]) for more details). |
|
|
||
| We solve the steady incompressible Navier–Stokes equations in a 2D channel using a **monolithic** formulation with **Crouzeix–Raviart (CR)** velocity and **DG-0** pressure. Dirichlet velocities on the left, top, and bottom are imposed **weakly via Nitsche’s method**, suitable for the nonconforming CR element. | ||
|
|
||
| --- |
| - **Linearization:** Picard (frozen convection); start from Stokes baseline | ||
| - **Stabilization:** small pressure mass $ \varepsilon_p \int_\Omega p\,q\,dx $ removes nullspace | ||
| - **Solver:** PETSc GMRES + Jacobi | ||
| - **Output:** XDMF for ParaView |
There was a problem hiding this comment.
| - **Output:** XDMF for ParaView | |
| - **Output:** {py:class}`VTXWriter<dolfinx.io.VTXWriter>` for Paraview |
As the velocity field is discontinuous this is a more suitable space.
| -\nu \Delta \mathbf{u} + (\mathbf{u}\cdot\nabla)\mathbf{u} + \nabla p = \mathbf{f}, | ||
| \quad \text{in } \Omega, | ||
| $$ | ||
|
|
||
| $$ | ||
| \nabla\cdot\mathbf{u} = 0, \quad \text{in } \Omega. |
There was a problem hiding this comment.
| -\nu \Delta \mathbf{u} + (\mathbf{u}\cdot\nabla)\mathbf{u} + \nabla p = \mathbf{f}, | |
| \quad \text{in } \Omega, | |
| $$ | |
| $$ | |
| \nabla\cdot\mathbf{u} = 0, \quad \text{in } \Omega. | |
| \begin{align} | |
| -\nu \Delta \mathbf{u} + (\mathbf{u}\cdot\nabla)\mathbf{u} + \nabla p &= \mathbf{f} && \text{in } \Omega,\\ | |
| \nabla\cdot\mathbf{u} &= 0 && \text{in } \Omega. | |
| \end{align} |
| a_{\text{core}}(\mathbf{u},p;\mathbf{v},q) | ||
| = \nu (\nabla\mathbf{u}, \nabla\mathbf{v})_\Omega | ||
| - (p, \nabla\cdot\mathbf{v})_\Omega | ||
| + (q, \nabla\cdot\mathbf{u})_\Omega, | ||
| $$ | ||
|
|
||
| $$ | ||
| a_{\text{conv}}(\mathbf{u};\mathbf{v}\mid\mathbf{u}_k) | ||
| = ((\mathbf{u}_k\cdot\nabla)\mathbf{u},\mathbf{v})_\Omega, | ||
| $$ | ||
|
|
||
| $$ | ||
| a_{\text{pm}}(p;q) | ||
| = \varepsilon_p (p,q)_\Omega. |
There was a problem hiding this comment.
| a_{\text{core}}(\mathbf{u},p;\mathbf{v},q) | |
| = \nu (\nabla\mathbf{u}, \nabla\mathbf{v})_\Omega | |
| - (p, \nabla\cdot\mathbf{v})_\Omega | |
| + (q, \nabla\cdot\mathbf{u})_\Omega, | |
| $$ | |
| $$ | |
| a_{\text{conv}}(\mathbf{u};\mathbf{v}\mid\mathbf{u}_k) | |
| = ((\mathbf{u}_k\cdot\nabla)\mathbf{u},\mathbf{v})_\Omega, | |
| $$ | |
| $$ | |
| a_{\text{pm}}(p;q) | |
| = \varepsilon_p (p,q)_\Omega. | |
| \begin{align} | |
| a_{\text{core}}(\mathbf{u},p;\mathbf{v},q) | |
| &= \nu (\nabla\mathbf{u}, \nabla\mathbf{v})_\Omega | |
| - (p, \nabla\cdot\mathbf{v})_\Omega | |
| + (q, \nabla\cdot\mathbf{u})_\Omega,\\ | |
| a_{\text{conv}}(\mathbf{u};\mathbf{v}\mid\mathbf{u}_k) | |
| &= ((\mathbf{u}_k\cdot\nabla)\mathbf{u},\mathbf{v})_\Omega\\ | |
| a_{\text{pm}}(p;q)&= \varepsilon_p (p,q)_\Omega. | |
| \end{align} |
There was a problem hiding this comment.
General comment here; Do you have a reference for this variational formulation?
| return ufl.lhs(F), ufl.rhs(F) | ||
|
|
||
|
|
||
| def solve_once(): |
There was a problem hiding this comment.
Is there a specific motivation behind not using dolfinx.fem.petsc.LinearProblem here?
It just means that one needs to do a lot of PETSc object management.
| for tag in (3, 4): | ||
| F += ( | ||
| - nu * ufl.inner(ufl.grad(u) * n, v) | ||
| - nu * ufl.inner(ufl.grad(v) * n, (u - zero_vec)) | ||
| + alpha * nu / h * ufl.inner(u - zero_vec, v) | ||
| - ufl.inner(p * n, v) | ||
| + q * ufl.inner(u - zero_vec, n) | ||
| ) * ds(tag) |
There was a problem hiding this comment.
| for tag in (3, 4): | |
| F += ( | |
| - nu * ufl.inner(ufl.grad(u) * n, v) | |
| - nu * ufl.inner(ufl.grad(v) * n, (u - zero_vec)) | |
| + alpha * nu / h * ufl.inner(u - zero_vec, v) | |
| - ufl.inner(p * n, v) | |
| + q * ufl.inner(u - zero_vec, n) | |
| ) * ds(tag) | |
| F += ( | |
| - nu * ufl.inner(ufl.grad(u) * n, v) | |
| - nu * ufl.inner(ufl.grad(v) * n, (u - zero_vec)) | |
| + alpha * nu / h * ufl.inner(u - zero_vec, v) | |
| - ufl.inner(p * n, v) | |
| + q * ufl.inner(u - zero_vec, n) | |
| ) * ds((3,4)) |
| F += ( | ||
| - nu * ufl.inner(ufl.grad(u) * n, v) # consistency | ||
| - nu * ufl.inner(ufl.grad(v) * n, (u - u_in)) # symmetry | ||
| + alpha * nu / h * ufl.inner(u - u_in, v) # penalty | ||
| - ufl.inner(p * n, v) # momentum BC consistency | ||
| + q * ufl.inner(u - u_in, n) # continuity BC consistency: q * ((u-u_in)·n) | ||
| ) * ds(1) |
There was a problem hiding this comment.
Since we are using the same formulation here and below, I think we should make a convenience function
def nitsche_terms(g, tags):
return (-nu * ufl.inner(ufl.grad(u)*n, v) - nu * ufl.inner(ufl.inner(grad(v) * n, (u - g)) + alpha * nu / h + ufl.inner(u - g, v) + .... ) * ds(tags)so that we can do
F += nitsche_terms(zero_vec, (3, 4))
F += nitsche_terms(u_in, (1, ))| def build_forms(): | ||
| n = ufl.FacetNormal(domain) | ||
| h = ufl.CellDiameter(domain) | ||
| alpha = PETSc.ScalarType(args.alpha) |
There was a problem hiding this comment.
Wrap alpha as a dolfinx.fem.Constant to avoid recompilation.
| zero = fem.Constant(domain, PETSc.ScalarType(0.0)) | ||
| f_vec = ufl.as_vector((zero, zero)) # zero body force |
There was a problem hiding this comment.
| zero = fem.Constant(domain, PETSc.ScalarType(0.0)) | |
| f_vec = ufl.as_vector((zero, zero)) # zero body force | |
| f_vec = fem.Constant(domain, np.zeros(mesh.geometry.dim, dtype=dolfinx.default_scalar_type)) |
2 participants
Summary
This pull request introduces an initial monolithic CR–P0 steady Navier–Stokes solver.
Current contents
chapter2/monolithic_navierstokes_crp0.py— runnable solver with:--stokes-onlymode for testing)Work in progress
.md) page will be added in a future commitmkdocs.ymlupdate pendingStatus
Code compiles and runs successfully
Tutorial and visualization assets pending
Opened as a work-in-progress (WIP) PR in line with project guidelines
I kindly request feedback from the maintainers to ensure that the implementation aligns with the project’s coding standards and tutorial structure. I am happy to revise and improve this contribution based on your suggestions before completing the documentation and visualization components.