Preconditioners are not available in mixed-dimensional branch

1

Dear all,
Working with mixed dimensional branch and demos like 2D-2D block assembly i’ve noticed that preconditioners are not available yet for mixed variational forms. For example:

#Create mesh and define function space
mesh = UnitSquareMesh(100, 100)

marker = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
for c in cells(mesh):
    marker[c] = c.midpoint().y() < 0.5
    #marker[c] = c.midpoint().x() < 0.5

submesh1 = MeshView.create(marker, 1)
submesh2 = MeshView.create(marker, 0)

# Define Dirichlet boundary
def boundarySub1(x):
    return x[1] < DOLFIN_EPS
    #return x[0] < DOLFIN_EPS

def boundarySub2(x):
    return x[1] > 1.0 - DOLFIN_EPS
    #return x[0] > 1.0 - DOLFIN_EPS

#element2D = FiniteElement("Lagrange", triangle, 1)
W1 = FunctionSpace(submesh1, "Lagrange", 1)
W2 = FunctionSpace(submesh2, "Lagrange", 1)
# Define the mixed function space
V = MixedFunctionSpace( W1, W2 )

# Define boundary conditions
u0 = Constant(0.0)
# Subdomain 1
bc1 = DirichletBC(V.sub_space(0), u0, boundarySub1)
# Subdomain 2
bc2 = DirichletBC(V.sub_space(1), u0, boundarySub2)
bcs = [bc1,bc2]

# Define variational problem
# Use directly TrialFunction and TestFunction on the product space
# Subdomain 1
u1 = TrialFunction(W1)
v1 = TestFunction(W1)
# Subdomain 2
u2 = TrialFunction(W2)
v2 = TestFunction(W2)
# Mixed
(u1_m,u2_m) = TrialFunctions(V)
(v1_m,v2_m) = TestFunctions(V)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)

# Mixed weak form
a = inner(grad(u1_m), grad(v1_m))*dx + inner(grad(u2_m), grad(v2_m))*dx
L = f*v1_m*dx + f*v2_m*dx
# Solve the problem
sol = Function(V)
solve(a == L, sol, bcs, solver_parameters={"linear_solver":"gmres","preconditioner":"ilu"})

returns the error:

 Error:   Unable to successfully call PETSc function 'KSPSolve'.

Similar errors occured for other preconditioners. Changing preconditioner to ‘‘none’’ solved the problem.Is there any workaround to this ;

2

Hello,

I was able to reproduce your error, and it seems that some preconditioners are not compatible with the PETSc MATNEST format (which is the format used to define the block shaped matrices).
The assembled matrix can be converted to the more standard MATAIJ type using the function :

A.convert_to_aij()

where A is your matrix of type PETScNestMatrix.

When using solve(...) as in your example, you don’t directly manipulate your matrix, since this function does the block extraction and the block-by-block assembly for you. Which corresponds to :

a_list = extract_blocks(a)
A_blocks = [assemble_mixed(a) for a in a_list]
A = PETScNestMatrix(A_blocks)

In the very last version of the framework (from today), the function convert_to_aij() is used by solve(...) when the preconditioner is not the default one as in your case (the latest Docker container is up to date). Hope that fixes the issue.

2 Likes
3

Thank you @cdaversin,
I pulled the latest docker container and is working perfectly. Is there any demo/tutorial that contains the process of manual assembly of block matrices, application of dirichlet bcs and PETSc solver call ;

4

Is there any demo/tutorial that contains the process of manual assembly of block matrices, application of dirichlet bcs and PETSc solver call ;

I believe you can find some information here

2 Likes
5

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?

6

Hello,

You can configure the non linear solver as you did, through solver.parameters.
I was able to reproduce your issue with the linear_solver configuration, and I updated the mixed-dimensional branch together with the latest Docker container so that this should be fixed.

2 Likes
7

Dear @cdaversin

Thank you very much for your help. Unfortunately, I have pulled the latest version of mixed-dimensional branch and this time I get new following error

*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 92 (See https://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html for possible LU and Cholesky solvers).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  53c3cc165c012a37dca95e6f90da7b6ee05d5156

And this time when I solve the problem using the direct method, which is solvable with previous version of the branch, I also get 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.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  f0a90b8fffc958ed612484cd3465e59994b695f9

The MWE that produce the latter error is

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'] = 1.0
solver.solve()