{ "cells": [ { "cell_type": "markdown", "id": "5cdf0eee", "metadata": {}, "source": [ "# L3B -- Hybridization and postprocessing\n", "\n", "Companion notebook to Lecture 3, second demonstration block (`Cell N` labels\n", "are cited on the slides). Tags as in L3A.\n", "\n", "Setting as in L3A: mixed Laplacian, unit square, natural boundary condition\n", "(no `DirichletBC` on the mixed problem; the trace space in cell 8 does carry\n", "one, because the multiplier represents the Dirichlet datum -- slide 27).\n", "\n", "> **Warning (Firedrake / FEniCS indexing).** Firedrake indexes the\n", "> Raviart--Thomas family shifted by one relative to the course notation\n", "> (Boffi--Brezzi--Fortin): `FunctionSpace(mesh, \"RT\", 1)` is $RT_0$.\n", "> The $BDM$ indices agree: `(\"BDM\", 1)` is $BDM_1$. Cross-check\n", "> element definitions on [DefElement](https://defelement.org).\n", "\n", "Provenance: expected numbers in the comments are pinned by\n", "`verify_rt0_2d.py` (equivalence to machine precision; condensed system\n", "symmetric positive definite, smallest eigenvalue 0.102 at $N=8$; multiplier\n", "second-order at edge midpoints). Timings and iteration counts are\n", "machine-dependent and are produced live, not quoted on the slides." ] }, { "cell_type": "code", "execution_count": null, "id": "5445fddd", "metadata": {}, "outputs": [], "source": [ "# Cell 1 [SETUP] -- Firedrake import guard (FEM on Colab).\n", "# On a local Firedrake installation this cell only imports and reports.\n", "try:\n", " import firedrake\n", "except ImportError:\n", " try:\n", " import google.colab # noqa: F401\n", " !wget \"https://fem-on-colab.github.io/releases/firedrake-install-release-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n", " import firedrake\n", " except ImportError:\n", " raise RuntimeError(\n", " \"Firedrake is not available. Install it locally or run this \"\n", " \"notebook on Google Colab (the cell above then installs it).\")\n", "try:\n", " from importlib.metadata import version\n", " print(\"Firedrake version:\", version(\"firedrake\"))\n", "except Exception:\n", " print(\"Firedrake imported (version string not available).\")\n", "print(\"Firedrake location:\", firedrake.__file__)" ] }, { "cell_type": "code", "execution_count": null, "id": "6d1491f3", "metadata": {}, "outputs": [], "source": [ "# Cell 2 [LECTURE] -- the monolithic solve (baseline).\n", "import time\n", "import numpy as np\n", "from firedrake import *\n", "import warnings\n", "warnings.filterwarnings(\n", " \"ignore\", message=\".*ndarray.shape.*\", category=DeprecationWarning)\n", "# \"from firedrake import *\" shadows the standard-library logging module;\n", "# use an alias (see L3A, cell 2, for the full rationale).\n", "import logging as pylogging\n", "pylogging.getLogger(\"tsfc\").setLevel(pylogging.ERROR)\n", "\n", "N = 64\n", "mesh = UnitSquareMesh(N, N)\n", "V = FunctionSpace(mesh, \"RT\", 1) # RT_0 (see the header warning)\n", "Q = FunctionSpace(mesh, \"DG\", 0)\n", "W = V * Q\n", "x, y = SpatialCoordinate(mesh)\n", "u_ex = sin(pi*x)*sin(pi*y)\n", "f = 2*pi**2*sin(pi*x)*sin(pi*y)\n", "sigma, u = TrialFunctions(W)\n", "tau, v = TestFunctions(W)\n", "a = (inner(sigma, tau) + div(tau)*u + div(sigma)*v)*dx\n", "L = -f*v*dx\n", "\n", "# Each configuration is solved twice: the first call includes Firedrake's\n", "# just-in-time kernel compilation; the second is the honest solver time.\n", "w_mono = Function(W)\n", "t_mono = []\n", "for _ in range(2):\n", " w_mono.assign(0)\n", " t0 = time.perf_counter()\n", " solve(a == L, w_mono, solver_parameters={\n", " \"ksp_type\": \"preonly\", \"pc_type\": \"lu\",\n", " \"pc_factor_mat_solver_type\": \"mumps\"})\n", " t_mono.append(time.perf_counter() - t0)\n", "print(f\"monolithic mixed system: {W.dim()} unknowns \"\n", " f\"({V.dim()} fluxes + {Q.dim()} scalars), MUMPS LU\")\n", "print(f\" first solve {t_mono[0]:.2f} s (includes code generation), \"\n", " f\"second {t_mono[1]:.2f} s\")\n", "sig_m, u_m = w_mono.subfunctions\n", "print(\"||u - u_h||_0 =\", sqrt(assemble((u_ex - u_m)**2*dx)))" ] }, { "cell_type": "markdown", "id": "adedba9a", "metadata": {}, "source": [ "## Demo 2a -- HybridizationPC (slide 30)\n", "\n", "The same discrete problem, solved through the hybridized route: broken fluxes, facet multipliers, elementwise condensation, and a conjugate-gradient solve of the symmetric positive definite trace system." ] }, { "cell_type": "code", "execution_count": null, "id": "af9edf88", "metadata": {}, "outputs": [], "source": [ "# Cell 4 [LECTURE] -- the same problem through HybridizationPC.\n", "# Firedrake assembles the hybridized (broken + multiplier) problem, condenses\n", "# to the trace system (Proposition 28: symmetric positive definite), solves\n", "# it by CG, and back-substitutes elementwise. Adjust the inner solver\n", "# parameters to the installed build if needed (pre-flight item).\n", "w_hyb = Function(W)\n", "params = {\n", " \"mat_type\": \"matfree\",\n", " \"ksp_type\": \"preonly\",\n", " \"pc_type\": \"python\",\n", " \"pc_python_type\": \"firedrake.HybridizationPC\",\n", " \"hybridization\": {\n", " \"ksp_type\": \"cg\",\n", " \"pc_type\": \"gamg\",\n", " \"ksp_rtol\": 1e-10,\n", " \"ksp_converged_reason\": None, # prints the CG iteration count\n", " },\n", "}\n", "t_hyb = []\n", "for _ in range(2):\n", " w_hyb.assign(0)\n", " t0 = time.perf_counter()\n", " solve(a == L, w_hyb, solver_parameters=params)\n", " t_hyb.append(time.perf_counter() - t0)\n", "sig_hy, u_hy = w_hyb.subfunctions\n", "\n", "ntrace = FunctionSpace(mesh, \"HDiv Trace\", 0).dim()\n", "nbdry = 4*N\n", "print(f\"trace unknowns: {ntrace} total facets, {ntrace - nbdry} interior \"\n", " f\"(condensed system size; cf. slide 28: 12160 at N=64)\")\n", "print(f\"hybridized solve: first {t_hyb[0]:.2f} s (code generation), \"\n", " f\"second {t_hyb[1]:.2f} s (monolithic, second: {t_mono[1]:.2f} s)\")\n", "print(\"At this size a 2D direct factorization is expected to win; the\")\n", "print(\"condensed route pays off through structure (SPD, CG/multigrid,\")\n", "print(\"memory) at large scale and in three dimensions, not at N = 64.\")\n", "print(\"consistency ||sigma_mono - sigma_hyb||_0 =\",\n", " sqrt(assemble(inner(sig_m - sig_hy, sig_m - sig_hy)*dx)))\n", "print(\"(Theorem 27: the two routes compute the same discrete solution;\")\n", "print(\" the reference implementation confirms this to 9.7e-15 at N=8.)\")" ] }, { "cell_type": "markdown", "id": "d1497803", "metadata": {}, "source": [ "## Demo 2b -- the postprocessed scalar (slides 29 and 31)" ] }, { "cell_type": "code", "execution_count": null, "id": "c306f0a5", "metadata": {}, "outputs": [], "source": [ "# Cell 6 [LECTURE] -- the postprocessed scalar u* (k = 0 construction).\n", "# On each element: grad u*|_K = elementwise mean of sigma_h (kappa = 1 here),\n", "# and the mean of u* equals u_h. Both are local; see slide 29.\n", "def postprocess(sig_h, u_h, mesh):\n", " DG1 = FunctionSpace(mesh, \"DG\", 1)\n", " VDG0 = VectorFunctionSpace(mesh, \"DG\", 0)\n", " s0 = Function(VDG0).project(sig_h) # elementwise means of sigma_h\n", " xc = Function(VDG0).interpolate(SpatialCoordinate(mesh)) # centroids\n", " xx = SpatialCoordinate(mesh)\n", " return Function(DG1).interpolate(u_h + dot(s0, xx - xc))\n", "\n", "print(\"N ||u-u_h|| ||u-u*||\")\n", "errs = []\n", "for NN in [8, 16, 32, 64]:\n", " mesh_ = UnitSquareMesh(NN, NN)\n", " V_ = FunctionSpace(mesh_, \"RT\", 1); Q_ = FunctionSpace(mesh_, \"DG\", 0)\n", " W_ = V_ * Q_\n", " x_, y_ = SpatialCoordinate(mesh_)\n", " u_e = sin(pi*x_)*sin(pi*y_)\n", " f_ = 2*pi**2*sin(pi*x_)*sin(pi*y_)\n", " s_, u_ = TrialFunctions(W_); t_, v_ = TestFunctions(W_)\n", " a_ = (inner(s_, t_) + div(t_)*u_ + div(s_)*v_)*dx\n", " w_ = Function(W_)\n", " solve(a_ == -f_*v_*dx, w_, solver_parameters={\n", " \"ksp_type\": \"preonly\", \"pc_type\": \"lu\",\n", " \"pc_factor_mat_solver_type\": \"mumps\"})\n", " sh_, uh_ = w_.subfunctions\n", " us_ = postprocess(sh_, uh_, mesh_)\n", " e1 = sqrt(assemble((u_e - uh_)**2*dx))\n", " e2 = sqrt(assemble((u_e - us_)**2*dx))\n", " errs.append((e1, e2))\n", " print(f\"{NN:<6d}{e1:14.4e}{e2:16.4e}\")\n", "errs = np.array(errs)\n", "r = np.log2(errs[:-1]/errs[1:])\n", "for rr in r:\n", " print(f\"rate {rr[0]:12.2f}{rr[1]:16.2f}\")\n", "print(\"\\nExpected (verify_rt0_2d.py): rates 1.00 and 2.00;\")\n", "print(\"||u - u*|| = 4.9e-4 at N = 32.\")" ] }, { "cell_type": "markdown", "id": "3cbd0bed", "metadata": {}, "source": [ "The multiplier $\\lambda_h$ produced inside `HybridizationPC` is internal to the preconditioner; cell 8 assembles the three-field problem explicitly with Slate, which exposes $\\lambda_h$ and reproduces the condensation of slide 28 line by line." ] }, { "cell_type": "code", "execution_count": null, "id": "56406bed", "metadata": {}, "outputs": [], "source": [ "# Cell 8 [SELF-STUDY S1 / ADD-BACK 1] -- hybridization by hand, with Slate.\n", "# The three-field problem of slide 27, condensed explicitly to the trace\n", "# system of slide 28. This cell follows the Firedrake static-condensation\n", "# demo; validate on the installed build (pre-flight item). The two-triangle\n", "# hand computation of Exercise 4 is the same algebra on one shared edge.\n", "Nh = 16\n", "mesh_h = UnitSquareMesh(Nh, Nh)\n", "el = FiniteElement(\"RT\", triangle, 1) # RT_0 in course indexing\n", "Vb = FunctionSpace(mesh_h, BrokenElement(el))\n", "Qh = FunctionSpace(mesh_h, \"DG\", 0)\n", "Th = FunctionSpace(mesh_h, \"HDiv Trace\", 0)\n", "Wh = Vb * Qh * Th\n", "xh, yh = SpatialCoordinate(mesh_h)\n", "u_e = sin(pi*xh)*sin(pi*yh)\n", "fh = 2*pi**2*sin(pi*xh)*sin(pi*yh)\n", "sigma, u, lam = TrialFunctions(Wh)\n", "tau, v, mu = TestFunctions(Wh)\n", "n = FacetNormal(mesh_h)\n", "a_h = (inner(sigma, tau) + div(tau)*u + div(sigma)*v)*dx \\\n", " - (jump(tau, n)*lam('+') + jump(sigma, n)*mu('+'))*dS\n", "L_h = -fh*v*dx\n", "# lambda = u_D = 0 on the boundary facets. The condensed operator S carries\n", "# arguments on the collapsed trace space Th, not on the indexed subspace\n", "# Wh.sub(2), so the boundary condition must be defined on Th (defining it on\n", "# Wh.sub(2) raises \"bc space does not match the test or trial function\n", "# space\" at assemble(S, bcs=bc) on current Firedrake).\n", "bc = DirichletBC(Th, 0, \"on_boundary\")\n", "\n", "A = Tensor(a_h)\n", "F = Tensor(L_h)\n", "# elementwise elimination of (sigma, u): the (0:2, 0:2) block is cell-local\n", "S = A.blocks[2, 2] - A.blocks[2, :2] * A.blocks[:2, :2].inv * A.blocks[:2, 2]\n", "E = -A.blocks[2, :2] * A.blocks[:2, :2].inv * F.blocks[:2]\n", "\n", "Smat = assemble(S, bcs=bc)\n", "Evec = assemble(E)\n", "lam_h = Function(Th)\n", "# The condensed matrix is symmetric positive definite (Proposition 28), so\n", "# at this size a direct factorization is the appropriate solver; the\n", "# iterative route (CG on the trace system) is what HybridizationPC runs in\n", "# cell 4, where Firedrake configures the preconditioner for this operator.\n", "# (An off-the-shelf \"cg\" + \"gamg\" here can diverge: the assembled matrix\n", "# carries identity rows from the boundary condition, and the outer SNES\n", "# masks the inner failure reason.)\n", "solve(Smat, lam_h, Evec, solver_parameters={\n", " \"ksp_type\": \"preonly\", \"pc_type\": \"lu\",\n", " \"pc_factor_mat_solver_type\": \"mumps\"})\n", "\n", "# back-substitution, elementwise\n", "sig_r = Function(Vb); u_r = Function(Qh)\n", "rec = A.blocks[:2, :2].inv * (F.blocks[:2]\n", " - A.blocks[:2, 2] * AssembledVector(lam_h))\n", "assemble(rec.blocks[0], tensor=sig_r)\n", "assemble(rec.blocks[1], tensor=u_r)\n", "print(\"||u_slate - u_exact||_0 =\", sqrt(assemble((u_r - u_e)**2*dx)))\n", "\n", "# Superconvergence of the multiplier: lambda_e approximates u on every\n", "# interior edge to second order. The skeleton L2 distance\n", "# sqrt(assemble((lam_h('+') - u_e('+'))**2*dS)) is NOT the right diagnostic:\n", "# it is first-order (the multiplier is constant per edge) over a skeleton of\n", "# growing total measure (about 2.5e-1 at N = 16). The second-order quantity\n", "# is the per-edge comparison. We compute it with the same dS assembly\n", "# machinery as the hybrid form: for the trace test function mu, assembling\n", "# g*mu('+')*dS returns the per-facet integrals of g; dividing by the facet\n", "# measures gives per-facet means.\n", "import numpy as np\n", "mu = TestFunction(Th)\n", "meas = assemble(Constant(1.0)*mu('+')*dS).dat.data_ro.copy() # facet lengths\n", "num = assemble((lam_h('+') - u_e('+'))*mu('+')*dS).dat.data_ro.copy()\n", "interior = meas > 1e-14 # boundary facets carry no dS entry\n", "lam_merr = np.max(np.abs(num[interior] / meas[interior]))\n", "print(f\"multiplier superconvergence: max_e |lambda_e - mean_e(u)| = \"\n", " f\"{lam_merr:.4e}\")\n", "print(\"(verify_rt0_2d.py: 4.20e-2 / 1.20e-2 / 3.14e-3 at N = 4/8/16,\")\n", "print(\" rates 1.81, 1.93; the midpoint comparison of the same record gives\")\n", "print(\" 5.81e-2 / 1.65e-2 / 4.24e-3, rates 1.82, 1.96.)\")" ] }, { "cell_type": "markdown", "id": "7247a907", "metadata": {}, "source": [ "Summary: the hybridized route computes the same $(\\boldsymbol\\sigma_h,u_h)$ as the monolithic solve (Theorem 27), through a smaller symmetric positive definite system (Proposition 28); the multiplier is a second-order approximation of $u$ on the mesh skeleton, and the postprocessed $u_h^*$ carries the scalar accuracy of a richer pair at lowest-order cost (Proposition 29)." ] } ], "metadata": { "kernelspec": { "display_name": "Firedrake", "language": "python", "name": "firedrake" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }