Hello @cdaversin
I am trying to solve a coupled problem and also face this error although I have updated the branch to the latest version.
*** 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.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: f0a90b8fffc958ed612484cd3465e59994b695f9
The MWE is as follwing,
from dolfin import *
import matplotlib.pyplot as plt
from ufl.algorithms import expand_derivatives
# define the edges of the domain
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
class LeftF(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and between(x[1], (0.7, 1))
class RightF(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0) and between(x[1], (0.7, 1))
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.7)
class LeftS(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and between(x[1], (0.0, 0.7))
class RightS(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0) and between(x[1], (0.0, 0.7))
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
# Create meshes on a unit square 2D domain
mesh = UnitSquareMesh(100, 100)
# Define the region for the fluid as the upper part of the domain
class Fluid(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.7
fluid = Fluid()
# Define the the region for the solid as the lower part of the domain
class Solid(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.7
solid = Solid()
# numbering the regions
regions = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
fluid.mark(regions, 1)
solid.mark(regions, 2)
# numbering the interface edge
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)
Interface().mark(boundaries,1)
# create submeshes for fluid, solid and interface
mesh_F = MeshView.create(regions, 1)
mesh_S = MeshView.create(regions, 2)
mesh_I = MeshView.create(boundaries, 1)
# Define Function spaces
U_F = VectorFunctionSpace(mesh_F, "CG", 2) # Velocity on the fluid domain
P_F = FunctionSpace(mesh_F, "CG", 1) # Pressure on the fluid domain
T_F = FunctionSpace(mesh_F, "CG", 1) # Temperature on the fluid domain
T_S = FunctionSpace(mesh_S, "CG", 1) # Temperature on the solid domain
L_I_temp = FunctionSpace(mesh_I, "CG", 1) # Lagrange variable for temperature on the interface
# Define the mixed function space W = W1 x W2
W = MixedFunctionSpace(U_F, P_F, T_F, T_S, L_I_temp)
# Define markers for Dirichlet and Neuman bcs
ff_F = MeshFunction("size_t", mesh_F, mesh_F.geometry().dim()-1, 0)
Top().mark(ff_F, 1) # Top edge of fluid mesh
LeftF().mark(ff_F, 2) # Left edge of fluid mesh
RightF().mark(ff_F, 3) # Right edge of fluid mesh
Interface().mark(ff_F, 4) # Bottom edge of fluid mesh (interface)
ff_S = MeshFunction("size_t", mesh_S, mesh_S.geometry().dim()-1, 0)
Interface().mark(ff_S, 5) # Top edge of solid mesh (interface)
LeftS().mark(ff_S, 6) # Left edge of solid mesh
RightS().mark(ff_S, 7) # Right edge of solid mesh
Bottom().mark(ff_S, 8) # Bottom edge of solid mesh
p_left = 40.0
p_right = 0.0
pressure_left = Constant(p_left)
pressure_right = Constant(p_right)
# Define boundary conditions for fluid part
bc1 = DirichletBC(U_F, Constant((0.0,0.0)), ff_F, 1)
bc2 = DirichletBC(U_F, Constant((0.0,0.0)), ff_F, 4)
bc3 = DirichletBC(T_F, Constant(300), ff_F, 2)
bcU = [bc1,bc2,bc3]
# Define boundary conditions for solid part
bc4 = DirichletBC(T_S, Constant(310), ff_S, 8)
bcL = [bc4]
# Define boundary conditions for the entire part
bcs = [bc1,bc2,bc3,bc4]
# define trial functions and test functions
(uf,pf,tf,ts,l_temp) = TrialFunctions(W)
(vf,qf,zf,zs,e_temp) = TestFunctions(W)
# define functions (unknowns in the weak form)
u = Function(W)
uf = u.sub(0)
pf = u.sub(1)
tf = u.sub(2)
ts = u.sub(3)
l_temp = u.sub(4)
# define the domains for integration
dx_F = Measure("dx", domain=W.sub_space(0).mesh())
ds_F = Measure("ds", domain=W.sub_space(0).mesh(), subdomain_data=ff_F)
dx_S = Measure("dx", domain=W.sub_space(3).mesh())
ds_S = Measure("ds", domain=W.sub_space(3).mesh(), subdomain_data=ff_S)
ds_I = Measure("dx", domain=W.sub_space(4).mesh())
# Fluid parameters for the variational problem
mu = 1.002
rho = 998
k = 0.5984
cp = 4181
alp = k/rho/cp
# Solid parameters for the variational problem
muS = 1
rhoS = 2705
kS = 227
cpS = 900
alpS = kS/rhoS/cpS
nf = FacetNormal(mesh_F) # unit normal vector on the surfaces of the fluid mesh
ns = FacetNormal(mesh_S) # unit normal vector on the surfaces of the solid mesh
# weak form for the fluid mesh (coupled equation of NSE, Continuity and Temperature Transport equations)
F1 = rho*inner(grad(uf)*uf, vf)*dx_F \
+ mu*inner(grad(uf), grad(vf))*dx_F \
- pf*div(vf)*dx_F \
+ pressure_left*inner(nf, vf)*ds_F(2) \
+ pressure_right*inner(nf, vf)*ds_F(3) \
+ qf*div(uf)*dx_F \
+ alp*inner(grad(tf), grad(zf))*dx_F \
+ inner(dot(uf,grad(tf)),zf)*dx_F \ # weak form for the solid mesh (Diffusion equation for temperature)
F2 = alpS*inner(grad(ts), grad(zs))*dx_S # diffusion term for temperature transport equation
# weak form for the whole problem (combine the weak form of fluid and solid and add lagrange multiplier at the interface)
F = F1 + F2 + inner(l_temp, zs-zf)*ds_I + inner(ts-tf, e_temp)*ds_I
# solving the problem
F_blocks = extract_blocks(F)
Jacobian = []
for Fi in F_blocks:
for ui in u.split():
d = derivative(Fi, ui)
d = expand_derivatives(d)
Jacobian.append(d)
problem = MixedNonlinearVariationalProblem(F_blocks, u.split(), bcs, J=Jacobian)
solver = MixedNonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-5
prm['newton_solver']['relative_tolerance'] = 1E-5
prm['newton_solver']['maximum_iterations'] = 50
prm['newton_solver']['relaxation_parameter'] = 0.8
prm['newton_solver']['linear_solver'] = 'gmres'
prm['newton_solver']['preconditioner'] = 'ilu'
prm['newton_solver']['krylov_solver']['absolute_tolerance'] = 1E-9
prm['newton_solver']['krylov_solver']['relative_tolerance'] = 1E-9
prm['newton_solver']['krylov_solver']['maximum_iterations'] = 1000
prm['newton_solver']['krylov_solver']['monitor_convergence'] = True
prm['newton_solver']['krylov_solver']['nonzero_initial_guess'] = False
solver.solve()
I am not sure what is the problem but this code works with the direct solver.
Also, If I need to manipulate the solver by myself e.g. like in this post, do you have any example of doing so in the mixed dimensional branch?