WIP( Work in progress): add monolithic CR–P0 Navier–Stokes solver

chapter2/monolithic_navierstokes_crp0.md

@@ -0,0 +1,97 @@
+# Monolithic Incompressible Navier–Stokes with CR–P0 (Steady Picard)
+
+We solve the steady incompressible Navier–Stokes equations in a 2D channel.
+
+We use a **monolithic** formulation with  
+[Crouzeix–Raviart (CR)](https://defelement.org/elements/crouzeix-raviart.html)  
+for the velocity and **DG–0** for the pressure.
+
+Dirichlet conditions for the velocity on the **top**, **bottom**, and **left** boundaries  
+are imposed **weakly using Nitsche’s method**  
+(see the tutorial on  
+[Weak imposition of Dirichlet conditions for the Poisson equation](../chapter1/nitsche)  
+for details).
+
+## Key Parameters
+
+- **Velocity space:** $V_h = [\mathrm{CR}_1]^2$ (nonconforming, facet-based)  
+- **Pressure space:** $Q_h = \mathrm{DG}_0$  
+- **Linearization:** Picard (frozen convection); start from Stokes baseline  
+- **Stabilization:** small pressure mass $\varepsilon_p \int_\Omega p\,q\,\mathrm{d}x$ removes the constant-pressure nullspace  
+- **Solver:** PETSc GMRES + Jacobi  
+- **Output:** {py:class}`VTKWriter<dolfinx.io.VTXWriter>` for ParaView
+
+## Governing Equations
+
+\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}
+
+### Boundary Conditions
+
+- Inlet (left): $\mathbf{u} = \mathbf{u}_{\text{in}}$  
+- Walls (top/bottom): $\mathbf{u} = \mathbf{0}$  
+- Outlet (right): natural (traction-free)
+
+## Weak Formulation
+
+For test functions $\mathbf{v}, q$, find $(\mathbf{u},p)$ such that
+
+\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}
+
+Right-hand side:
+\[
+\ell(\mathbf{v}) = (\mathbf{f},\mathbf{v})_\Omega.
+\]
+
+> **Reference**  
+> The CR–P0 mixed finite element formulation follows the framework in  
+> *F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, 1991.*  
+> The weak enforcement of velocity Dirichlet conditions is based on the symmetric  
+> Nitsche method (*J. Nitsche, 1971*).
+
+### Nitsche Boundary Terms
+
+\[
+\begin{aligned}
+& -\nu\int_{\Gamma_D} (\nabla\mathbf{u}\,\mathbf{n})\cdot\mathbf{v}\,\mathrm{d}s
+ -\nu\int_{\Gamma_D} (\nabla\mathbf{v}\,\mathbf{n})\cdot(\mathbf{u}-\mathbf{u}_D)\,\mathrm{d}s
+ + \alpha\frac{\nu}{h}\int_{\Gamma_D} (\mathbf{u}-\mathbf{u}_D)\cdot\mathbf{v}\,\mathrm{d}s \\
+&\quad - \int_{\Gamma_D} p\,\mathbf{n}\cdot\mathbf{v}\,\mathrm{d}s
+ + \int_{\Gamma_D} q\,\mathbf{n}\cdot(\mathbf{u}-\mathbf{u}_D)\,\mathrm{d}s.
+\end{aligned}
+\]
+
+### Under-relaxed Picard Update
+
+\[
+\mathbf{u}_{k+1}^{(\text{lin})}
+= \theta\,\mathbf{u}_{k}^{(\text{new})}
+ + (1-\theta)\,\mathbf{u}_{k}^{(\text{lin})}, \quad 0<\theta\le1.
+\]
+
+## Example Output (ParaView)
+
+![Velocity field](open_cavity_velocity_glyphs.png)
+
+**Figure:** Velocity magnitude field with CR–P0 (Stokes baseline).  
+Color shows $|\mathbf{u}|$ and arrows indicate flow direction and speed.  
+Inlet on the left, natural outlet on the right.
+
+## Run Instructions (Docker)
+
+```bash
+docker run --rm -it -v "$PWD":"$PWD" -w "$PWD" \
+  ghcr.io/fenics/dolfinx/dolfinx:stable \
+  python3 chapter2/monolithic_navierstokes_crp0.py \
+  --uin 0.5 --nu 1e-2 --theta 0.5 --alpha 300

chapter2/monolithic_navierstokes_crp0.py

@@ -0,0 +1,272 @@
+from __future__ import annotations
+import argparse
+import numpy as np
+from mpi4py import MPI
+from dolfinx import fem, io, mesh
+from dolfinx.fem import petsc
+from petsc4py import PETSc
+import ufl
+from basix.ufl import element, mixed_element
+
+
+def parse_args():
+    p = argparse.ArgumentParser(
+        description="CR–P0 steady Navier–Stokes with Nitsche BCs (monolithic)"
+    )
+    p.add_argument("--nx", type=int, default=64)
+    p.add_argument("--ny", type=int, default=32)
+    p.add_argument("--Lx", type=float, default=2.0)
+    p.add_argument("--Ly", type=float, default=1.0)
+    p.add_argument("--nu", type=float, default=1e-2)
+    p.add_argument("--uin", type=float, default=1.0)
+    p.add_argument("--picard-it", type=int, default=20)
+    p.add_argument("--picard-tol", type=float, default=1e-8)
+    p.add_argument("--theta", type=float, default=0.5)
+    p.add_argument(
+        "--eps-p",
+        type=float,
+        default=1e-12,
+        help="Tiny pressure mass to remove nullspace",
+    )
+    p.add_argument("--stokes-only", action="store_true")
+    p.add_argument("--no-output", action="store_true")
+    p.add_argument(
+        "--alpha",
+        type=float,
+        default=400.0,
+        help="Nitsche penalty (100–1000 recommended)",
+    )
+    # Outfile stem; we append _u.bp / _p.bp
+    p.add_argument(
+        "--outfile",
+        type=str,
+        default="chapter2/open_cavity_crp0",
+        help="Output file stem (without extension)",
+    )
+    return p.parse_args()
+
+
+args = parse_args()
+comm = MPI.COMM_WORLD
+rank = comm.rank
+
+# --- Mesh --------------------------------------------------------------
+domain = mesh.create_rectangle(
+    comm,
+    [np.array([0.0, 0.0]), np.array([args.Lx, args.Ly])],
+    [args.nx, args.ny],
+    cell_type=mesh.CellType.triangle,
+)
+tdim = domain.topology.dim
+fdim = tdim - 1
+gdim = domain.geometry.dim
+cell = domain.basix_cell()
+
+domain.topology.create_connectivity(fdim, tdim)
+domain.topology.create_connectivity(tdim, tdim)
+
+# --- Boundary facet tags (for ds & Nitsche) ---------------------------
+left_facets = mesh.locate_entities_boundary(
+    domain, fdim, lambda x: np.isclose(x[0], 0.0)
+)
+right_facets = mesh.locate_entities_boundary(
+    domain, fdim, lambda x: np.isclose(x[0], args.Lx)
+)
+bottom_facets = mesh.locate_entities_boundary(
+    domain, fdim, lambda x: np.isclose(x[1], 0.0)
+)
+top_facets = mesh.locate_entities_boundary(
+    domain, fdim, lambda x: np.isclose(x[1], args.Ly)
+)
+
+# Tag: 1=left, 2=right, 3=bottom, 4=top
+facet_indices = np.concatenate(
+    [left_facets, right_facets, bottom_facets, top_facets]
+).astype(np.int32)
+facet_tags = np.concatenate(
+    [
+        np.full_like(left_facets, 1, dtype=np.int32),
+        np.full_like(right_facets, 2, dtype=np.int32),
+        np.full_like(bottom_facets, 3, dtype=np.int32),
+        np.full_like(top_facets, 4, dtype=np.int32),
+    ]
+)
+
+if facet_indices.size == 0:
+    ft = mesh.meshtags(
+        domain,
+        fdim,
+        np.array([], dtype=np.int32),
+        np.array([], dtype=np.int32),
+    )
+else:
+    order = np.argsort(facet_indices)
+    ft = mesh.meshtags(domain, fdim, facet_indices[order], facet_tags[order])
+
+dx = ufl.Measure("dx", domain=domain)
+ds = ufl.Measure("ds", domain=domain, subdomain_data=ft)
+
+# --- Spaces: CR–P0 (mixed) --------------------------------------------
+V_el = element("CR", cell, 1, shape=(gdim,))
+Q_el = element("DG", cell, 0)
+W = fem.functionspace(domain, mixed_element([V_el, Q_el]))
+(u, p) = ufl.TrialFunctions(W)
+(v, q) = ufl.TestFunctions(W)
+
+# --- Parameters / RHS -------------------------------------------------
+nu = fem.Constant(domain, PETSc.ScalarType(args.nu))
+eps_p = fem.Constant(domain, PETSc.ScalarType(args.eps_p))
+zero = fem.Constant(domain, PETSc.ScalarType(0.0))
+f_vec = ufl.as_vector((zero, zero))  # zero body force
+
+# --- Picard state (for convection) -----------------------------------
+w = fem.Function(W)  # holds (u_k, p_k)
+u_k, p_k = w.split()
+
+
+def build_forms():
+    n = ufl.FacetNormal(domain)
+    h = ufl.CellDiameter(domain)
+    alpha = PETSc.ScalarType(args.alpha)
+
+    u_in = ufl.as_vector(
+        (PETSc.ScalarType(args.uin), PETSc.ScalarType(0.0))
+    )
+    zero_vec = ufl.as_vector(
+        (PETSc.ScalarType(0.0), PETSc.ScalarType(0.0))
+    )
+
+    # Core Stokes + tiny pressure mass
+    F = (
+        nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
+        - ufl.div(v) * p * dx
+        + q * ufl.div(u) * dx
+        + eps_p * p * q * dx
+        - ufl.inner(f_vec, v) * dx
+    )
+
+    # Convection (Picard) if not Stokes-only
+    if not args.stokes_only:
+        F += ufl.inner(ufl.dot(u_k, ufl.nabla_grad(u)), v) * dx
+
+    # --- Nitsche / consistency terms on Dirichlet boundaries ---------
+    # Tag 1: left inlet, u = u_in
+    F += (
+        -nu * ufl.inner(ufl.grad(u) * n, v)
+        -nu * ufl.inner(ufl.grad(v) * n, (u - u_in))
+        + alpha * nu / h * ufl.inner(u - u_in, v)
+        - ufl.inner(p * n, v)
+        + q * ufl.dot(u - u_in, n)
+    ) * ds(1)
+
+    # Tags 3 and 4: bottom + top, u = 0
+    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.dot(u - zero_vec, n)
+        ) * ds(tag)
+
+    # Tag 2 (right) kept as natural outlet (do-nothing)
+    a = ufl.lhs(F)
+    L = ufl.rhs(F)
+    return a, L
+
+
+def solve_once():
+    a, L = build_forms()
+    # No strong BCs with CR → Nitsche handles boundaries
+    problem = petsc.LinearProblem(
+        a,
+        L,
+        u=w,
+        bcs=[],
+        petsc_options={
+            "ksp_type": "gmres",
+            "pc_type": "jacobi",
+            "ksp_rtol": 1.0e-8,
+            "ksp_max_it": 1000,
+        },
+        petsc_options_prefix="ns_",
+    )
+    problem.solve()
+
+
+# --- Picard loop ------------------------------------------------------
+theta = float(args.theta)
+tol = float(args.picard_tol)
+max_it = int(args.picard_it)
+
+V_sub = W.sub(0)
+V0, _ = V_sub.collapse()
+u_prev_fun = fem.Function(V0)
+u_prev_arr = np.zeros_like(u_prev_fun.x.array)
+
+for it in range(1, max_it + 1):
+    solve_once()
+
+    # Velocity from mixed solution
+    u_view = w.sub(0).collapse()
+    u_curr = u_view[0] if isinstance(u_view, tuple) else u_view
+    u_curr_arr = u_curr.x.array.copy()
+
+    diff = u_curr_arr - u_prev_arr
+    err = float(np.linalg.norm(diff)) if np.all(
+        np.isfinite(diff)
+    ) else float("inf")
+
+    if rank == 0:
+        PETSc.Sys.Print(f"Picard {it:02d}: ||u - u_prev|| = {err:.3e}")
+        PETSc.Sys.Print(
+            f"  |u| range: [{np.abs(u_curr_arr).min():.3e}, "
+            f"{np.abs(u_curr_arr).max():.3e}]"
+        )
+
+    if not np.isfinite(err) or err < tol or args.stokes_only:
+        break
+
+    # Under-relax: u_k := θ u_curr + (1-θ) u_prev
+    relaxed = theta * u_curr_arr + (1.0 - theta) * u_prev_arr
+    u_k.x.array[: len(relaxed)] = relaxed
+    u_prev_arr = u_curr_arr.copy()
+
+# --- Output (BP4 / VTKWriter) ----------------------------------------
+if not args.no_output:
+    # Extract velocity and pressure from mixed solution
+    u_view = w.sub(0).collapse()
+    p_view = w.sub(1).collapse()
+    u_fun = u_view[0] if isinstance(u_view, tuple) else u_view
+    p_fun = p_view[0] if isinstance(p_view, tuple) else p_view
+
+    # Interpolate to CG1 for smoother visualization
+    Vout_u = fem.functionspace(
+        domain, element("Lagrange", cell, 1, shape=(gdim,))
+    )
+    Vout_p = fem.functionspace(domain, element("Lagrange", cell, 1))
+
+    u_out = fem.Function(Vout_u, name="u")
+    u_out.interpolate(u_fun)
+
+    p_out = fem.Function(Vout_p, name="p")
+    p_out.interpolate(p_fun)
+
+    outfile = args.outfile
+
+    # Write velocity
+    with io.VTXWriter(
+        domain.comm, outfile + "_u.bp", u_out, engine="BP4"
+    ) as vtx_u:
+        vtx_u.write(0.0)
+
+    # Write pressure
+    with io.VTXWriter(
+        domain.comm, outfile + "_p.bp", p_out, engine="BP4"
+    ) as vtx_p:
+        vtx_p.write(0.0)
+
+    if rank == 0:
+        PETSc.Sys.Print(
+            f"Wrote {outfile}_u.bp and {outfile}_p.bp"
+        )