Help to solve a contour problem: PETSc error 73

variational formulation
1

Hello guys. I really appreciate with someone could help me.

I’m trying to solve the following contour problem:

\varsigma \int_{\Gamma_s} \nabla : \sigma^s (u) ds = \int_{\Gamma_s} w P.n dx

using the following structure: (giving just the essential)

class Interface(SubDomain):
def inside(self, x, on_boundary):
dy = x[1]
dz = x[2]
dr = (dydy + dzdz)**0.5
zeta = DOLFIN_EPS + 1.0E-3
return dr >( ri - zeta) and between(x[0], (0.2,4.0))

parts = MeshFunction(“size_t”, mesh, mesh.topology().dim()-1, 0)
Interface().mark(parts,4)
ds = Measure(‘ds’, domain=mesh,subdomain_data=parts)

timeseries_p = TimeSeries(‘pressure_series’)
timeseries_p.retrieve(pr.vector(), t)

pavg = (assemble(pr*ds(4))/(assemble(Constant(1.0)*ds(4)))) #mean pressure in the interface
w1 = TestFunction(W.sub(0))
u1 = TrialFunction(W.sub(0))
us = Function(W.sub(0).collapse())

Fhp = inner(pavg*n, w1)ds(4) - qsiinner(sigma_s(u1),epsilon_s(w1))*ds(4)
ahp, Lhp = system(Fhp)
Ahp, bhp = assemble(ahp), assemble(Lhp)
solve(Ahp, us.vector(), bhp, ‘gmres’, ‘ilu’)

and I’m stuck with the following error:

Error: Unable to successfully call PETSc function ‘KSPSolve’.
*** Reason: PETSc error code is: 73 (Object is in wrong state).
*** Where: This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.

Someone could help me?

2

Please provide enough code to make reproducing the error possible.

3

I Think I’ve found the problem.
I guess it is a bug, because I’m using a mixed function space.

V = VectorElement(“Lagrange”, tetrahedron, 1) #Velocity vector element
P = FiniteElement(“Lagrange”, tetrahedron, 1) #Pressure scalar element
VP = MixedElement([V,P]) #mixed element
W = FunctionSpace(mesh, VP) #mixed function space

So I created a new VectorFunctionSpace to be equal to the subspace(0)

A = VectorFunctionSpace(mesh, ‘P’, 1) = W.sub(0)

and using

w1 = TestFunction(A)
u1 = TrialFunction(A)
us = Function(A)

I was able to solve

Fhp = inner(pavg*n, w1) ds(4) - qsi inner(sigma_s(u1),epsilon_s(w1))*ds(4)
ahp, Lhp = system(Fhp)
Ahp, bhp = assemble(ahp), assemble(Lhp)
solve(Ahp, us.vector(), bhp, ‘gmres’, ‘ilu’)

with no errors.

Should I report this bug? Or I’m missing something?

Thx