Error using PETSc SNES

General
1

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.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.vector)
        self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        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)
        b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        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))
        F = ufl.variable(I + ufl.grad(u))
        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))
    F = ufl.variable(I + ufl.grad(u)) 
    det_F  = ufl.variable(ufl.det(F))     
        
    metadata = {"quadrature_degree": 4}
    ds = ufl.Measure('ds', domain=domain, metadata=metadata)
    dx = ufl.Measure("dx", domain=domain, metadata=metadata)
    
    #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()
2

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.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.vector)
        self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        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)
        b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        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()
3

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.

4

See for instance: oasisx/src/oasisx/fracstep.py at 3e7f8ea2c229bfdd1d8d73901d94aad9472477fb · ComputationalPhysiology/oasisx · GitHub

It would be insert and forward.

5

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])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)

As = assemble_matrix(as)
As.assemble()
bs = assemble_vector(Ls)
bs.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

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

delb = b.copy()
delb.zeroEntries()
I.multTranspose(bs, delb)
delb.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
b += delb
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

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])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
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)
bs.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
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)
6

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