 # Custom NewtonSolver using the mixed-dimensional branch and PETScNestMatrix

Hello everyone,

I’m trying to use the PETScNestMatrix to solve a double-saddle point problem. I precompute one of the blocks to a PETScMatrix. I can’t use directly the mixed-dimensional branch’s solvers because I need to add this precomputed matrix directly to the system. I started by trying to solve an incompressible elasticity problem. I did something similar to this to define a Problem where I could manipulate the assembly of the matrix and create the PETScNestMatrix.

``````class Problem(d.NonlinearProblem):
def __init__(self, Ji, Fi):
self.bilinear_forms = Ji
self.linear_forms = Fi
d.NonlinearProblem.__init__(self)

def F(self, b, x):
bs = []
for Gi in self.linear_forms:
bs_ = d.PETScVector()
d.assemble(Gi, tensor = bs_)
bs.append(bs_)

b_ = bs

self.A.init_vectors(b, b_)

def J(self):
Amat = d.PETScMatrix()
Ai = []
for dg in self.bilinear_forms:
Amat = d.assemble(dg)
Amat.apply('insert')
Ai.append(Amat)

Amat = d.PETScNestMatrix(Ai)
Amat.apply('insert')

self.A = Amat
return Amat
``````

It’s usual that I passed the matrix `A` to the function `J`, however, I couldn’t use the option `set_nest` that it’s used in the C++ code of the mixed-dimensional branch. So my first question is: Is it possible to create an empty PETScNestMatrix and then pass to it the list with the PETScMatrix from Python?

I tried to solve it anyway and I created a NewtonSolver following the C++ codes of the mixed-dimensional branch and I tested in a simple problem.

``````class NewtonSolver(d.NewtonSolver):
def __init__(self, mesh):
factory = d.PETScFactory.instance()
comm = mesh.mpi_comm()
d.NewtonSolver.__init__(self, comm, d.PETScLUSolver(),
factory)

self.solver = self.linear_solver()

def solve(self, problem, u):
matP = d.PETScMatrix()

# Extract parameters
convergence_criterion = self.parameters['convergence_criterion']
maxiter = self.parameters['maximum_iterations']
relaxation_parameter = 1.

# Reset iteration counts
newton_iteration = 0
krylov_iterations = 0

# Compute J to initialize the vectors
matA = problem.J()

x = d.PETScVector()
matA.init_vectors(x, u)

b = d.PETScVector()
problem.F(b, x)
dx = b.copy()
dx.zero()

# Check convergence
newton_converged = False
if convergence_criterion == 'residual':
newton_converged = self.converged(b, problem, 0)
elif convergence_criterion == 'incremental':
newton_converged = False
else:
raise ValueError("The convergence criterion %s is unknown, known criteria are 'residual' or 'incremental'")

# Start iterations
while (not newton_converged and newton_iteration < maxiter):
# compute Jacobian
matA = problem.J()
matA.convert_to_aij()

# Setup (linear) solver (including set operators)
self.solver_setup(matA, matP, problem, newton_iteration);

# Perform linear solve and update total number of Krylov iterations
if not dx.empty():
dx.zero()

dx0 = np.linalg.solve(matA.array(), b.get_local())
print('numpy solver:', np.max(dx0))

krylov_iterations += self.solver.solve(dx, b);
print('PETSc solver:', np.max(dx.get_local()))

# Update solution
self.update_solution(x, dx, relaxation_parameter, problem, newton_iteration);
newton_iteration += 1

# Compute F
matA = problem.J()
problem.F(b, x);

# Test for convergence
newton_converged = self.converged(b, problem, newton_iteration)

if newton_converged:
if self.mpi_comm.rank == 0:
d.info("Newton solver finished in %d" % newton_iteration +
" iterations and %d" % krylov_iterations + " linear solver iterations.")
else:
error_on_nonconvergence = self.parameters["error_on_nonconvergence"]
if error_on_nonconvergence:
if newton_iteration == maxiter:
raise ValueError("NewtonSolver.cpp",
"solve nonlinear system with NewtonSolver",
"Newton solver did not converge because maximum number of iterations reached")

else:
raise ValueError("NewtonSolver.cpp",
"solve nonlinear system with NewtonSolver",
"Newton solver did not converge")
else:
print("Newton solver did not converge.")

