From fb40d690ba8de962ef2eb3049c9ae557ccb7e43a Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Wed, 19 Nov 2025 18:04:08 -0800 Subject: [PATCH 01/12] =?UTF-8?q?WIP:=20add=20monolithic=20CR=E2=80=93P0?= =?UTF-8?q?=20Navier=E2=80=93Stokes=20solver?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- chapter2/monolithic_navierstokes_crp0.py | 219 +++++++++++++++++++++++ 1 file changed, 219 insertions(+) create mode 100644 chapter2/monolithic_navierstokes_crp0.py diff --git a/chapter2/monolithic_navierstokes_crp0.py b/chapter2/monolithic_navierstokes_crp0.py new file mode 100644 index 00000000..a423c680 --- /dev/null +++ b/chapter2/monolithic_navierstokes_crp0.py @@ -0,0 +1,219 @@ +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 (robust)") + 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) # 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)") + p.add_argument("--outfile", type=str, default="chapter2/open_cavity_crp0.xdmf") + 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() + +# Ensure needed connectivities +domain.topology.create_connectivity(fdim, tdim) +domain.topology.create_connectivity(tdim, tdim) + +# --- boundary facet tags (needed for ds with 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), +]) + +# Create facet meshtags +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 + ) + + # Add convection if not stokes-only + if not args.stokes_only: + F += ufl.inner(ufl.dot(u_k, ufl.nabla_grad(u)), v) * dx + + # -------- Nitsche + boundary consistency terms on Dirichlet parts -------- + # Left inlet (tag=1): u = u_in + 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) + + # Bottom (3) and Top (4): no-slip 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.inner(u - zero_vec, n) + ) * ds(tag) + + # Right (2) left as natural outlet (do-nothing) + return ufl.lhs(F), ufl.rhs(F) + + +def solve_once(): + a, L = build_forms() + # No strong BCs with CR → Nitsche handles boundaries + A = petsc.assemble_matrix(fem.form(a), bcs=[]) + A.assemble() + b = A.createVecRight(); b.set(0.0) + petsc.assemble_vector(b, fem.form(L)) + # (No lifting/BC since bcs=[]) + b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) + + if rank == 0: + PETSc.Sys.Print(f" ||b||_2 = {b.norm():.3e}") + + # Robust Krylov (no factorization) + ksp = PETSc.KSP().create(comm) + ksp.setOperators(A) + ksp.setType("gmres") + ksp.setTolerances(rtol=1e-8, max_it=1000) + ksp.getPC().setType("jacobi") + + # wrap dolfinx vector, solve, copy back + x = PETSc.Vec().createWithArray(w.x.array, comm=comm) + ksp.solve(b, x) + w.x.array[:] = x.getArray() + w.x.scatter_forward() + if rank == 0: + PETSc.Sys.Print(f" ||x||_2 = {x.norm():.3e}") + + +# --- Picard loop --- +theta = float(args.theta) +tol = float(args.picard_tol) +max_it = int(args.picard_it) + +# for residual check (velocity only) +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}, {np.abs(u_curr_arr).max():.3e}]") + + if not np.isfinite(err) or err < tol or args.stokes_only: + break + + # under-relax linearization state 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 --- +if not args.no_output: + 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 cleaner XDMF viewing + Vout = fem.functionspace(domain, element("Lagrange", cell, 1, shape=(gdim,))) + Qout = fem.functionspace(domain, element("Lagrange", cell, 1)) + u_out = fem.Function(Vout, name="u"); u_out.interpolate(u_fun) + p_out = fem.Function(Qout, name="p"); p_out.interpolate(p_fun) + + with io.XDMFFile(domain.comm, args.outfile, "w") as xdmf: + xdmf.write_mesh(domain) + xdmf.write_function(u_out) + xdmf.write_function(p_out) + if rank == 0: + PETSc.Sys.Print(f"Wrote {args.outfile}") From 2dddd178ab83a0cd1391bc4b5330e6906f2e03c7 Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Tue, 25 Nov 2025 00:55:05 -0800 Subject: [PATCH 02/12] =?UTF-8?q?Add=20CR=E2=80=93P0=20Navier=E2=80=93Stok?= =?UTF-8?q?es=20tutorial=20markdown=20file?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- chapter2/monolithic_navierstokes_crp0.md | 103 +++++++++++++++++++++++ 1 file changed, 103 insertions(+) create mode 100644 chapter2/monolithic_navierstokes_crp0.md diff --git a/chapter2/monolithic_navierstokes_crp0.md b/chapter2/monolithic_navierstokes_crp0.md new file mode 100644 index 00000000..68d91405 --- /dev/null +++ b/chapter2/monolithic_navierstokes_crp0.md @@ -0,0 +1,103 @@ +# Monolithic Incompressible Navier–Stokes with CR–P0 (Steady Picard) + +--- + +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. + +--- + +## 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\,dx $ removes nullspace +- **Solver:** PETSc GMRES + Jacobi +- **Output:** XDMF for ParaView + +--- + +## Governing Equations + +$$ +-\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. +$$ + +### 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: + +$$ +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. +$$ + +Right-hand side: +$$ +\ell(\mathbf{v}) = (\mathbf{f},\mathbf{v})_\Omega. +$$ + +--- + +### Nitsche Boundary Terms + +$$ +\begin{aligned} +& -\nu\int_{\Gamma_D} (\nabla\mathbf{u}\,\mathbf{n})\cdot\mathbf{v}\,ds + -\nu\int_{\Gamma_D} (\nabla\mathbf{v}\,\mathbf{n})\cdot(\mathbf{u}-\mathbf{u}_D)\,ds + + \alpha\frac{\nu}{h}\int_{\Gamma_D} (\mathbf{u}-\mathbf{u}_D)\cdot\mathbf{v}\,ds \\ +&\quad - \int_{\Gamma_D} p\,\mathbf{n}\cdot\mathbf{v}\,ds + + \int_{\Gamma_D} q\,\mathbf{n}\cdot(\mathbf{u}-\mathbf{u}_D)\,ds. +\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](chapter2/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 From 526a5556b99bd5bbfeb481df30faf555190e40cd Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Sat, 6 Dec 2025 00:31:20 -0800 Subject: [PATCH 03/12] Update tutorial .md with corrected image path and equations --- chapter2/monolithic_navierstokes_crp0.md | 70 +++++++++++------------- 1 file changed, 32 insertions(+), 38 deletions(-) diff --git a/chapter2/monolithic_navierstokes_crp0.md b/chapter2/monolithic_navierstokes_crp0.md index 68d91405..dcdcfa8f 100644 --- a/chapter2/monolithic_navierstokes_crp0.md +++ b/chapter2/monolithic_navierstokes_crp0.md @@ -1,66 +1,64 @@ # Monolithic Incompressible Navier–Stokes with CR–P0 (Steady Picard) ---- +We solve the steady incompressible Navier–Stokes equations in a 2D channel. -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 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\,dx $ removes nullspace - **Solver:** PETSc GMRES + Jacobi -- **Output:** XDMF for ParaView - ---- +- **Output:** {py:class}`VTKWriter` for ParaView ## Governing Equations -$$ --\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} ### Boundary Conditions -- Inlet (left): $ \mathbf{u} = \mathbf{u}_{\text{in}} $ -- Walls (top/bottom): $ \mathbf{u} = \mathbf{0} $ -- Outlet (right): natural (traction-free) ---- +- 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: -$$ -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, -$$ - -$$ +\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, -$$ - -$$ + &= ((\mathbf{u}_k\cdot\nabla)\mathbf{u},\mathbf{v})_\Omega,\\ a_{\text{pm}}(p;q) - = \varepsilon_p (p,q)_\Omega. -$$ + &= \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 @@ -82,18 +80,14 @@ $$ + (1-\theta)\,\mathbf{u}_{k}^{(\text{lin})}, \quad 0<\theta\le1. $$ ---- - ## Example Output (ParaView) -![Velocity field](chapter2/open_cavity_velocity_glyphs.png) +![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 From 9adcad80d71d1927733303b22877cb3715ca6147 Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Sat, 6 Dec 2025 01:46:45 -0800 Subject: [PATCH 04/12] =?UTF-8?q?Tidy=20CR=E2=80=93P0=20Navier=E2=80=93Sto?= =?UTF-8?q?kes=20tutorial=20text?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- chapter2/monolithic_navierstokes_crp0.md | 38 ++++++++++++------------ 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/chapter2/monolithic_navierstokes_crp0.md b/chapter2/monolithic_navierstokes_crp0.md index dcdcfa8f..d0785d9c 100644 --- a/chapter2/monolithic_navierstokes_crp0.md +++ b/chapter2/monolithic_navierstokes_crp0.md @@ -14,10 +14,10 @@ for details). ## Key Parameters -- **Velocity space:** $ V_h = [\mathrm{CR}_1]^2 $ (nonconforming, facet-based) -- **Pressure space:** $ Q_h = \mathrm{DG}_0 $ +- **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\,dx $ removes nullspace +- **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` for ParaView @@ -30,13 +30,13 @@ for details). ### Boundary Conditions -- Inlet (left): $ \mathbf{u} = \mathbf{u}_{\text{in}} $ -- Walls (top/bottom): $ \mathbf{u} = \mathbf{0} $ +- 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: +For test functions $\mathbf{v}, q$, find $(\mathbf{u},p)$ such that \begin{align} a_{\text{core}}((\mathbf{u},p);(\mathbf{v},q)) @@ -50,11 +50,11 @@ a_{\text{pm}}(p;q) \end{align} Right-hand side: -$$ +\[ \ell(\mathbf{v}) = (\mathbf{f},\mathbf{v})_\Omega. -$$ +\] -> **Reference:** +> **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 @@ -62,30 +62,30 @@ $$ ### Nitsche Boundary Terms -$$ +\[ \begin{aligned} -& -\nu\int_{\Gamma_D} (\nabla\mathbf{u}\,\mathbf{n})\cdot\mathbf{v}\,ds - -\nu\int_{\Gamma_D} (\nabla\mathbf{v}\,\mathbf{n})\cdot(\mathbf{u}-\mathbf{u}_D)\,ds - + \alpha\frac{\nu}{h}\int_{\Gamma_D} (\mathbf{u}-\mathbf{u}_D)\cdot\mathbf{v}\,ds \\ -&\quad - \int_{\Gamma_D} p\,\mathbf{n}\cdot\mathbf{v}\,ds - + \int_{\Gamma_D} q\,\mathbf{n}\cdot(\mathbf{u}-\mathbf{u}_D)\,ds. +& -\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. +Color shows $|\mathbf{u}|$ and arrows indicate flow direction and speed. Inlet on the left, natural outlet on the right. ## Run Instructions (Docker) From aa123d32b2c1a48c41e5b486f8560a5c7d0b7e81 Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Sat, 6 Dec 2025 15:12:01 -0800 Subject: [PATCH 05/12] Update chapter2/monolithic_navierstokes_crp0.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jørgen Schartum Dokken --- chapter2/monolithic_navierstokes_crp0.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/chapter2/monolithic_navierstokes_crp0.py b/chapter2/monolithic_navierstokes_crp0.py index a423c680..3fe29232 100644 --- a/chapter2/monolithic_navierstokes_crp0.py +++ b/chapter2/monolithic_navierstokes_crp0.py @@ -206,14 +206,12 @@ def solve_once(): p_fun = p_view[0] if isinstance(p_view, tuple) else p_view # interpolate to CG1 for cleaner XDMF viewing - Vout = fem.functionspace(domain, element("Lagrange", cell, 1, shape=(gdim,))) + Vout = fem.functionspace(domain, element("DG", cell, 1, shape=(gdim,))) Qout = fem.functionspace(domain, element("Lagrange", cell, 1)) - u_out = fem.Function(Vout, name="u"); u_out.interpolate(u_fun) - p_out = fem.Function(Qout, name="p"); p_out.interpolate(p_fun) + u_out = fem.Function(Vout, name="u") + u_out.interpolate(u_fun) - with io.XDMFFile(domain.comm, args.outfile, "w") as xdmf: - xdmf.write_mesh(domain) - xdmf.write_function(u_out) - xdmf.write_function(p_out) + with io.VTXWriter(domain.comm, args.outfile, [u_out, p_fun]) as writer: + writer.write(0.0) if rank == 0: PETSc.Sys.Print(f"Wrote {args.outfile}") From 72ccb23cac89f92070c72070ec1a457d8a54f6b6 Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Sat, 6 Dec 2025 15:12:23 -0800 Subject: [PATCH 06/12] Update chapter2/monolithic_navierstokes_crp0.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jørgen Schartum Dokken --- chapter2/monolithic_navierstokes_crp0.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/chapter2/monolithic_navierstokes_crp0.py b/chapter2/monolithic_navierstokes_crp0.py index 3fe29232..5d7647bb 100644 --- a/chapter2/monolithic_navierstokes_crp0.py +++ b/chapter2/monolithic_navierstokes_crp0.py @@ -122,14 +122,13 @@ def build_forms(): ) * ds(1) # Bottom (3) and Top (4): no-slip 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.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)) # Right (2) left as natural outlet (do-nothing) return ufl.lhs(F), ufl.rhs(F) From 71ff533d58a1580596834c2f9d5572f4902994f8 Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Sat, 6 Dec 2025 15:13:20 -0800 Subject: [PATCH 07/12] Update chapter2/monolithic_navierstokes_crp0.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jørgen Schartum Dokken --- chapter2/monolithic_navierstokes_crp0.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/chapter2/monolithic_navierstokes_crp0.py b/chapter2/monolithic_navierstokes_crp0.py index 5d7647bb..369c1cea 100644 --- a/chapter2/monolithic_navierstokes_crp0.py +++ b/chapter2/monolithic_navierstokes_crp0.py @@ -83,8 +83,7 @@ def parse_args(): # --- 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 +f_vec = fem.Constant(domain, np.zeros(mesh.geometry.dim, dtype=dolfinx.default_scalar_type)) # --- Picard state (for convection) --- w = fem.Function(W) # holds (u_k, p_k) From 7a18b68cc235341355e4146ec6aed8085e4f4953 Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Sat, 6 Dec 2025 15:14:03 -0800 Subject: [PATCH 08/12] Update chapter2/monolithic_navierstokes_crp0.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jørgen Schartum Dokken --- chapter2/monolithic_navierstokes_crp0.py | 1 + 1 file changed, 1 insertion(+) diff --git a/chapter2/monolithic_navierstokes_crp0.py b/chapter2/monolithic_navierstokes_crp0.py index 369c1cea..fb770de6 100644 --- a/chapter2/monolithic_navierstokes_crp0.py +++ b/chapter2/monolithic_navierstokes_crp0.py @@ -152,6 +152,7 @@ def solve_once(): ksp.setType("gmres") ksp.setTolerances(rtol=1e-8, max_it=1000) ksp.getPC().setType("jacobi") + ksp.setErrorIfNotConverged(True) # wrap dolfinx vector, solve, copy back x = PETSc.Vec().createWithArray(w.x.array, comm=comm) From 5341048f539dadd268f9ad5413e61fca544ee220 Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Sat, 6 Dec 2025 15:14:56 -0800 Subject: [PATCH 09/12] Update chapter2/monolithic_navierstokes_crp0.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jørgen Schartum Dokken --- chapter2/monolithic_navierstokes_crp0.py | 27 ++++++++++++------------ 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/chapter2/monolithic_navierstokes_crp0.py b/chapter2/monolithic_navierstokes_crp0.py index fb770de6..bb53d8b2 100644 --- a/chapter2/monolithic_navierstokes_crp0.py +++ b/chapter2/monolithic_navierstokes_crp0.py @@ -55,20 +55,19 @@ def parse_args(): 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), -]) - -# Create facet meshtags -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]) +num_facets_local = ( + domain.topology.index_map(fdim).size_local + + domain.topology.index_map(fdim).num_ghosts +) +facet_markers = np.full(num_facets_local, 0, dtype=np.int32) +facet_markers[left_facets] = 1 +facet_markers[right_facets] = 2 +facet_markers[bottom_facets] = 3 +facet_markers[top_facets] = 4 +facet_indices = np.flatnonzero(facet_markers) + +ft = mesh.meshtags(domain, fdim, facet_indices, facet_markers[facet_indices]) + dx = ufl.Measure("dx", domain=domain) ds = ufl.Measure("ds", domain=domain, subdomain_data=ft) From 97c45f6b08f2082f6f549cce2f2fe295d270a56a Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Sat, 6 Dec 2025 15:16:03 -0800 Subject: [PATCH 10/12] Update chapter2/monolithic_navierstokes_crp0.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jørgen Schartum Dokken --- chapter2/monolithic_navierstokes_crp0.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/chapter2/monolithic_navierstokes_crp0.py b/chapter2/monolithic_navierstokes_crp0.py index bb53d8b2..dfcc4804 100644 --- a/chapter2/monolithic_navierstokes_crp0.py +++ b/chapter2/monolithic_navierstokes_crp0.py @@ -80,7 +80,7 @@ def parse_args(): (v, q) = ufl.TestFunctions(W) # --- parameters / RHS --- -nu = fem.Constant(domain, PETSc.ScalarType(args.nu)) +nu = fem.Constant(domain, dolfinx.default_scalar_type(args.nu)) eps_p = fem.Constant(domain, PETSc.ScalarType(args.eps_p)) f_vec = fem.Constant(domain, np.zeros(mesh.geometry.dim, dtype=dolfinx.default_scalar_type)) From 44d23f3135b6bd64760c5dbf1563c782a0353887 Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Sat, 6 Dec 2025 15:16:30 -0800 Subject: [PATCH 11/12] Update chapter2/monolithic_navierstokes_crp0.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jørgen Schartum Dokken --- chapter2/monolithic_navierstokes_crp0.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/chapter2/monolithic_navierstokes_crp0.py b/chapter2/monolithic_navierstokes_crp0.py index dfcc4804..334d4687 100644 --- a/chapter2/monolithic_navierstokes_crp0.py +++ b/chapter2/monolithic_navierstokes_crp0.py @@ -129,7 +129,7 @@ def build_forms(): ) * ds((3,4)) # Right (2) left as natural outlet (do-nothing) - return ufl.lhs(F), ufl.rhs(F) + return ufl.system(F) def solve_once(): From 8e965a3fae98970b54bbf8b6b1d446a13bf79181 Mon Sep 17 00:00:00 2001 From: Anusha Karve Date: Sat, 6 Dec 2025 15:41:35 -0800 Subject: [PATCH 12/12] =?UTF-8?q?Apply=20reviewer=20feedback=20and=20updat?= =?UTF-8?q?e=20Navier=E2=80=93Stokes=20CR=E2=80=93P0=20solver?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- chapter2/monolithic_navierstokes_crp0.py | 235 ++++++++++++++--------- 1 file changed, 146 insertions(+), 89 deletions(-) diff --git a/chapter2/monolithic_navierstokes_crp0.py b/chapter2/monolithic_navierstokes_crp0.py index 334d4687..8221fe9a 100644 --- a/chapter2/monolithic_navierstokes_crp0.py +++ b/chapter2/monolithic_navierstokes_crp0.py @@ -10,21 +10,39 @@ def parse_args(): - p = argparse.ArgumentParser(description="CR–P0 steady Navier–Stokes with Nitsche BCs (robust)") + 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("--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) # tiny pressure mass to remove nullspace + 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)") - p.add_argument("--outfile", type=str, default="chapter2/open_cavity_crp0.xdmf") + 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() @@ -32,7 +50,7 @@ def parse_args(): comm = MPI.COMM_WORLD rank = comm.rank -# --- mesh --- +# --- Mesh -------------------------------------------------------------- domain = mesh.create_rectangle( comm, [np.array([0.0, 0.0]), np.array([args.Lx, args.Ly])], @@ -44,57 +62,79 @@ def parse_args(): gdim = domain.geometry.dim cell = domain.basix_cell() -# Ensure needed connectivities domain.topology.create_connectivity(fdim, tdim) domain.topology.create_connectivity(tdim, tdim) -# --- boundary facet tags (needed for ds with 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)) +# --- 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 -num_facets_local = ( - domain.topology.index_map(fdim).size_local - + domain.topology.index_map(fdim).num_ghosts +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), + ] ) -facet_markers = np.full(num_facets_local, 0, dtype=np.int32) -facet_markers[left_facets] = 1 -facet_markers[right_facets] = 2 -facet_markers[bottom_facets] = 3 -facet_markers[top_facets] = 4 -facet_indices = np.flatnonzero(facet_markers) - -ft = mesh.meshtags(domain, fdim, facet_indices, facet_markers[facet_indices]) +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) --- +# --- 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, dolfinx.default_scalar_type(args.nu)) +# --- Parameters / RHS ------------------------------------------------- +nu = fem.Constant(domain, PETSc.ScalarType(args.nu)) eps_p = fem.Constant(domain, PETSc.ScalarType(args.eps_p)) -f_vec = fem.Constant(domain, np.zeros(mesh.geometry.dim, dtype=dolfinx.default_scalar_type)) +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) +# --- 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))) + 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 = ( @@ -105,69 +145,60 @@ def build_forms(): - ufl.inner(f_vec, v) * dx ) - # Add convection if not stokes-only + # 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 + boundary consistency terms on Dirichlet parts -------- - # Left inlet (tag=1): u = u_in + # --- Nitsche / consistency terms on Dirichlet boundaries --------- + # Tag 1: left inlet, u = u_in 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) + -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) - # Bottom (3) and Top (4): no-slip u = 0 - 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)) + # 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) - # Right (2) left as natural outlet (do-nothing) - return ufl.system(F) + # 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 - A = petsc.assemble_matrix(fem.form(a), bcs=[]) - A.assemble() - b = A.createVecRight(); b.set(0.0) - petsc.assemble_vector(b, fem.form(L)) - # (No lifting/BC since bcs=[]) - b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) - - if rank == 0: - PETSc.Sys.Print(f" ||b||_2 = {b.norm():.3e}") - - # Robust Krylov (no factorization) - ksp = PETSc.KSP().create(comm) - ksp.setOperators(A) - ksp.setType("gmres") - ksp.setTolerances(rtol=1e-8, max_it=1000) - ksp.getPC().setType("jacobi") - ksp.setErrorIfNotConverged(True) - - # wrap dolfinx vector, solve, copy back - x = PETSc.Vec().createWithArray(w.x.array, comm=comm) - ksp.solve(b, x) - w.x.array[:] = x.getArray() - w.x.scatter_forward() - if rank == 0: - PETSc.Sys.Print(f" ||x||_2 = {x.norm():.3e}") + 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 --- +# --- Picard loop ------------------------------------------------------ theta = float(args.theta) tol = float(args.picard_tol) max_it = int(args.picard_it) -# for residual check (velocity only) V_sub = W.sub(0) V0, _ = V_sub.collapse() u_prev_fun = fem.Function(V0) @@ -176,40 +207,66 @@ def solve_once(): for it in range(1, max_it + 1): solve_once() - # velocity from mixed solution + # 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") + 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}, {np.abs(u_curr_arr).max():.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 linearization state u_k := θ u_curr + (1-θ) u_prev + # 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_k.x.array[: len(relaxed)] = relaxed u_prev_arr = u_curr_arr.copy() -# --- output --- +# --- 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 cleaner XDMF viewing - Vout = fem.functionspace(domain, element("DG", cell, 1, shape=(gdim,))) - Qout = fem.functionspace(domain, element("Lagrange", cell, 1)) - u_out = fem.Function(Vout, name="u") + # 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) - with io.VTXWriter(domain.comm, args.outfile, [u_out, p_fun]) as writer: - writer.write(0.0) + 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 {args.outfile}") + PETSc.Sys.Print( + f"Wrote {outfile}_u.bp and {outfile}_p.bp" + )