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[2], 1.)
class Bottom(d.SubDomain):
def inside(self, x, on_boundary):
return d.near(x[2], 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
F = I + d.grad(u) # Deformation gradient
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,