PETSc python API update

Changelog.md

@@ -1,7 +1,9 @@
 # Changelog
 
 ## main
-- No changes
+- `dolfinx.geometry.BoundingBoxTree` has been changed to `dolfinx.geometry.bb_tree`
+- `create_mesh` with Meshio has been modified. Note that you now need to pass dtype `np.int32` to the cell_data.
+- Update dolfinx petsc API. Now one needs to explicitly import `dolfinx.fem.petsc` and `dolfinx.fem.nls`, as PETSc is no longer a strict requirement. Replace `petsc4py.PETSc.ScalarType` with `dolfinx.default_scalar_type` in demos where we do not use `petsc4py` explicitly.
 
 ## v0.6.0
 - Remove `ipygany` and `pythreejs` as plotting backends. Using `panel`.

chapter1/complex_mode.py

@@ -62,6 +62,7 @@
 # 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'
 
@@ -89,7 +90,7 @@
 
 J = u_c**2 * ufl.dx
 F = ufl.derivative(J, u_c, ufl.conj(v))
-residual = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(F))
+residual = assemble_vector(dolfinx.fem.form(F))
 print(residual.array)
 
 # We define our Dirichlet condition and setup and solve the variational problem.

chapter1/fundamentals.md

@@ -120,7 +120,7 @@ For the Poisson equation, we have:
 a(u,v) &= \int_{\Omega} \nabla u \cdot \nabla v ~\mathrm{d} x,\\
 L(v) &= \int_{\Omega} fv ~\mathrm{d} x.
 \end{align}
-From literature $a(u,v)$ is known as the _bilinear form_ and $L(V)$ as a _linear form_. 
+From literature $a(u,v)$ is known as the _bilinear form_ and $L(v)$ as a _linear form_. 
 For every linear problem, we will identify all terms with the unknown $u$ and collect them in $a(u,v)$, and collect all terms with only  known functions in $L(v)$.
 
 To solve a linear PDE in FEniCSx, such as the Poisson equation, a user thus needs to perform two steps:

chapter1/fundamentals_code.py

@@ -89,24 +89,6 @@
 # + vscode={"languageId": "python"}
 from dolfinx.fem import FunctionSpace
 V = FunctionSpace(domain, ("Lagrange", 1))
-# -
-
-# The second argument is the tuple containing the type of finite element, and the element degree. The type of element here is "Lagrange", which implies the standard Lagrange family of elements. 
-# DOLFINx supports a large variety on elements on simplices 
-# (triangles and tetrahedra) and non-simplices (quadrilaterals
-# and hexahedra). For an overview, see:
-# *FIXME: Add link to all the elements we support*
-#
-# The element degree in the code is 1. This means that we are choosing the standard $P_1$ linear Lagrange element, which has degrees of freedom at the vertices. 
-# The computed solution will be continuous across elements and linearly varying in $x$ and $y$ inside each element. Higher degree polynomial approximations are obtained by increasing the degree argument. 
-#
-# ## Defining the boundary conditions
-#
-# The next step is to specify the boundary condition $u=u_D$ on $\partial\Omega_D$, which is done by over several steps. 
-# The first step is to define the function $u_D$. Into this function, we would like to interpolate the boundary condition $1 + x^2+2y^2$.
-# We do this by first defining a `dolfinx.fem.Function`, and then using a [lambda-function](https://docs.python.org/3/tutorial/controlflow.html#lambda-expressions) in Python to define the 
-# spatially varying function.
-#
 
 # + vscode={"languageId": "python"}
 from dolfinx import fem
@@ -158,8 +140,8 @@
 # As the source term is constant over the domain, we use `dolfinx.Constant`
 
 # + vscode={"languageId": "python"}
-from petsc4py.PETSc import ScalarType
-f = fem.Constant(domain, ScalarType(-6))
+from dolfinx import default_scalar_type
+f = fem.Constant(domain, default_scalar_type(-6))
 # -
 
 # ```{admonition} Compilation speed-up
@@ -194,10 +176,13 @@
 # ## Forming and solving the linear system
 #
 # Having defined the finite element variational problem and boundary condition, we can create our `dolfinx.fem.petsc.LinearProblem`, as class for solving 
-# the variational problem: Find $u_h\in V$ such that $a(u_h, v)==L(v) \quad \forall v \in \hat{V}$. We will use PETSc as our linear algebra backend, using a direct solver (LU-factorization).  See the [PETSc-documentation](https://petsc.org/main/docs/manual/ksp/?highlight=ksp#ksp-linear-system-solvers) of the method for more information.
+# the variational problem: Find $u_h\in V$ such that $a(u_h, v)==L(v) \quad \forall v \in \hat{V}$. We will use PETSc as our linear algebra backend, using a direct solver (LU-factorization).
+# See the [PETSc-documentation](https://petsc.org/main/docs/manual/ksp/?highlight=ksp#ksp-linear-system-solvers) of the method for more information.
+# PETSc is not a required dependency of DOLFINx, and therefore we explicitly import the DOLFINx wrapper for interfacing with PETSc.
 
 # + vscode={"languageId": "python"}
-problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
+from dolfinx.fem.petsc import LinearProblem
+problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
 uh = problem.solve()
 # -
 

chapter1/membrane_code.py

@@ -54,6 +54,7 @@
 
 # +
 from dolfinx.io import gmshio
+from dolfinx.fem.petsc import LinearProblem
 from mpi4py import MPI
 
 gmsh_model_rank = 0
@@ -70,10 +71,10 @@
 # The right hand side pressure function is represented using `ufl.SpatialCoordinate` and two constants, one for $\beta$ and one for $R_0$.
 
 import ufl
-from petsc4py.PETSc import ScalarType
+from dolfinx import default_scalar_type
 x = ufl.SpatialCoordinate(domain)
-beta = fem.Constant(domain, ScalarType(12))
-R0 = fem.Constant(domain, ScalarType(0.3))
+beta = fem.Constant(domain, default_scalar_type(12))
+R0 = fem.Constant(domain, default_scalar_type(0.3))
 p = 4 * ufl.exp(-beta**2 * (x[0]**2 + (x[1] - R0)**2))
 
 # ## Create a Dirichlet boundary condition using geometrical conditions
@@ -92,7 +93,7 @@ def on_boundary(x):
 
 # As our Dirichlet condition is homogeneous (`u=0` on the whole boundary), we can initialize the `dolfinx.fem.dirichletbc` with a constant value, the degrees of freedom and the function space to apply the boundary condition on.
 
-bc = fem.dirichletbc(ScalarType(0), boundary_dofs, V)
+bc = fem.dirichletbc(default_scalar_type(0), boundary_dofs, V)
 
 # ## Defining the variational problem
 # The variational problem is the same as in our first Poisson problem, where `f` is replaced by `p`.
@@ -101,7 +102,7 @@ def on_boundary(x):
 v = ufl.TestFunction(V)
 a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
 L = p * v * ufl.dx
-problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
+problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
 uh = problem.solve()
 
 # ## Interpolation of a `ufl`-expression

chapter1/nitsche.py

@@ -22,10 +22,10 @@
 # We start by importing the required modules and creating the mesh and function space for our solution
 
 # +
-from dolfinx import fem, mesh, plot
+from dolfinx import fem, mesh, plot, default_scalar_type
+from dolfinx.fem.petsc import LinearProblem
 import numpy
 from mpi4py import MPI
-from petsc4py.PETSc import ScalarType 
 from ufl import (Circumradius, FacetNormal, SpatialCoordinate, TrialFunction, TestFunction,
                  div, dx, ds, grad, inner)
 
@@ -64,15 +64,15 @@
 v = TestFunction(V)
 n = FacetNormal(domain)
 h = 2 * Circumradius(domain)
-alpha = fem.Constant(domain, ScalarType(10))
+alpha = fem.Constant(domain, default_scalar_type(10))
 a = inner(grad(u), grad(v)) * dx - inner(n, grad(u)) * v * ds
 a += - inner(n, grad(v)) * u * ds + alpha / h * inner(u, v) * ds
 L = inner(f, v) * dx 
 L += - inner(n, grad(v)) * uD * ds + alpha / h * inner(uD, v) * ds
 
 # As we now have the variational form, we can solve the linear problem
 
-problem = fem.petsc.LinearProblem(a, L)
+problem = LinearProblem(a, L)
 uh = problem.solve()
 
 # We compute the error of the computation by comparing it to the analytical solution

chapter2/diffusion_code.ipynb

@@ -28,10 +28,11 @@
     "import ufl\n",
     "import numpy as np\n",
     "\n",
-    "from mpi4py import MPI\n",
     "from petsc4py import PETSc\n",
+    "from mpi4py import MPI\n",
     "\n",
     "from dolfinx import fem, mesh, io, plot\n",
+    "from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc\n",
     "\n",
     "# Define temporal parameters\n",
     "t = 0  # Start time\n",
@@ -153,9 +154,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "A = fem.petsc.assemble_matrix(bilinear_form, bcs=[bc])\n",
+    "A = assemble_matrix(bilinear_form, bcs=[bc])\n",
     "A.assemble()\n",