return newton_iteration, newton_converged

def update_solution(self, x, dx, relaxation_parameter, problem, iteration):
if relaxation_parameter == 1.0:
x -= dx
else:
x.axpy(-relaxation_parameter, dx)

#%% Parameters
mu = d.Constant(1)
lmbda = d.Constant(1)

h1 = d.Constant((0, 0, 0.1))
h2 = d.Constant((0, 0, -0.1))

n = 3
mesh = d.UnitCubeMesh(n, n, n)

#%% Faces
class Top(d.SubDomain):
def inside(self, x, on_boundary):
return d.near(x, 1.)

class Bottom(d.SubDomain):
def inside(self, x, on_boundary):
return d.near(x, 0.)

top = Top()
bot = Bottom()

marker = d.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
bot.mark(marker, 1)
top.mark(marker, 2)

#%% Spaces
V = d.VectorFunctionSpace(mesh, 'CG', 2)
P = d.FunctionSpace(mesh, 'DG', 1)

W = d.MixedFunctionSpace(V, P)

w = d.Function(W)
u, p = w.sub(0), w.sub(1)
(du, dus) = d.TrialFunctions(W)
(delta_u, delta_us) = d.TestFunctions(W)

dS = d.Measure("ds", domain = W.sub_space(0), subdomain_data = marker)

#%% Kinematics
dim = len(u)
I = d.Identity(dim)             # Identity tensor

J  = d.det(F)
FiT = d.inv(F.T)
C = F.T*F/J**(2./3.)
I1C = d.tr(C)

# Strain energy function
Psi = mu/2.*(I1C - 3.) - p*(J-1)
Pi = Psi*d.dx + d.inner(h1, u)*dS(1) + d.inner(h2, u)*dS(2)
Gi = [d.derivative(Pi, u), d.derivative(Pi, p)]

dGi = []
for gi in Gi:
for wi in w.split():
dGi.append(d.derivative(gi, wi))

#%% Set solver
problem = Problem(dGi, Gi)
solver = NewtonSolver(mesh)
solver.parameters['maximum_iterations'] = 1
solver.parameters['error_on_nonconvergence'] = False
solver.solve(problem, [u.vector(), p.vector()])
``````

When I try to solve this problem, in the first iteration (and the followings) the PETScLUSolver gave me a different value than, if for example, I solved it with numpy:

``````numpy solver: 29.3059696871
PETSc solver: 10856.9383564
``````

If I use a PETScKrylovSolver results in this 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
``````

Any idea of why this is happening? I’m using the latest docker container of the mixed-dimensional branch.

Thank you,

Hello,

There is indeed a `set_nest` function in the C++ PETScMatrix class, taking as parameters a list of PETSc matrices, which will create or update the corresponding PETScNestMatrix. This C++ function doesn’t seem to be accessible from the python interface.

Regarding the NewtonSolver and PETScKrylovSolver, it is difficult to figure out what is the problem without a MWE.

Hello javijv4,

maybe you can setup your nested matrix directly using petsc4py.
I am using something like this:

``````M = PETSc.Mat().createNest([[M00,M01], [M10,M11]], comm=MPI.COMM_WORLD)
``````

where all the `M` as well as all sub-matrices are backend_type matrices (created with e.g. `as_backend_type(X).mat`)

best wishes
Florian

@cdaversin, thank you for the fast answer. Here is a MEW:

``````import dolfin as d
import numpy as np

mu = d.Constant(1)
lmbda = d.Constant(1)
h1 = d.Constant((0, 0, 0.1))
h2 = d.Constant((0, 0, -0.1))

mesh = d.UnitCubeMesh(3, 3, 3)
class Top(d.SubDomain):
def inside(self, x, on_boundary):
return d.near(x, 1.)
class Bottom(d.SubDomain):
def inside(self, x, on_boundary):
return d.near(x, 0.)
top, bot = Top(), Bottom()

marker = d.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
bot.mark(marker, 1)
top.mark(marker, 2)

# Spaces
V = d.VectorFunctionSpace(mesh, 'CG', 2)
P = d.FunctionSpace(mesh, 'DG', 1)
W = d.MixedFunctionSpace(V, P)

