I am working on a problem consisting of multiple equations that need to be solved on different subdomains. Everything works fine when I solve the coupled equations on the whole domain. As soon as I try to solve the equations on particular subdomains, I get an error that I cannot solve by myself.
To reproduce this error, I have written a code for solving the pure elasticity problem of a rectangle consisting of an outer and an inner subdomain. Here the equation is to be solved on the outer subdomain. See the following figure:
Consider the following the MWE:
from dolfinx import la
from dolfinx.fem import (dirichletbc, form, Function, FunctionSpace, locate_dofs_geometrical)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, create_matrix, create_vector, set_bc)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (create_rectangle, locate_entities, meshtags)
from mpi4py import MPI
import numpy as np
from petsc4py import PETSc
from ufl import (derivative, dx, grad, Identity, inner, lhs, Measure, rhs, sym, TestFunction, tr, TrialFunction, VectorElement)
#--------------------------Create mesh--------------------------
nxy = 50
mesh = create_rectangle(comm=MPI.COMM_WORLD, points=[np.array([-2, -2]), np.array([2, 2])], n=[nxy,nxy])
#--------------------------Boundary conditions--------------------------
def Top(x):
return np.isclose(x[1], 2)
def Bottom(x):
return np.isclose(x[1], -2)
def bc_l(x):
result1 = np.isclose(x[0], -2)
result2 = np.isclose(x[1], -2)
return np.logical_and(result1, result2)
u_el = VectorElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, u_el)
U_0, submap0= V.sub(0).collapse()
U_1, submap1 = V.sub(1).collapse()
dofs_u_top = locate_dofs_geometrical((V.sub(1), U_1), Top)
dofs_u_bot = locate_dofs_geometrical((V.sub(1), U_1), Bottom)
dofs_u_left = locate_dofs_geometrical((V.sub(0), U_0), bc_l)
u_Ly = Function(U_1)
with u_Ly.vector.localForm() as bc_local:
bc_local.set(0.0)
zero_u = Function(U_1)
with zero_u.vector.localForm() as bc_local:
bc_local.set(0.0)
bc_u_top = dirichletbc(u_Ly, dofs_u_top, V.sub(1))
bc_u_bot = dirichletbc(zero_u, dofs_u_bot, V.sub(1))
bc_u_left = dirichletbc(zero_u, dofs_u_left, V.sub(0))
bc_u = [bc_u_top, bc_u_bot, bc_u_left]
#--------------------------measure dx--------------------------
Tol = 1e-1
def outersub(x):
result1 = x[0] <= -1+Tol
result2 = x[0] >= 1-Tol
result3 = x[1] <= -1+Tol
result4 = x[1] >= 1-Tol
return np.logical_or(result1, np.logical_or(result2, np.logical_or(result3, result4)))
def innersub(x):
result1 = x[0] > -1
result2 = x[0] < 1
result3 = x[1] > -1
result4 = x[1] < 1
return np.logical_and(result1, np.logical_and(result2, np.logical_and(result3, result4)))
boundaries = [(0, outersub),
(1, innersub)]
cell_indices, cell_markers = [], []
for (marker, locator) in boundaries:
cells = locate_entities(mesh, mesh.topology.dim, locator)
cell_indices.append(cells)
cell_markers.append(np.full(len(cells), marker))
cell_indices = np.array(np.hstack(cell_indices), dtype=np.int32)
cell_markers = np.array(np.hstack(cell_markers), dtype=np.int32)
sorted_cells = np.argsort(cell_indices)
cell_tag = meshtags(mesh, mesh.topology.dim, cell_indices[sorted_cells], cell_markers[sorted_cells])
dx = Measure("dx", domain=mesh, subdomain_data=cell_tag)
#--------------------------Variational form--------------------------
#~~~~~~~~Introduce manually the material parameters
E = PETSc.ScalarType(15.96)
nu = PETSc.ScalarType(0.2)
lamda = PETSc.ScalarType((nu*E)/((1+nu)*(1-2*nu)))
mu = PETSc.ScalarType(E/(2*(1+nu)))
#~~~~~~~~Constitutive functions
def epsilon(u_):
return sym(grad(u_))
def sigma(u_):
return 2.0*mu*epsilon(u_)+lamda*tr(epsilon(u_))*Identity(len(u_))
#~~~~~~~~
u = Function(V, name="displacement")
v = TestFunction(V)
E_u = (inner(sigma(u), epsilon(v)))*dx(0)
#--------------------------create solver--------------------------
class NonlinearPDE_SNESProblem():
def __init__(self, F, J, soln_vars, bcs, P=None):
self.L = form(F)
self.a = form(J)
self.a_precon = P
self.bcs = bcs
self.soln_vars = soln_vars
def F(self, snes, x, F):
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
with x.localForm() as _x:
self.soln_vars.x.array[:] = _x.array_r
with F.localForm() as f_local:
f_local.set(0.0)
assemble_vector(F, self.L)
apply_lifting(F, [self.a], bcs=[self.bcs], x0=[x], scale=-1.0)
F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(F, self.bcs, x, -1.0)
def J(self, snes, x, J, P):
J.zeroEntries()
assemble_matrix(J, self.a, bcs=self.bcs, diagonal=1.0)
J.assemble()
if self.a_precon is not None:
P.zeroEntries()
assemble_matrix(P, self.a_precon, bcs=self.bcs, diagonal=1.0)
P.assemble()
#~~~~~~~~displacement problem
dU = TrialFunction(V)
a_u = lhs(E_u)
L_u = rhs(E_u)
F_u=a_u-L_u
J_u = derivative(F_u, u, dU)
F_u, J_u = form(F_u), form(J_u)
Jmat_u = create_matrix(J_u)
Fvec_u = create_vector(F_u)
problem_u = NonlinearPDE_SNESProblem(F=F_u, J=J_u, soln_vars=u, bcs=bc_u)
#~~~~~~~~KSP solver
b_u = la.create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs)
J_u = create_matrix(problem_u.a)
solver_u_snes = PETSc.SNES().create(MPI.COMM_WORLD)
solver_u_snes.setFunction(problem_u.F, Fvec_u)
solver_u_snes.setJacobian(problem_u.J, Jmat_u, P=None)
solver_u_snes.getKSP().setType("gmres")
solver_u_snes.getKSP().getPC().setType("lu")
#--------------------------solve--------------------------
with u_Ly.vector.localForm() as bc_local:
bc_local.set(-0.1)
solver_u_snes.solve(None, u.vector)
xdmf = XDMFFile(mesh.comm, "./MWE.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(u)
xdmf.close()
The source of the solver is: dolfinx/test_nonlinear_assembler.py at 611361eb24a6b8f6a9480b2034f53dc67feecf7e · FEniCS/dolfinx · GitHub
The output/error message reads:
Form has no parts with arity 2.
WARNING:UFL:Form has no parts with arity 2.
Traceback (most recent call last):
File "/root/shared/PhaseField/Different H's/Fenicsx/2D/MWEx.py", line 157, in <module>
solver_u_snes.solve(None, u.vector)
File "PETSc/SNES.pyx", line 590, in petsc4py.PETSc.SNES.solve
petsc4py.PETSc.Error: error code 73
[0] SNESSolve() at /usr/local/petsc/src/snes/interface/snes.c:4756
[0] SNESSolve_NEWTONLS() at /usr/local/petsc/src/snes/impls/ls/ls.c:222
[0] KSPSolve() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:1078
[0] KSPSolve_Private() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:843
[0] KSPSetUp() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:407
[0] PCSetUp() at /usr/local/petsc/src/ksp/pc/interface/precon.c:990
[0] PCSetUp_LU() at /usr/local/petsc/src/ksp/pc/impls/factor/lu/lu.c:93
[0] MatLUFactorSymbolic() at /usr/local/petsc/src/mat/interface/matrix.c:2966
[0] MatLUFactorSymbolic_SeqAIJ() at /usr/local/petsc/src/mat/impls/aij/seq/aijfact.c:304
[0] Object is in wrong state
[0] Matrix is missing diagonal entry 838
Dolfinx version (installed via Docker dolfinx/dolfinx): 0.3.1.0
petsc4py version: 3.17.0
Thank you in advance!
With kind regards
Dejan