Issue with PETScSENS solver in dolfinx for two field problem

Hello everyone,

I am trying to implement the phase field method (two field problem) using PETScSENS solver in dolfinx environment as done on the following link:

This is my code

from dolfinx import fem, plot, mesh, cpp , la           # dolfinx library
from mpi4py import MPI
import matplotlib.pyplot as plt    # Plotting Library
import numpy as np
from ufl import replace, TestFunction, TrialFunction, nabla_div, Measure, grad, inner, dot, derivative, lhs, rhs 

import dolfinx

import ufl

from mpi4py import MPI
from petsc4py import PETSc
import sys, os, shutil, math
import numpy as np
from math import pi
from dolfinx.mesh import locate_entities, meshtags
from dolfinx.fem import FunctionSpace, Function, locate_dofs_geometrical, dirichletbc, Constant
from dolfinx.fem.petsc import LinearProblem
from petsc4py.PETSc import ScalarType
import petsc4py
from snes_problem import SNESProblem
length = 200.  
num_elem = 400 
mesh = mesh.create_interval(MPI.COMM_WORLD,num_elem,[0., length])

V_u = FunctionSpace(mesh, ('Lagrange', 1)) 
V_alpha = FunctionSpace(mesh, ('Lagrange', 1) )

u, du, v = Function(V_u), TrialFunction(V_u), TestFunction(V_u) 
alpha, dalpha, beta = Function(V_alpha), TrialFunction(V_alpha), TestFunction(V_alpha)
lb = Function(V_alpha)
ub = Function(V_alpha)
u_R = Function(V_u)
disp = Function(V_u)
dx = Measure("dx",domain=mesh)
ds = Measure("ds",domain=mesh)

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

dof_left_u = locate_dofs_geometrical(V_u, on_boundary_L)
bcl = dirichletbc(ScalarType(0), dof_left_u, V_u)

E = Constant(mesh,30000.0)
nu = Constant(mesh,0.2)
lmbda = (E*nu)/((1+nu)*(1-2*nu))
mu = E/(2*(1+nu))

ell = Constant(mesh,5.)
kapa = 1e-7 
Gc = Constant(mesh,0.12)
c_w = 2.0                                    

def w(alpha):
    return alpha**2

def a(alpha):
    k_ell = Constant(mesh,1.e-6)
    return (1-alpha)**2 +k_ell
def eps(u):
    """Strain tensor as a function of the displacement"""
    return grad(u)

def sigma_0(u):
    """Stress tensor of the undamaged material as a function of the displacement"""
    mu    = E / (2.0*(1.0 + nu))
    lmbda = (E * nu) / ((1.0 - 2.0*nu)*(1.0 + nu))
    return E*(eps(u))

def sigma(u,alpha):
    """Stress tensor of the damaged material as a function of the displacement and the damage"""
    return (a(alpha))*sigma_0(u)

# Apply the limit on the phase field
zero_alpha = Function(V_alpha)
with zero_alpha.vector.localForm() as bc_local:

one_alpha = Function(V_alpha)
with one_alpha.vector.localForm() as bc_local:


# Governing Equations
f_ext = Constant(mesh,0.0)
external_work = dot(f_ext, u) * dx 
elastic_energy = 0.5*inner(sigma(u,alpha), eps(u))*dx 
dissipated_energy = Gc/float(c_w)*(w(alpha)/ell + ell*dot(grad(alpha), grad(alpha)))*dx 
total_energy = elastic_energy + dissipated_energy - external_work

E_u = derivative(total_energy,u,v)

E_alpha = derivative(total_energy,alpha,beta)
E_alpha_alpha = derivative(E_alpha,alpha,dalpha)

bc_alpha_L = dirichletbc(Constant(mesh,0.0), dof_left_u, V_alpha)
bcs_alpha = [bc_alpha_L]
dolfinx.fem.set_bc(ub.vector, bcs_alpha)
damage_problem = SNESProblem(E_alpha, E_alpha_alpha, alpha, bcs_alpha)

b = la.create_petsc_vector(V_alpha.dofmap.index_map, V_alpha.dofmap.index_map_bs)
J = dolfinx.fem.petsc.create_matrix(damage_problem.a)

solver_alpha_snes = PETSc.SNES().create()
solver_alpha_snes.setFunction(damage_problem.F, b)
solver_alpha_snes.setJacobian(damage_problem.J, J)
solver_alpha_snes.setTolerances(rtol=1.0e-9, max_it=50)


