from petsc4py import PETSc
from mpi4py import MPI
import dolfinx as df
import ufl
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, set_bc
def project(e, target_func, bcs=[]):
"""Project UFL expression.
Note
----
This method solves a linear system (using KSP defaults).
"""
# Ensure we have a mesh and attach to measure
V = target_func.function_space
dx = ufl.dx(V.mesh)
# Define variational problem for projection
w = ufl.TestFunction(V)
v = ufl.TrialFunction(V)
a = df.fem.form(ufl.inner(v, w) * dx)
L = df.fem.form(ufl.inner(e, w) * dx)
# Assemble linear system
A = assemble_matrix(a, bcs)
A.assemble()
b = assemble_vector(L)
apply_lifting(b, [a], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)
# Solve linear system
solver = PETSc.KSP().create(A.getComm())
solver.setType("bcgs")
solver.getPC().setType("bjacobi")
solver.rtol = 1.0e-05
solver.setOperators(A)
solver.solve(b, target_func.x.petsc_vec) # ERROR: Changed from target_func.vector
assert solver.reason > 0
target_func.x.scatter_forward()
# Destroy PETSc linear algebra objects and solver
solver.destroy()
A.destroy()
b.destroy()
mesh = df.mesh.create_unit_square(MPI.COMM_WORLD,8,8) #Mesh is whatever you want it to be, in my use case it was a sphere of radius 1
V = df.fem.functionspace(mesh, ("DG", 1))
x = ufl.SpatialCoordinate(mesh)
f = ufl.exp(-((-x[1]-1)/0.1)**2)
D0 = df.fem.Function(V)
project(f,D0)
Does this look correct for projecting function f on function space V? Also, is there maybe a simpler way to do this?
I would suggest making the projector reusable, as explained in:
which has the code in dropdown:
from typing import Optional
import numpy as np
import pyvista
class Projector:
"""
Projector for a given function.
Solves Ax=b, where
.. highlight:: python
.. code-block:: python
u, v = ufl.TrialFunction(Space), ufl.TestFunction(space)
dx = ufl.Measure("dx", metadata=metadata)
A = inner(u, v) * dx
b = inner(function, v) * dx(metadata=metadata)
Args:
function: UFL expression of function to project
space: Space to project function into
petsc_options: Options to pass to PETSc
jit_options: Options to pass to just in time compiler
form_compiler_options: Options to pass to the form compiler
metadata: Data to pass to the integration measure
"""
_A: PETSc.Mat # The mass matrix
_b: PETSc.Vec # The rhs vector
_lhs: dolfinx.fem.Form # The compiled form for the mass matrix
_ksp: PETSc.KSP # The PETSc solver
_x: dolfinx.fem.Function # The solution vector
_dx: ufl.Measure # Integration measure
def __init__(
self,
space: dolfinx.fem.FunctionSpace,
petsc_options: Optional[dict] = None,
jit_options: Optional[dict] = None,
form_compiler_options: Optional[dict] = None,
metadata: Optional[dict] = None,
):
petsc_options = {} if petsc_options is None else petsc_options
jit_options = {} if jit_options is None else jit_options
form_compiler_options = {} if form_compiler_options is None else form_compiler_options
# Assemble projection matrix once
u = ufl.TrialFunction(space)
v = ufl.TestFunction(space)
self._dx = ufl.Measure("dx", domain=space.mesh, metadata=metadata)
a = ufl.inner(u, v) * self._dx(metadata=metadata)
self._lhs = dolfinx.fem.form(a, jit_options=jit_options, form_compiler_options=form_compiler_options)
self._A = dolfinx.fem.petsc.assemble_matrix(self._lhs)
self._A.assemble()
# Create vectors to store right hand side and the solution
self._x = dolfinx.fem.Function(space)
self._b = dolfinx.fem.Function(space)
# Create Krylov Subspace solver
self._ksp = PETSc.KSP().create(space.mesh.comm)
self._ksp.setOperators(self._A)
# Set PETSc options
prefix = f"projector_{id(self)}"
opts = PETSc.Options()
opts.prefixPush(prefix)
for k, v in petsc_options.items():
opts[k] = v
opts.prefixPop()
self._ksp.setFromOptions()
for opt in opts.getAll().keys():
del opts[opt]
# Set matrix and vector PETSc options
self._A.setOptionsPrefix(prefix)
self._A.setFromOptions()
self._b.x.petsc_vec.setOptionsPrefix(prefix)
self._b.x.petsc_vec.setFromOptions()
def reassemble_lhs(self):
dolfinx.fem.petsc.assemble_matrix(self._A, self._lhs)
self._A.assemble()
def assemble_rhs(self, h: ufl.core.expr.Expr):
"""
Assemble the right hand side of the problem
"""
v = ufl.TestFunction(self._b.function_space)
rhs = ufl.inner(h, v) * self._dx
rhs_compiled = dolfinx.fem.form(rhs)
self._b.x.array[:] = 0.0
dolfinx.fem.petsc.assemble_vector(self._b.x.petsc_vec, rhs_compiled)
self._b.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
self._b.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
def project(self, h: ufl.core.expr.Expr) -> dolfinx.fem.Function:
"""
Compute projection using a PETSc KSP solver
Args:
assemble_rhs: Re-assemble RHS and re-apply boundary conditions if true
"""
self.assemble_rhs(h)
self._ksp.solve(self._b.x.petsc_vec, self._x.x.petsc_vec)
return self._x
def __del__(self):
self._A.destroy()
self._ksp.destroy()
Thank you! That answers both questions. I’ve created a separate file .py file with the class inside of it to call when needed so that it can be instantiated from other files for future use. The code you’ve provide works just as well without issue and provides the same results, so it looks like what I had was correct.