# Error using PETSc SNES

Hi,
I am trying to solve a incompressible hyperelasticity Problem with a SNES solver. I always got an errormessage and i dont know why. Maybe someone of you can help me please. The Error is:
`src/petsc4py.PETSc.c:352134: __Pyx_PyCFunction_FastCall: Assertion !PyErr_Occurred()' failed.`

Here is my code:

``````"""Unit tests for Newton solver assembly"""

import numpy as np
import math
import ufl
from dolfinx import cpp as _cpp
from dolfinx import fem, la, mesh, plot, io, nls, log, cpp
from mpi4py import MPI
from petsc4py import PETSc

class NonlinearPDE_SNESProblem:
def __init__(self, F, J, u, bcs):
self.L = F
self.a = J
self.bcs = bcs
self.u = u

# Create matrix and vector to be used for assembly
# of the non-linear problem
self.A = fem.petsc.create_matrix(self.a)
self.b = fem.petsc.create_vector(self.L)

def F(self, snes, x, b):
"""Assemble residual vector."""
x.copy(self.u.vector)

with b.localForm() as b_local:
b_local.set(0.0)
fem.assemble_vector(b, self.L)
fem.apply_lifting(b, [self.a], [self.bcs], [x],-1.0)
fem.set_bc(b, self.bcs, x, -1.0)

def J(self, snes, x, A, P):
"""Assemble Jacobian matrix."""
A.zeroEntries()
fem.assemble_matrix(A, self.a, self.bcs)
A.assemble()

def nonlinear_pde_snes():
"""Test Newton solver for a simple nonlinear PDE"""
c1 = PETSc.ScalarType(0.15)

def P(u,p):
# Kinematics
d = len(u)
I = ufl.variable(ufl.Identity(d))
C = ufl.variable(F.T * F)
I1 = ufl.variable(ufl.tr(C))
I2 = ufl.variable(ufl.tr(C*C))
J  = ufl.variable(ufl.det(F))
# Stored strain energy density
psi = c1*(I1-3)
# Hyper-elasticity
return ufl.diff(psi, F) + p * J * ufl.inv(F.T)

# Create mesh and function space
comm = MPI.COMM_WORLD
rank = comm.rank
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [1, 1, 1]], [10,10,10], mesh.CellType.hexahedron)
P2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)

state_space = fem.FunctionSpace(domain, P2 * P1)
V, _ = state_space.sub(0).collapse()

def left(x):
return np.isclose(x[0], 0)

def right(x):
return np.isclose(x[0], 1)

fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

u0_bc = fem.Function(V)
u0 = lambda x: np.zeros_like(x, dtype=PETSc.ScalarType)
u0_bc.interpolate(u0)

class incremented_displacement_expression:
def __init__(self):
self.t = 0.0

def eval(self, x):
# linearly incremented displacement
return np.stack((np.full(x.shape[1], self.t * 0.01), np.zeros(x.shape[1]), np.zeros(x.shape[1])))

incremented_displacement_expr = incremented_displacement_expression()
incremented_displacement_expr.t = 0
incremented_displacement = fem.Function(V)
incremented_displacement.interpolate(incremented_displacement_expr.eval)

left_dofs = fem.locate_dofs_topological((state_space.sub(0), V), facet_tag.dim, facet_tag.find(1))
right_dofs = fem.locate_dofs_topological((state_space.sub(0), V), facet_tag.dim, facet_tag.find(2))
bc_left = fem.dirichletbc(u0_bc, left_dofs, state_space.sub(0))
bc_right = fem.dirichletbc(incremented_displacement, right_dofs, state_space.sub(0))
bcs = [bc_left,bc_right]

#Define TestFunction
state = fem.Function(state_space, name="state")
test_state = ufl.TestFunctions(state_space)
u, p = ufl.split(state)
v, q = test_state

d = len(u)
I = ufl.variable(ufl.Identity(d))
det_F  = ufl.variable(ufl.det(F))

#Variational Problem
R = ufl.inner(ufl.grad(v), P(u,p))*dx + q * (det_F - 1) * dx

residual = fem.form(R)
Jac = ufl.derivative(R, state)
jacobian = fem.form(Jac)

# Create nonlinear problem
problem = NonlinearPDE_SNESProblem(residual,jacobian, state, bcs)
b = la.create_petsc_vector(state_space .dofmap.index_map, state_space .dofmap.index_map_bs)
J = fem.petsc.create_matrix(problem.a)

# Create Newton solver and solve
snes = PETSc.SNES().create(domain.comm)
opts = PETSc.Options() # read the prefix of solver
opts['snes_linesearch_type'] = 'bt'
opts['snes_monitor'] = None
opts['snes_linesearch_monitor'] = None
snes.setFromOptions() # make above settings effective

