I want to solve a Petrov-Galerkin Variational formulation using Fenicsx, i.e. a formulation where trial and test space are not equal, typically also including different boundary conditions. Small example:
I implemented that with f(x) = pi^2 sin(pi*x), b = 0, c = 1, i.e. the exact solution reads u(x) = sin(pi*x) with the following code:
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem
from dolfinx.fem.petsc import LinearProblem
import ufl
np.set_printoptions(edgeitems=30, linewidth=100000, precision=2)
# define mesh
msh = mesh.create_unit_interval(MPI.COMM_WORLD, 5)
msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
tdim = msh.topology.dim-1
# define problem
x = ufl.SpatialCoordinate(msh)
n = ufl.FacetNormal(msh)
u_exact = lambda x: np.sin(np.pi * x)
b = u_exact(0)
c = np.pi * np.cos(np.pi*0)
f = ufl.pi**2 * ufl.sin(ufl.pi * x[0])
# get the two boundary parts {0} and {1}
gamma0_facets = mesh.locate_entities_boundary(msh, dim=tdim, marker=lambda tx: np.isclose(tx[0], 0.0)) # Gamma_0 = {0}
gamma1_facets = mesh.locate_entities_boundary(msh, dim=tdim, marker=lambda tx: np.isclose(tx[0], 1.0)) # Gamma_1 = {1}
# construct the boundary measure over gamma_0
gamma0 = mesh.meshtags(msh, tdim, gamma0_facets, np.full_like(gamma0_facets, 0))
ds = ufl.Measure("ds", domain=msh, subdomain_data=gamma0)
ds_gamma0 = ds(0)
# define the two (!) function spaces -> Petrov-Galerkin
U = fem.functionspace(msh, ("Lagrange", 1))
V = fem.functionspace(msh, ("Lagrange", 1))
# setup the different (!) Dirichlet BCs for each space
gamma0_dofs_U = fem.locate_dofs_topological(U, tdim, gamma0_facets)
gamma1_dofs_V = fem.locate_dofs_topological(V, tdim, gamma1_facets)
bcs = [fem.dirichletbc(b, gamma0_dofs_U, U),
fem.dirichletbc(0.0, gamma1_dofs_V, V)]
# define variational problem
u = ufl.TrialFunction(U)
v = ufl.TestFunction(V)
A = ufl.dot(ufl.grad(u),ufl.grad(v)) * ufl.dx
l = f*v*ufl.dx - c*v*ds_gamma0
# print discrete system
print("A:\n" + str(fem.assemble_matrix(fem.form(A), bcs=bcs).to_scipy().todense()))
print("l:\n" + str(fem.assemble_vector(fem.form(l)).array))
# try to solve the problem -> failure
problem = LinearProblem(A,l, bcs=bcs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "ksp_error_if_not_converged" : True},
petsc_options_prefix="petrov_galerkin_"
)
uh = problem.solve()
But Fenics is not able to solve that problem, so the solving throws an error, and when we look at the linear system, it is clear why:
A:
[[ 0. -5. 0. 0. 0. 0.]
[ 0. 10. -5. 0. 0. 0.]
[ 0. -5. 10. -5. 0. 0.]
[ 0. 0. -5. 10. -5. 0.]
[ 0. 0. 0. -5. 10. -5.]
[ 0. 0. 0. 0. 0. 0.]]
l:
[-2.94 1.12 1.82 1.82 1.12 0.2 ]
The boundary conditions are not handled properly, resulting an all zero row and column. If i am not mistaking, Fenicsx handles Dirichlet BCs by 1) setting the corresponding rows and columns to zero and 2) placing a 1.0 on the diagonal, instead of 1) deleting these rows and columns completely and 2) modifying the right hand side. Of course, using this approach, the solver can not be able to solve a Petrov-Galerkin type problem, as it can no longer place a 1.0 at the diagonal (the fundamental assumption for that is that trial and test space dofs have a 1one1 relation, which is no longer guaranteed in the Petrov-Galerkin setting and is in fact not the case for the above example). Fenicsx / PETSc seems to realize that there is a problem, as it doesn’t place a 1.0 at the diagonal, but on the same time still sets the rows and columns to zero, thus leaving the system in an unfinished state (where the BCs are not handled properly and the resulting matrix is singular)
I can of course fix this on my own using numpy with the following code
A = fem.assemble_matrix(fem.form(A)).to_scipy().todense()
l = fem.assemble_vector(fem.form(l)).array
# identify free dofs (by removing dirichlet dofs)
dofs_U = np.ones(A.shape[1], dtype=bool)
dofs_U[dirichlet_dofs_U] = False
dofs_V = np.ones(A.shape[0], dtype=bool)
dofs_V[dirichlet_dofs_V] = False
# delete rows and columns corresponding to dirichlet dofs and modify right hand side
A_ = A[dofs_V,:][:,dofs_U]
l_ = l[dofs_V] + A[dofs_V, dirichlet_dofs_U] * b
# solve the reduced system and insert dirichlet values into the solution
u = np.zeros((A.shape[1],1))
u[dofs_U] = np.linalg.solve(A_,l_.T)
u[dirichlet_dofs_U] = b
import matplotlib.pyplot as plt
plt.plot(u)
plt.show()
This then gives a perfect sin function and hence reproduces the exact solution.
Now, finally to my question:
Is it possible to solve this problem whitin Fenicsx / PETSc, without having to manually extract the matricies to SciPy or NumPy and handle the dirichlet BCs manually?