Update main branch to work with nightly.

Changelog.md

@@ -2,6 +2,8 @@
 
 ## v0.10.0
 
+- Full refactoring of `dolfinx.fem.petsc.NonlinearProblem`, which now uses the PETSc SNES backend. See [the non-linear poisson demo](./chapter2/nonlinpoisson_code.ipynb) for details.
+- `dolfinx.fem.petsc.LinearProblem` now requires an additional argument, `petsc_options_prefix`. This should be a unique string identifier for each `LinearProblem` that is created.
 - Change how one reads in GMSH data with `gmshio`. See [the membrane code](./chapter1/membrane_code.ipynb) for more details.
 - `dolfinx.fem.FiniteElement.interpolation_points()` -> `dolfinx.fem.FiniteElement.interpolation_points`.
 

README.md

@@ -21,7 +21,7 @@ Alternatively, if you want to add a separate chapter, a Jupyter notebook can be
 Any code added to the tutorial should work in parallel. If any changes are made to `ipynb` files, please ensure that these changes are reflected in the corresponding `py` files by using [`jupytext`](https://jupytext.readthedocs.io/en/latest/faq.html#can-i-use-jupytext-with-jupyterhub-binder-nteract-colab-saturn-or-azure):
 
 ```bash
-python3 -m jupytext --sync  */*.ipynb
+python3 -m jupytext --sync  */*.ipynb --set-formats ipynb,py:light
 ```
 
 Any code added to the tutorial should work in parallel.
@@ -61,4 +61,3 @@ from the root of this repository, and run
 ```
 
 from the main directory.
-

_config.yml

@@ -68,3 +68,4 @@ html:
     </div>
 
 exclude_patterns: [README.md, chapter2/advdiffreac.md]
+only_build_toc_files: true

chapter1/complex_mode.ipynb

@@ -11,7 +11,7 @@
     "\n",
     "Many PDEs, such as the [Helmholtz equation](https://docs.fenicsproject.org/dolfinx/v0.4.1/python/demos/demo_helmholtz.html) require complex-valued fields.\n",
     "\n",
-    "For simplicity, let us consider a Poisson equation of the form: \n",
+    "For simplicity, let us consider a Poisson equation of the form:\n",
     "\n",
     "$$-\\Delta u = f \\text{ in } \\Omega,$$\n",
     "$$ f = -1 - 2j \\text{ in } \\Omega,$$\n",
@@ -47,12 +47,13 @@
     "from mpi4py import MPI\n",
     "import dolfinx\n",
     "import numpy as np\n",
+    "\n",
     "mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)\n",
     "V = dolfinx.fem.functionspace(mesh, (\"Lagrange\", 1))\n",
-    "u_r = dolfinx.fem.Function(V, dtype=np.float64) \n",
+    "u_r = dolfinx.fem.Function(V, dtype=np.float64)\n",
     "u_r.interpolate(lambda x: x[0])\n",
     "u_c = dolfinx.fem.Function(V, dtype=np.complex128)\n",
-    "u_c.interpolate(lambda x:0.5*x[0]**2 + 1j*x[1]**2)\n",
+    "u_c.interpolate(lambda x: 0.5 * x[0] ** 2 + 1j * x[1] ** 2)\n",
     "print(u_r.x.array.dtype)\n",
     "print(u_c.x.array.dtype)"
    ]
@@ -78,8 +79,9 @@
    "source": [
     "from petsc4py import PETSc\n",
     "from dolfinx.fem.petsc import assemble_vector\n",
+    "\n",
     "print(PETSc.ScalarType)\n",
-    "assert np.dtype(PETSc.ScalarType).kind == 'c'"
+    "assert np.dtype(PETSc.ScalarType).kind == \"c\""
    ]
   },
   {
@@ -99,6 +101,7 @@
    "outputs": [],
    "source": [
     "import ufl\n",
+    "\n",
     "u = ufl.TrialFunction(V)\n",
     "v = ufl.TestFunction(V)\n",
     "f = dolfinx.fem.Constant(mesh, PETSc.ScalarType(-1 - 2j))\n",
@@ -167,11 +170,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)\n",
+    "mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)\n",
     "boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)\n",
-    "boundary_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, boundary_facets)\n",
+    "boundary_dofs = dolfinx.fem.locate_dofs_topological(\n",
+    "    V, mesh.topology.dim - 1, boundary_facets\n",
+    ")\n",
     "bc = dolfinx.fem.dirichletbc(u_c, boundary_dofs)\n",
-    "problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc])\n",
+    "problem = dolfinx.fem.petsc.LinearProblem(\n",
+    "    a, L, bcs=[bc], petsc_options_prefix=\"complex_poisson\"\n",
+    ")\n",
     "uh = problem.solve()"
    ]
   },
@@ -193,11 +200,13 @@
    "outputs": [],
    "source": [
     "x = ufl.SpatialCoordinate(mesh)\n",
-    "u_ex = 0.5 * x[0]**2 + 1j*x[1]**2\n",
-    "L2_error = dolfinx.fem.form(ufl.dot(uh-u_ex, uh-u_ex) * ufl.dx(metadata={\"quadrature_degree\": 5}))\n",
+    "u_ex = 0.5 * x[0] ** 2 + 1j * x[1] ** 2\n",
+    "L2_error = dolfinx.fem.form(\n",
+    "    ufl.dot(uh - u_ex, uh - u_ex) * ufl.dx(metadata={\"quadrature_degree\": 5})\n",
+    ")\n",
     "local_error = dolfinx.fem.assemble_scalar(L2_error)\n",
     "global_error = np.sqrt(mesh.comm.allreduce(local_error, op=MPI.SUM))\n",
-    "max_error = mesh.comm.allreduce(np.max(np.abs(u_c.x.array-uh.x.array)))\n",
+    "max_error = mesh.comm.allreduce(np.max(np.abs(u_c.x.array - uh.x.array)))\n",
     "print(global_error, max_error)"
    ]
   },
@@ -219,38 +228,23 @@
    "outputs": [],
    "source": [
     "import pyvista\n",
-    "pyvista.start_xvfb()\n",
+    "\n",
+    "pyvista.start_xvfb(0.1)\n",
     "mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)\n",
     "p_mesh = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh, mesh.topology.dim))\n",
     "pyvista_cells, cell_types, geometry = dolfinx.plot.vtk_mesh(V)\n",
     "grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry)\n",
     "grid.point_data[\"u_real\"] = uh.x.array.real\n",
     "grid.point_data[\"u_imag\"] = uh.x.array.imag\n",
-    "_ = grid.set_active_scalars(\"u_real\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "16",
-   "metadata": {},
-   "outputs": [],
-   "source": [
+    "_ = grid.set_active_scalars(\"u_real\")\n",
+    "\n",
     "p_real = pyvista.Plotter()\n",
     "p_real.add_text(\"uh real\", position=\"upper_edge\", font_size=14, color=\"black\")\n",
     "p_real.add_mesh(grid, show_edges=True)\n",
     "p_real.view_xy()\n",
     "if not pyvista.OFF_SCREEN:\n",
-    "    p_real.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "17",
-   "metadata": {},
-   "outputs": [],
-   "source": [
+    "    p_real.show()\n",
+    "\n",
     "grid.set_active_scalars(\"u_imag\")\n",
     "p_imag = pyvista.Plotter()\n",
     "p_imag.add_text(\"uh imag\", position=\"upper_edge\", font_size=14, color=\"black\")\n",
@@ -259,14 +253,6 @@
     "if not pyvista.OFF_SCREEN:\n",
     "    p_imag.show()"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "18",
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {

chapter1/complex_mode.py

@@ -6,7 +6,7 @@
 #       extension: .py
 #       format_name: light
 #       format_version: '1.5'
-#       jupytext_version: 1.16.5
+#       jupytext_version: 1.17.2
 #   kernelspec:
 #     display_name: Python 3 (DOLFINx complex)
 #     language: python
@@ -19,7 +19,7 @@
 #
 # Many PDEs, such as the [Helmholtz equation](https://docs.fenicsproject.org/dolfinx/v0.4.1/python/demos/demo_helmholtz.html) require complex-valued fields.
 #
-# For simplicity, let us consider a Poisson equation of the form: 
+# For simplicity, let us consider a Poisson equation of the form:
 #
 # $$-\Delta u = f \text{ in } \Omega,$$
 # $$ f = -1 - 2j \text{ in } \Omega,$$
@@ -45,38 +45,47 @@
 # FEniCSx supports both real and complex numbers, so we can create a function space with real valued or complex valued coefficients.
 #
 
+# +
 from mpi4py import MPI
 import dolfinx
 import numpy as np
+
 mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
 V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
-u_r = dolfinx.fem.Function(V, dtype=np.float64) 
+u_r = dolfinx.fem.Function(V, dtype=np.float64)
 u_r.interpolate(lambda x: x[0])
 u_c = dolfinx.fem.Function(V, dtype=np.complex128)
-u_c.interpolate(lambda x:0.5*x[0]**2 + 1j*x[1]**2)
+u_c.interpolate(lambda x: 0.5 * x[0] ** 2 + 1j * x[1] ** 2)
 print(u_r.x.array.dtype)
 print(u_c.x.array.dtype)
+# -
 
 # However, as we would like to solve linear algebra problems of the form $Ax=b$, we need to be able to use matrices and vectors that support real and complex numbers. As [PETSc](https://petsc.org/release/) is one of the most popular interfaces to linear algebra packages, we need to be able to work with their matrix and vector structures.
 #
 # Unfortunately, PETSc only supports one floating type in their matrices, thus we need to install two versions of PETSc, one that supports `float64` and one that supports `complex128`. In the [docker images](https://hub.docker.com/r/dolfinx/dolfinx) for DOLFINx, both versions are installed, and one can switch between them by calling `source dolfinx-real-mode` or `source dolfinx-complex-mode`. For the `dolfinx/lab` images, one can change the Python kernel to be either the real or complex mode, by going to `Kernel->Change Kernel...` and choosing `Python3 (ipykernel)` (for real mode) or `Python3 (DOLFINx complex)` (for complex mode).
 #
 # We check that we are using the correct installation of PETSc by inspecting the scalar type.
 
+# +
 from petsc4py import PETSc
 from dolfinx.fem.petsc import assemble_vector
+
 print(PETSc.ScalarType)
-assert np.dtype(PETSc.ScalarType).kind == 'c'
+assert np.dtype(PETSc.ScalarType).kind == "c"
+# -
 
 # ## Variational problem
 # We are now ready to define our variational problem
 
+# +
 import ufl
+
 u = ufl.TrialFunction(V)
 v = ufl.TestFunction(V)
 f = dolfinx.fem.Constant(mesh, PETSc.ScalarType(-1 - 2j))
 a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
 L = ufl.inner(f, v) * ufl.dx
+# -
 
 # Note that we have used the `PETSc.ScalarType` to wrap the constant source on the right hand side. This is because we want the integration kernels to assemble into the correct floating type.
 #
@@ -99,11 +108,15 @@
 # We define our Dirichlet condition and setup and solve the variational problem.
 # ## Solve variational problem
 
-mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
+mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
 boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
-boundary_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, boundary_facets)
+boundary_dofs = dolfinx.fem.locate_dofs_topological(
+    V, mesh.topology.dim - 1, boundary_facets
+)
 bc = dolfinx.fem.dirichletbc(u_c, boundary_dofs)
-problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc])
+problem = dolfinx.fem.petsc.LinearProblem(
+    a, L, bcs=[bc], petsc_options_prefix="complex_poisson"
+)
 uh = problem.solve()
 
 # We compute the $L^2$ error and the max error.
@@ -112,20 +125,24 @@
 #
 
 x = ufl.SpatialCoordinate(mesh)
-u_ex = 0.5 * x[0]**2 + 1j*x[1]**2
-L2_error = dolfinx.fem.form(ufl.dot(uh-u_ex, uh-u_ex) * ufl.dx(metadata={"quadrature_degree": 5}))
+u_ex = 0.5 * x[0] ** 2 + 1j * x[1] ** 2
+L2_error = dolfinx.fem.form(
+    ufl.dot(uh - u_ex, uh - u_ex) * ufl.dx(metadata={"quadrature_degree": 5})
+)
 local_error = dolfinx.fem.assemble_scalar(L2_error)
 global_error = np.sqrt(mesh.comm.allreduce(local_error, op=MPI.SUM))
-max_error = mesh.comm.allreduce(np.max(np.abs(u_c.x.array-uh.x.array)))
+max_error = mesh.comm.allreduce(np.max(np.abs(u_c.x.array - uh.x.array)))
 print(global_error, max_error)
 
 # ## Plotting
 #
 # Finally, we plot the real and imaginary solutions.
 #
 
+# +
 import pyvista
-pyvista.start_xvfb()
+
+pyvista.start_xvfb(0.1)
 mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
 p_mesh = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh, mesh.topology.dim))
 pyvista_cells, cell_types, geometry = dolfinx.plot.vtk_mesh(V)
@@ -148,5 +165,3 @@
 p_imag.view_xy()
 if not pyvista.OFF_SCREEN:
     p_imag.show()
-
-