snes.setFunction(problem.F,b)
snes.setJacobian(problem.J,J)
snes.setTolerances(rtol=1.0e-9, max_it=10)
snes.getKSP().setType("preonly")
snes.getKSP().setTolerances(rtol=1.0e-9)
snes.getKSP().getPC().setType("lu")

for n in range(0,50):
incremented_displacement_expr.t = n
incremented_displacement.interpolate(incremented_displacement_expr.eval)
snes.solve(None,state.vector)
state.x.scatter_forward()

nonlinear_pde_snes()
``````

Instead of using the `fem.assemble*`, `fem.apply_lifting`, you should use `fem.petsc.*`, i.e.

``````
def F(self, snes, x, b):
"""Assemble residual vector."""
x.copy(self.u.vector)

with b.localForm() as b_local:
b_local.set(0.0)
fem.petsc.assemble_vector(b, self.L)
fem.petsc.apply_lifting(b, [self.a], [self.bcs], [x],-1.0)
fem.petsc.set_bc(b, self.bcs, x, -1.0)

def J(self, snes, x, A, P):
"""Assemble Jacobian matrix."""
A.zeroEntries()
fem.petsc.assemble_matrix(A, self.a, self.bcs)
A.assemble()
``````

When a `PETSc.Vec` finished its computation (e.g. matrix-vector multiplication and it is `x`, `b` and `self.u.vector` in the demo), it needs to update ghost with `ghostUpdate`. What `addv` and `mode` are to be used as you used different options here.

It would be `insert` and `forward`.

Tutorial of @jackhale DOLFINx in Parallel with MPI — NewFrac FEniCSx Training introduces some about two modes of `ghostUpdate`, but I believe there is still some misunderstanding on the `ghostUpdate` in my FSI implementation.

Here I have two overlapped fluid and solid meshes where functionspaces are built respectively. `I` is my parallel interpolation matrix (`PETSc.Mat`), where the interpolation between two functions can be accomplished by

``````I.mult(f.vector, s.vector)
``````

In my implementation, I assemble matrices and vectors on two meshes respectively with the names `A`, `b`, `As`, `bs`. Transformation can be done and solver is called by

``````A += I.T * As * I       # 3 matrices multiplication
b += I.T * bs

solver.setOperators(A)
solver.solve(b, U.vector)
``````

The following is the part of my attempt but with some issue as the whole code is too long. It can run but the result does not meet the expectation. I can sense that is due to `ghostUpdate` and look for help with the setting here:

``````from dolfinx.fem.petsc import (assemble_matrix, assemble_vector,
apply_lifting, set_bc)

A = assemble_matrix(a, bcs)
A.assemble()
b = assemble_vector(L)
apply_lifting(b, [a], bcs=[bcs])
set_bc(b, bcs)

As = assemble_matrix(as)
As.assemble()
bs = assemble_vector(Ls)

delA = I.transposeMatMult(As).matMult(I)
A += delA

delb = b.copy()
delb.zeroEntries()
I.multTranspose(bs, delb)
b += delb

solver.setOperators(A)
solver.solve(b, U.vector)
U.x.scatter_forward()
``````

I think I can workaround the above one (I think?) by wrapping the `PETSc.Vec` into a function like this:

``````from dolfinx.fem.petsc import (assemble_matrix, assemble_vector,
apply_lifting, set_bc)

A = assemble_matrix(a, bcs)
A.assemble()
b = assemble_vector(L)
apply_lifting(b, [a], bcs=[bcs])
set_bc(b, bcs)

b_fun = fem.Function(mixed_function_space)
b_fun.x.array[:mixed_function_space.dofmap.index_map.size_local] = b.getArray()
b_fun.x.scatter_forward()

As = assemble_matrix(as)
As.assemble()
bs = assemble_vector(Ls)
delA = I.transposeMatMult(As).matMult(I)
A += delA
delb = b.copy()
delb.zeroEntries()
I.multTranspose(bs, delb)

delb_fun = fem.Function(mixed_function_space)
delb_fun.x.array[:mixed_function_space.dofmap.index_map.size_local] = delb.getArray()
delb_fun.x.scatter_forward()

b_fun.x.array[:mixed_function_space.dofmap.index_map.size_local] += delb_fun.x.array[:mixed_function_space.dofmap.index_map.size_local]
b_fun.x.scatter_forward()
solver.setOperators(A)
solver.solve(b_fun.vector, U.vector)
``````

Try to answer my own question: `PETSc.Vec.ghostUpdate` should be called after vector assemble and transformation, i.e. matrix-vector multiplication. In my implementation, the ghost of `b`, `bs` and `delb` are update-to-date. To fix this, delete the duplicated one

``````b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
``````

of