Greetings! I’d like to thank you all for the question and answers, it gave me a good idea for a problem I’m working now. I tried to reproduce the method, specially with the jkrokowski code, but I’m currently struggling with an issue when solving the system with nested matrices with PETSc.
I was able to assemble the matrix and vectors, but when I try to solve the system with a CG type PETSc KSP solver it returns the following error:
Error: error code 60
[0] KSPSolve() at /tmp/petsc-src/src/ksp/ksp/interface/itfunc.c:1082
[0] KSPSolve_Private() at /tmp/petsc-src/src/ksp/ksp/interface/itfunc.c:910
[0] KSPSolve_CG() at /tmp/petsc-src/src/ksp/ksp/impls/cg/cg.c:248
[0] KSP_MatMult() at /tmp/petsc-src/include/petsc/private/kspimpl.h:347
[0] MatMult() at /tmp/petsc-src/src/mat/interface/matrix.c:2608
[0] MatMult_Nest() at /tmp/petsc-src/src/mat/impls/nest/matnest.c:51
[0] MatMultAdd() at /tmp/petsc-src/src/mat/interface/matrix.c:2770
[0] Nonconforming object sizes
[0] Mat mat,Vec v1: global dim 1056 16
It looks like is having some trouble with the nested structure when using the iterative methods, and I’m not seeing what exactly is missing. Below a sample code in Python I made which returns this error. I also let a Colab file with the code, in case it helps to confirm the versions in use Google Colab .
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 20, 20, dolfinx.mesh.CellType.quadrilateral)
#setting boundaries markers and indicator functions
#Electrodes e_j are displace in the middle of each boundary
boundaries = [
(0, lambda x: np.logical_and( np.isclose(x[0],0), np.logical_and(x[1]>=0.25,x[1]<=0.75))),
(1, lambda x: np.logical_and( np.isclose(x[0],1), np.logical_and(x[1]>=0.25,x[1]<=0.75))),
(2, lambda x: np.logical_and( np.isclose(x[1],0), np.logical_and(x[0]>=0.25,x[0]<=0.75))),
(3, lambda x: np.logical_and( np.isclose(x[1],1), np.logical_and(x[0]>=0.25,x[0]<=0.75)))
]
#creating facet tags
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = dolfinx.mesh.locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
#Creating measure object
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
V = dolfinx.fem.FunctionSpace(mesh, ("CG", 1)) #Lagrange space for u
V0 = dolfinx.fem.FunctionSpace(mesh, ("DG",0)) #Discontinuous Garlekin for gamma
N = V.dofmap.index_map.size_global; L = 4; z_value = 1
gamma = dolfinx.fem.Constant(mesh, PETSc.ScalarType(1))
z = dolfinx.fem.Constant(mesh, PETSc.ScalarType(z_value))
I = np.zeros(L)
I[0] = 1; I[1] = -1
u = ufl.TrialFunction(V) #function u_h
v = ufl.TestFunction(V)
#Assembling A1
A1_sum_array = [(1/z)*ufl.inner(u,v)*ds(i) for i in range(L)]
A1_form = (gamma * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + np.sum(A1_sum_array))
A1_form_obj = dolfinx.fem.form(A1_form)
A1_matrix = dolfinx.fem.petsc.assemble_matrix(A1_form_obj)
###
#Assembling A2
A2_form_array = [-ufl.inner((1/z), v) * ds(i) for i in range(L)]
A2_form_obj = [dolfinx.fem.form(form) for form in A2_form_array]
A2_vector_array= [dolfinx.fem.petsc.assemble_vector(form) for form in A2_form_obj]
A2_matrix = PETSc.Mat().create(comm=MPI.COMM_WORLD)
A2_matrix.setSizes([N,L])
for j in range(L):
A2_matrix.setValues(range(N),j,A2_vector_array[j].getValues(range(N)))
A2_matrix.assemble()
A2T_matrix = A2_matrix.transpose()
A2T_matrix.assemble()
###
#Assembling A4
cell_area_list = []
v = ufl.TrialFunction(V0) #trial object
for i in range(L):
s = v*ds(i) # creates object of the form \int_{e_i} 1*ds
cell_area_form = dolfinx.fem.form(s)
cell_area = dolfinx.fem.assemble_scalar(cell_area_form)
cell_area_list.append(cell_area.real)
cell_area_array = np.array(cell_area_list)
A4_diagonal = np.diag(cell_area_array*(1/z_value))
A4_matrix = PETSc.Mat().create(comm=MPI.COMM_WORLD)
A4_matrix.setSizes(L)
A4_matrix.setValues(range(L),range(L),A4_diagonal.ravel())
A4_matrix.assemble()
# Assembling matrix A
A = PETSc.Mat()
A.createNest([[A1_matrix, A2_matrix],
[A2T_matrix, A4_matrix]])
A.setUp()
A.assemble()
#Assembling b vector
zero_vec = PETSc.Vec().create(comm=MPI.COMM_WORLD)
zero_vec.setSizes(N)
zero_vec.setUp()
I_vec = PETSc.Vec().create(comm=MPI.COMM_WORLD)
I_vec.setSizes(L)
I_vec.setUp()
I_vec.setArray(I)
I_vec.assemble()
b_nested = PETSc.Vec()
b_nested.createNest([zero_vec,I_vec])
b_nested.setUp()
# Solving problem
#Create solution vector
u_vec = PETSc.Vec().create(comm=MPI.COMM_WORLD)
u_vec.setSizes(N)
U_vec = PETSc.Vec().create(comm=MPI.COMM_WORLD)
U_vec.setSizes(L)
u_vec.setUp()
U_vec.setUp()
u_nested = PETSc.Vec().create(comm=MPI.COMM_WORLD)
u_nested.createNest([u_vec, U_vec]) #x0 and x1 are petsc vectors for your solution
u_nested.setUp()
#Solving
ksp = PETSc.KSP().create(mesh.comm)
ksp.setType("cg")
ksp.setTolerances(rtol=1e-14)
ksp.setOperators(A)
ksp.setFromOptions()
ksp.solve(b_nested,u_nested)
For context, I’m working in the Electrical Impedance Tomography problem, more specific in the Complete Electrode Method. Here, we are looking to solve for a potential u:\Omega \to \mathbb R and respective noised measures U_1,...,U_L\in \mathbb R on the boundary. This model form looks like:
B((u,U),(v,V)) = \sum_{l=1}^L I_l V_l \\
where
B((u,U),(v,V)) = \iint_\Omega \gamma \nabla u \cdot \nabla v dx + \sum_{l=1}^L \frac{1}{z_l}\int_{E_l} (u-U_l)(v - V_l)d s ,
beeing u,v \in H^1(\Omega), U,V \in \mathbb R^L, \gamma \in L^\infty(\Omega), I_j, z_j \in \mathbb R and E_1,..., E_L a set of L integration domains on the boundary representing electrodes. As it needs to solve for a (u,U) \in H^1(\Omega)\times \mathbb R^L vector and there are no available spaces for \mathbb R^L domain, I tried to manually assemble the matrix in a block structure as suggested in this topic.
Given the discretized domain with N dofs and L integration domains, the assembled matrix has the block form
A = \begin{bmatrix}
A^{(1)} & A^{(2)} \\
(A^{(2)})^T & A^{(4)}
\end{bmatrix} \in \mathbb R^{(N+L) \times (N+L)}
with
A^{(1)}_{p,n} = \begin{bmatrix}
\left( \sum_{j=1}^{K} \gamma_j \int_{T_j} \nabla \Upsilon_p \cdot \nabla \Upsilon_n dx \right) + \left( \sum_{j=1}^{L} \frac{1}{z_j} \int_{E_j} \Upsilon_p \Upsilon_n dS \right)
\end{bmatrix} \in \mathbb{R}^{N \times N}
A^{(2)}_{p,n} = \begin{bmatrix}
\left( \frac{1}{z_n} \int_{E_n} \Upsilon_p dS \right)
\end{bmatrix}\in \mathbb{R}^{N \times L}
A^{(4)}_{p,k} = \begin{bmatrix}
\left( \frac{1}{z_p} \int_{E_p} dS \right)
\end{bmatrix}\in \mathbb{R}^{L \times L}
beeing \Upsilon the function element basis.
Thank you very much for the help!
