Skip to content

Mixed Poisson with efficient pre-conditioning #264

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jorgensd merged 4 commits into from
Aug 12, 2025
Merged

Mixed Poisson with efficient pre-conditioning #264

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7874f69
Start on mixed poisson demo
jorgensd Aug 9, 2025
920a854
Finalize implementation
jorgensd Aug 9, 2025
73a7a05
Add newline
jorgensd Aug 9, 2025
ad68e80
Remove adding fieldsplit as it is now part of main
jorgensd Aug 12, 2025
File filter

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Finalize implementation
  • Loading branch information
@jorgensd
jorgensd committed Aug 9, 2025
commit 920a854c312a5a7d2dfb401c969aeeb6e5c79bcd
1 change: 1 addition & 0 deletions _toc.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ parts:
- file: chapter3/em
- caption: Improving your FEniCSx code
chapters:
- file: chapter4/mixed_poisson
- file: chapter4/solvers
- file: chapter4/compiler_parameters
- file: chapter4/convergence
Expand Down
257 changes: 232 additions & 25 deletions chapter4/mixed_poisson.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,11 @@
"outputs": [],
"source": [
"from mpi4py import MPI\n",
"from petsc4py import PETSc\n",
"import dolfinx\n",
"\n",
"mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 50, 50)"
"N = 400\n",
"mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N)"
]
},
{
Expand Down Expand Up @@ -223,19 +225,20 @@
"metadata": {},
"source": [
"This can be split into a saddle point problem, with discretized matrices $A$ and $B$ and discretized\n",
"right-hand side $\\tilde f$.\n",
"right-hand side $\\mathbf{b}$.\n",
"\\begin{align}\n",
"\\begin{pmatrix}\n",
"A & B^T\\\\\n",
"B & 0\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
"u_h\\\\\n",
"sigma_h\n",
"\\sigma_h\n",
"\\end{pmatrix}\n",
"= \\begin{pmatrix}\n",
"0\\\\\n",
"\\tilde f\n",
"b_0\\\\\n",
"b_1\n",
"\\end{pmatrix}\n",
"\\end{align}\n",
"We can extract the block structure of the bilinear form using `ufl.extract_blocks`, which returns a nested list of bilinear forms.\n",
"You can also build this nested list by hand if you want to, but it is usually more error-prone."
Expand All @@ -261,13 +264,15 @@
{
"cell_type": "markdown",
"id": "3a66ff5c",
"metadata": {
"lines_to_next_cell": 2
},
"metadata": {},
"source": [
"Next we create the Dirichlet boundary condition for $\\sigma$.\n",
"For `sigma`, we have a manufactured solution that depends on the normal vector on the boundary,\n",
"which makes it slightly more complicated to implement."
"As we are using manufactured solutions for this problem, we could manually derive the explicit expression\n",
"for $\\sigma$ on the boundary $\\Gamma_N$.\n",
"However, in general this is not possible (especially for curved boundaries), and we have to use a more generic approach.\n",
"For this we will use the `dolfinx.fem.Expression` class to interpolate the expression into the function space $Q$.\n",
"This is done by evaluating the expression at the physical interpolation points of the mesh.\n",
"A convenience function for this is provided in the `interpolate_facet_expression` function below."
]
},
{
Expand Down Expand Up @@ -448,42 +453,244 @@
"metadata": {},
"outputs": [],
"source": [
"(sigma_h, u_h) = problem.solve()"
"import time"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5cea57f4",
"metadata": {
"lines_to_next_cell": 2
},
"metadata": {},
"outputs": [],
"source": [
"start = time.perf_counter()\n",
"(sigma_h, u_h) = problem.solve()\n",
"end = time.perf_counter()\n",
"print(f\"Direct solver took {end - start:.2f} seconds.\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5e61d522",
"metadata": {},
"outputs": [],
"source": [
"L2_u = dolfinx.fem.form(ufl.inner(u_h - u_exact, u_h - u_exact) * ufl.dx)\n",
"L2_sigma = dolfinx.fem.form(\n",
" ufl.inner(sigma_h - ufl.grad(u_exact), sigma_h - ufl.grad(u_exact)) * ufl.dx\n",
"Hdiv_sigma = dolfinx.fem.form(\n",
" ufl.inner(\n",
" ufl.div(sigma_h) - ufl.div(ufl.grad(u_exact)),\n",
" ufl.div(sigma_h) - ufl.div(ufl.grad(u_exact)),\n",
" )\n",
" * ufl.dx\n",
")\n",
"local_u_error = dolfinx.fem.assemble_scalar(L2_u)\n",
"local_sigma_error = dolfinx.fem.assemble_scalar(L2_sigma)\n",
"local_sigma_error = dolfinx.fem.assemble_scalar(Hdiv_sigma)\n",
"u_error = np.sqrt(mesh.comm.allreduce(local_u_error, op=MPI.SUM))\n",
"sigma_error = np.sqrt(mesh.comm.allreduce(local_sigma_error, op=MPI.SUM))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5e61d522",
"id": "94997937",
"metadata": {},
"outputs": [],
"source": [
"print(f\"Direct solver, L2(u): {u_error:.2e}, H(div)(sigma): {sigma_error:.2e}\")"
]
},
{
"cell_type": "markdown",
"id": "c9e08cf0",
"metadata": {},
"source": [
"## Iterative solver with Schur complement preconditioner\n",
"As mentioned earlier, there are more efficient ways of solving this problem, than using a direct solver.\n",
"Especially with the saddle point structure of the problem, we can use a Schur complement preconditioner.\n",
"As described in [FEniCSx PCTools: Mixed Poisson](https://rafinex-external-rifle.gitlab.io/fenicsx-pctools/demo/demo_mixed-poisson.html),\n",
"Instead of wrapping the matrices in a custom wrapper, we can use `dolfinx.fem.petsc.LinearProblem` to solve the problem."
]
},
{
"cell_type": "markdown",
"id": "a83cc5c8",
"metadata": {},
"source": [
"We start by defining the $S$ matrix in the Schur complement (see the aforementioned link for details on the variational formulation)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "edf158fb",
"metadata": {},
"outputs": [],
"source": [
"alpha = dolfinx.fem.Constant(mesh, 4.0)\n",
"gamma = dolfinx.fem.Constant(mesh, 9.0)\n",
"h = ufl.CellDiameter(mesh)\n",
"s = -(\n",
" ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx\n",
" - ufl.inner(ufl.avg(ufl.grad(v)), ufl.jump(u, n)) * ufl.dS\n",
" - ufl.inner(ufl.jump(u, n), ufl.avg(ufl.grad(v))) * ufl.dS\n",
" + (alpha / ufl.avg(h)) * ufl.inner(ufl.jump(u, n), ufl.jump(v, n)) * ufl.dS\n",
" - ufl.inner(ufl.grad(u), v * n) * dGammaD\n",
" - ufl.inner(u * n, ufl.grad(v)) * dGammaD\n",
" + (gamma / h) * u * v * dGammaD\n",
")\n",
"\n",
"S = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(s))\n",
"S.assemble()\n",
"\n",
"\n",
"class SchurInv:\n",
" def setUp(self, pc):\n",
" self.ksp = PETSc.KSP().create(mesh.comm)\n",
" self.ksp.setOptionsPrefix(pc.getOptionsPrefix() + \"SchurInv_\")\n",
" self.ksp.setOperators(S)\n",
" self.ksp.setTolerances(atol=1e-10, rtol=1e-10)\n",
" self.ksp.setFromOptions()\n",
"\n",
" def apply(self, pc, x, y):\n",
" self.ksp.solve(x, y)\n",
"\n",
" def __del__(self):\n",
" self.ksp.destroy()"
]
},
{
"cell_type": "markdown",
"id": "be27eab0",
"metadata": {},
"source": [
"Next we can create the linear problem instance with all the required options"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2dd5b523",
"metadata": {},
"outputs": [],
"source": [
"u_it = dolfinx.fem.Function(V, name=\"u_it\")\n",
"sigma_it = dolfinx.fem.Function(Q, name=\"sigma_it\")\n",
"petsc_options = {\n",
" \"ksp_error_if_not_converged\": True,\n",
" \"ksp_type\": \"gmres\",\n",
" \"ksp_rtol\": 1e-10,\n",
" \"ksp_atol\": 1e-10,\n",
" \"pc_type\": \"fieldsplit\",\n",
" \"pc_fieldsplit_type\": \"schur\",\n",
" \"pc_fieldsplit_schur_fact_type\": \"upper\",\n",
" \"pc_fieldsplit_schur_precondition\": \"user\",\n",
" f\"fieldsplit_{sigma_it.name}_0_ksp_type\": \"preonly\",\n",
" f\"fieldsplit_{sigma_it.name}_0_pc_type\": \"bjacobi\",\n",
" f\"fieldsplit_{u_it.name}_1_ksp_type\": \"preonly\",\n",
" f\"fieldsplit_{u_it.name}_1_pc_type\": \"python\",\n",
" f\"fieldsplit_{u_it.name}_1_pc_python_type\": __name__ + \".SchurInv\",\n",
" f\"fieldsplit_{u_it.name}_1_SchurInv_ksp_type\": \"preonly\",\n",
" f\"fieldsplit_{u_it.name}_1_SchurInv_pc_type\": \"hypre\",\n",
"}\n",
"w_it = (sigma_it, u_it)\n",
"problem = dolfinx.fem.petsc.LinearProblem(\n",
" a_blocked,\n",
" L_blocked,\n",
" u=w_it,\n",
" bcs=[bc_sigma],\n",
" petsc_options=petsc_options,\n",
" petsc_options_prefix=\"mp_\",\n",
" kind=\"nest\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "1a4c11f9",
"metadata": {},
"source": [
"```{admonition} NEST matrices\n",
"Note that instead of using `kind=\"mpi\"` we use `kind=\"nest\"` to indicate that we want to use a nested matrix structure\n",
"and employ the power of [PETSc fieldsplit](https://petsc.org/release/manual/ksp/#solving-block-matrices-with-pcfieldsplit).\n",
"```\n",
"The only modification required to linear problem is to attach the correct fieldsplit indexset."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7ef8ada6",
"metadata": {},
"outputs": [],
"source": [
"nest_IS = problem.A.getNestISs()\n",
"fieldsplit_IS = tuple(\n",
" [\n",
" (f\"{u.name + '_' if u.name != 'f' else ''}{i}\", IS)\n",
" for i, (u, IS) in enumerate(zip(problem.u, nest_IS[0]))\n",
" ]\n",
")\n",
"problem.solver.getPC().setFieldSplitIS(*fieldsplit_IS)\n",
"start_it = time.perf_counter()\n",
"problem.solve()\n",
"end_it = time.perf_counter()\n",
"print(\n",
" f\"Iterative solver took {end_it - start_it:.2f} seconds\"\n",
" + f\" in {problem.solver.getIterationNumber()} iterations\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "b79e329c",
"metadata": {},
"source": [
"We compute the error norms for the iterative solution"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2db0ee6f",
"metadata": {},
"outputs": [],
"source": [
"L2_u_it = dolfinx.fem.form(ufl.inner(u_it - u_exact, u_it - u_exact) * ufl.dx)\n",
"Hdiv_sigma_it = dolfinx.fem.form(\n",
" ufl.inner(\n",
" ufl.div(sigma_it) - ufl.div(ufl.grad(u_exact)),\n",
" ufl.div(sigma_it) - ufl.div(ufl.grad(u_exact)),\n",
" )\n",
" * ufl.dx\n",
")\n",
"local_u_error_it = dolfinx.fem.assemble_scalar(L2_u_it)\n",
"local_sigma_error_it = dolfinx.fem.assemble_scalar(Hdiv_sigma_it)\n",
"u_error_it = np.sqrt(mesh.comm.allreduce(local_u_error_it, op=MPI.SUM))\n",
"sigma_error_it = np.sqrt(mesh.comm.allreduce(local_sigma_error_it, op=MPI.SUM))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "478187c3",
"metadata": {},
"outputs": [],
"source": [
"print(f\"Iterative solver, L2(u): {u_error_it:.2e}, H(div)(sigma): {sigma_error_it:.2e}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "904e4b6a",
"metadata": {
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"print(f\"u error: {u_error:.3e}\")\n",
"print(f\"sigma error: {sigma_error:.3e}\")\n",
"u_D.name = \"u_ex\"\n",
"with dolfinx.io.XDMFFile(mesh.comm, \"mixed_poisson.xdmf\", \"w\") as xdmf:\n",
" xdmf.write_mesh(mesh)\n",
" xdmf.write_function(u_h, 0.0)\n",
" xdmf.write_function(u_D, 0.0)"
"np.testing.assert_allclose(u_h.x.array, u_it.x.array, rtol=1e-7, atol=1e-7)\n",
"np.testing.assert_allclose(sigma_h.x.array, sigma_it.x.array, rtol=1e-7, atol=1e-7)"
]
},
{
Expand Down
Loading