-    "b = fem.petsc.create_vector(linear_form)"
+    "b = create_vector(linear_form)"
    ]
   },
   {
@@ -241,12 +242,12 @@
     "    # Update the right hand side reusing the initial vector\n",
     "    with b.localForm() as loc_b:\n",
     "        loc_b.set(0)\n",
-    "    fem.petsc.assemble_vector(b, linear_form)\n",
+    "    assemble_vector(b, linear_form)\n",
     "\n",
     "    # Apply Dirichlet boundary condition to the vector\n",
-    "    fem.petsc.apply_lifting(b, [bilinear_form], [[bc]])\n",
+    "    apply_lifting(b, [bilinear_form], [[bc]])\n",
     "    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)\n",
-    "    fem.petsc.set_bc(b, [bc])\n",
+    "    set_bc(b, [bc])\n",
     "\n",
     "    # Solve linear problem\n",
     "    solver.solve(b, uh.vector)\n",

chapter2/diffusion_code.py

@@ -31,10 +31,11 @@
 import ufl
 import numpy as np
 
-from mpi4py import MPI
 from petsc4py import PETSc
+from mpi4py import MPI
 
 from dolfinx import fem, mesh, io, plot
+from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
 
 # Define temporal parameters
 t = 0  # Start time
@@ -103,9 +104,9 @@ def initial_condition(x, a=5):
 
 # We observe that the left hand side of the system, the matrix $A$ does not change from one time step to another, thus we only need to assemble it once. However, the right hand side, which is dependent on the previous time step `u_n`, we have to assemble it every time step. Therefore, we only create a vector `b` based on `L`, which we will reuse at every time step.
 
-A = fem.petsc.assemble_matrix(bilinear_form, bcs=[bc])
+A = assemble_matrix(bilinear_form, bcs=[bc])
 A.assemble()
-b = fem.petsc.create_vector(linear_form)
+b = create_vector(linear_form)
 
 # ## Using petsc4py to create a linear solver
 # As we have already assembled `a` into the matrix `A`, we can no longer use the `dolfinx.fem.petsc.LinearProblem` class to solve the problem. Therefore, we create a linear algebra solver using PETSc, and assign the matrix `A` to the solver, and choose the solution strategy.
@@ -153,12 +154,12 @@ def initial_condition(x, a=5):
     # Update the right hand side reusing the initial vector
     with b.localForm() as loc_b:
         loc_b.set(0)
-    fem.petsc.assemble_vector(b, linear_form)
+    assemble_vector(b, linear_form)
 
     # Apply Dirichlet boundary condition to the vector
-    fem.petsc.apply_lifting(b, [bilinear_form], [[bc]])
+    apply_lifting(b, [bilinear_form], [[bc]])
     b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
-    fem.petsc.set_bc(b, [bc])
+    set_bc(b, [bc])
 
     # Solve linear problem
     solver.solve(b, uh.vector)

chapter2/heat_code.ipynb

@@ -29,6 +29,7 @@
     "from mpi4py import MPI\n",
     "import ufl\n",
     "from dolfinx import mesh, fem\n",
+    "from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc\n",
     "import numpy\n",
     "t = 0  # Start time\n",
     "T = 2  # End time\n",
@@ -180,9 +181,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "A = fem.petsc.assemble_matrix(a, bcs=[bc])\n",
+    "A = assemble_matrix(a, bcs=[bc])\n",
     "A.assemble()\n",
-    "b = fem.petsc.create_vector(L)\n",
+    "b = create_vector(L)\n",
     "uh = fem.Function(V)"
    ]
   },
@@ -233,12 +234,12 @@
     "    # Update the right hand side reusing the initial vector\n",
     "    with b.localForm() as loc_b:\n",
     "        loc_b.set(0)\n",
-    "    fem.petsc.assemble_vector(b, L)\n",
+    "    assemble_vector(b, L)\n",
     "\n",
     "    # Apply Dirichlet boundary condition to the vector\n",
-    "    fem.petsc.apply_lifting(b, [a], [[bc]])\n",
+    "    apply_lifting(b, [a], [[bc]])\n",
     "    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)\n",
-    "    fem.petsc.set_bc(b, [bc])\n",
+    "    set_bc(b, [bc])\n",
     "\n",
     "    # Solve linear problem\n",
     "    solver.solve(b, uh.vector)\n",

chapter2/heat_code.py

@@ -31,6 +31,7 @@
 from mpi4py import MPI
 import ufl
 from dolfinx import mesh, fem
+from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
 import numpy
 t = 0  # Start time
 T = 2  # End time
@@ -96,9 +97,9 @@ def __call__(self, x):
 # ## Create the matrix and vector for the linear problem
 # To ensure that we are solving the variational problem efficiently, we will create several structures which can reuse data, such as matrix sparisty patterns. Especially note as the bilinear form `a` is independent of time, we only need to assemble the matrix once.
 
-A = fem.petsc.assemble_matrix(a, bcs=[bc])
+A = assemble_matrix(a, bcs=[bc])
 A.assemble()
-b = fem.petsc.create_vector(L)
+b = create_vector(L)
 uh = fem.Function(V)
 
 # ## Define a linear variational solver
@@ -125,12 +126,12 @@ def __call__(self, x):
     # Update the right hand side reusing the initial vector
     with b.localForm() as loc_b:
         loc_b.set(0)
-    fem.petsc.assemble_vector(b, L)
+    assemble_vector(b, L)
 
     # Apply Dirichlet boundary condition to the vector
-    fem.petsc.apply_lifting(b, [a], [[bc]])
+    apply_lifting(b, [a], [[bc]])
     b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
-    fem.petsc.set_bc(b, [bc])
+    set_bc(b, [bc])
 
     # Solve linear problem
     solver.solve(b, uh.vector)