w = d.Function(W)
u, p = w.sub(0), w.sub(1)
(delta_u, delta_us) = d.TestFunctions(W)

dS = d.Measure("ds", domain = W.sub_space(0), subdomain_data = marker)

# Kinematics
I = d.Identity(3)             # Identity tensor
J  = d.det(F)
I1C = d.tr(F.T*F/J**(2./3.))

Psi = mu/2.*(I1C - 3.) - p*(J-1)        # Strain energy function
Pi = Psi*d.dx + d.inner(h1, u)*dS(1) + d.inner(h2, u)*dS(2)
Gi = [d.derivative(Pi, u), d.derivative(Pi, p)]

dGi = []
for gi in Gi:
for wi in w.split():
dGi.append(d.derivative(gi, wi))

# Assemble bilinear form
Ai = []
for dg in dGi:
Ai.append(d.assemble(dg))
Amat = d.PETScNestMatrix(Ai)
Amat.apply('insert')

# Assemble linear form
bs = []
for g in Gi:
bs.append(d.assemble(g))

b = d.PETScVector()
Amat.init_vectors(b, bs)

# Set solver
solver_type = 'lu'
if solver_type == 'krylov':
solver = d.PETScKrylovSolver()
elif solver_type == 'lu':
solver = d.PETScLUSolver()

Amat.convert_to_aij()
solver.set_operator(Amat)

# Solve
dx0 = np.linalg.solve(Amat.array(), b.get_local())
print('numpy solver:', np.max(dx0))

dx = b.copy()
dx.zero()
solver.solve(dx, b);
print('PETSc solver:', np.max(dx.get_local()))
``````

Like I said before, if I set `solver_type = 'krylov'`, then I get the PESTc error 73. If I set `solver_type = 'lu'`, then I get a different result depending if I solve using numpy or the PETSc solver.

@Florian_Bruckner yes, I know I can do that. However, if I want to create the matrix and the use a function to give it the submatrices it doesn’t work. For example:

``````import dolfin as d
from petsc4py import PETSc

def Jmatrix(A, dGi):
Ai = []
for dg in dGi:
Amat = d.assemble(dg)
Amat.apply('insert')
Ai.append(d.as_backend_type(Amat).mat())

A, B, Bt, C = Ai
Amat = PETSc.Mat().createNest([[A,B], [Bt,C]])

mu = d.Constant(1)
lmbda = d.Constant(1)

mesh = d.UnitCubeMesh(3, 3, 3)
V = d.VectorFunctionSpace(mesh, 'CG', 2)
P = d.FunctionSpace(mesh, 'DG', 1)

W = d.MixedFunctionSpace(V, P)

w = d.Function(W)
u, p = w.sub(0), w.sub(1)
(du, dus) = d.TrialFunctions(W)
(delta_u, delta_us) = d.TestFunctions(W)

dim = len(u)
I = d.Identity(dim)             # Identity tensor

J = d.det(F)
FiT = d.inv(F.T)
C = F.T*F/J**(2./3.)
I1C = d.tr(C)

# Strain energy function
Psi = mu/2.*(I1C - 3.) - p*(J-1)
Pi = Psi*d.dx
Gi = [d.derivative(Pi, u), d.derivative(Pi, p)]

dGi = []
for gi in Gi:
for wi in w.split():
dGi.append(d.derivative(gi, wi))

A = d.PETScNestMatrix()
Jmatrix(A, dGi)

print(A.array().shape)
``````

This will print out (0,0). It is useful to do give the values to `A` inside a Function because this way I only create the matrix once and then I just update its values. Anyway, this is convenient but I can definitely work around that.

My guess is that the PETSc solver and the numpy solver you are using are different, and might then give different results. I am not a PETSc expert, but I did a quick test from your MWE with a different PETSc solver configuration. For example, using `superlu` :

``````solver = d.PETScLUSolver()
ksp = solver.ksp()
ksp.setType('preonly')
pc = ksp.getPC()
pc.setType('lu')
pc.setFactorSolverPackage('superlu')
ksp.setFromOptions()
``````

gives :

``````numpy solver: 13.312943321351039
PETSc solver: 287.7411886846228
``````

I am not sure what is behind the numpy solver, and I think your question is closely related to PETSc so it could be worth posting it to the PETSc mailing list.