solver_alpha_snes.solve(None, alpha.vector) 

code for module sens_problem is

from typing import List

from petsc4py import PETSc

import dolfinx
import ufl
class SNESProblem:
    """Nonlinear problem class compatible with PETSC.SNES solver.

    def __init__(self, F: ufl.form.Form, J: ufl.form.Form, u: dolfinx.fem.Function, bcs: List[dolfinx.fem.dirichletbc]):
        """This class set up structures for solving a non-linear problem using Newton's method.

        F: Residual.
        J: Jacobian.
        u: Solution.
        bcs: Dirichlet boundary conditions.
        F, J = dolfinx.fem.form(F), dolfinx.fem.form(J)
        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 = dolfinx.fem.petsc.create_matrix(self.a)
        self.b = dolfinx.fem.petsc.create_vector(self.L)

    def F(self, snes: PETSc.SNES, x: PETSc.Vec, b: PETSc.Vec):
        """Assemble the residual F into the vector b. 

        snes: the snes object
        x: Vector containing the latest solution.
        b: Vector to assemble the residual into.
        # We need to assign the vector to the function

        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        # Zero the residual vector
        with b.localForm() as b_local:
        dolfinx.fem.assemble_vector(b, self.L)
        # Apply boundary conditions
        dolfinx.fem.apply_lifting(b, [self.a], [self.bcs], [x], -1.0)
        dolfinx.fem.set_bc(b, self.bcs, x, -1.0)
    def J(self, snes, x: PETSc.Vec, A: PETSc.Mat, P: PETSc.Mat):
        """Assemble the Jacobian matrix.

        x: Vector containing the latest solution.
        A: Matrix to assemble the Jacobian into.
        dolfinx.fem.assemble_matrix(A, self.a, self.bcs)

when I run the code, I got the following error. I am not able to understand this issue and how to fix it.

python3: src/petsc4py.PETSc.c:319281: __Pyx_PyCFunction_FastCall: Assertion `!PyErr_Occurred()' failed.
$ *** Process received signal ***
$ Signal: Aborted (6)
$ Signal code:  (-6)
$ [ 0] /lib/x86_64-linux-gnu/[0x7f7b3e936090]
$ [ 1] /lib/x86_64-linux-gnu/[0x7f7b3e93600b]
$ [ 2] /lib/x86_64-linux-gnu/[0x7f7b3e915859]
$ [ 3] /lib/x86_64-linux-gnu/[0x7f7b3e915729]
$ [ 4] /lib/x86_64-linux-gnu/[0x7f7b3e926fd6]
$ [ 5] /usr/lib/python3/dist-packages/petsc4py/lib/[0x7f7b2f9ce9f6]
$ [ 6] /usr/lib/python3/dist-packages/petsc4py/lib/[0x7f7b2fc559cd]
$ [ 7] /lib/x86_64-linux-gnu/[0x7f7b3c6d5430]
$ [ 8] /lib/x86_64-linux-gnu/[0x7f7b3d2479be]
$ [ 9] /lib/x86_64-linux-gnu/[0x7f7b3d1ed3ee]
$ [10] /lib/x86_64-linux-gnu/[0x7f7b3d2559b9]
$ [11] /usr/lib/python3/dist-packages/petsc4py/lib/[0x7f7b2fa127d8]
$ [12] python3[0x504d49]
$ [13] python3(_PyEval_EvalFrameDefault+0x85a)[0x56bbfa]
$ [14] python3(_PyEval_EvalCodeWithName+0x26a)[0x569dba]
$ [15] python3(PyEval_EvalCode+0x27)[0x6902a7]
$ [16] python3[0x67f951]
$ [17] python3[0x67f9cf]
$ [18] python3[0x67fa71]
$ [19] python3(PyRun_SimpleFileExFlags+0x197)[0x681b97]
$ [20] python3(Py_RunMain+0x212)[0x6b9d32]
$ [21] python3(Py_BytesMain+0x2d)[0x6ba0bd]
$ [22] /lib/x86_64-linux-gnu/[0x7f7b3e917083]
$ [23] python3(_start+0x2e)[0x5fc5fe]
$ *** End of error message ***
Aborted (core dumped)

Can any anyone please help me to fix the issue.

Thank you

A subtle issue throughout this code is that you are using
and not
This should be changed for all functions that interfaces with petsc objects, i.e, assemble_*, apply lifting, set_bc etc.

1 Like

Dear Prof. Dokken,

Thank you for your prompt response and time. I am able to solve this issue with your solution.