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:
https://newfrac.gitlab.io/newfrac-fenicsx-training/04-phase-field/phase-field.html
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:
bc_local.set(0.0)
one_alpha = Function(V_alpha)
with one_alpha.vector.localForm() as bc_local:
bc_local.set(1.0)
ub.interpolate(one_alpha)
# 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.setType("vinewtonrsls")
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.getKSP().setType("preonly")
solver_alpha_snes.getKSP().setTolerances(rtol=1.0e-9)
solver_alpha_snes.getKSP().getPC().setType("lu")
solver_alpha_snes.setVariableBounds(lb.vector,ub.vector)
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.
Parameters
==========
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.
Parameters
==========
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)
x.copy(self.u.vector)
self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
# Zero the residual vector
with b.localForm() as b_local:
b_local.set(0.0)
dolfinx.fem.assemble_vector(b, self.L)
# Apply boundary conditions
dolfinx.fem.apply_lifting(b, [self.a], [self.bcs], [x], -1.0)
b.ghostUpdate(addv=PETSc.InsertMode.ADD,mode=PETSc.ScatterMode.REVERSE)
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.
Parameters
==========
x: Vector containing the latest solution.
A: Matrix to assemble the Jacobian into.
"""
A.zeroEntries()
dolfinx.fem.assemble_matrix(A, self.a, self.bcs)
A.assemble()
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/libc.so.6(+0x43090)[0x7f7b3e936090]
$ [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0xcb)[0x7f7b3e93600b]
$ [ 2] /lib/x86_64-linux-gnu/libc.so.6(abort+0x12b)[0x7f7b3e915859]
$ [ 3] /lib/x86_64-linux-gnu/libc.so.6(+0x22729)[0x7f7b3e915729]
$ [ 4] /lib/x86_64-linux-gnu/libc.so.6(+0x33fd6)[0x7f7b3e926fd6]
$ [ 5] /usr/lib/python3/dist-packages/petsc4py/lib/PETSc.cpython-38-x86_64-linux-gnu.so(+0xa39f6)[0x7f7b2f9ce9f6]
$ [ 6] /usr/lib/python3/dist-packages/petsc4py/lib/PETSc.cpython-38-x86_64-linux-gnu.so(+0x32a9cd)[0x7f7b2fc559cd]
$ [ 7] /lib/x86_64-linux-gnu/libpetsc_real.so.3.12(PetscError+0x140)[0x7f7b3c6d5430]
$ [ 8] /lib/x86_64-linux-gnu/libpetsc_real.so.3.12(SNESComputeFunction+0x35e)[0x7f7b3d2479be]
$ [ 9] /lib/x86_64-linux-gnu/libpetsc_real.so.3.12(SNESSolve_VINEWTONRSLS+0x15e)[0x7f7b3d1ed3ee]
$ [10] /lib/x86_64-linux-gnu/libpetsc_real.so.3.12(SNESSolve+0x849)[0x7f7b3d2559b9]
$ [11] /usr/lib/python3/dist-packages/petsc4py/lib/PETSc.cpython-38-x86_64-linux-gnu.so(+0xe77d8)[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/libc.so.6(__libc_start_main+0xf3)[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
Manish