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 = (Enu)/((1+nu)(1-2nu))
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.0nu)(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.5inner(sigma(u,alpha), eps(u))dx
dissipated_energy = Gc/float(c_w)(w(alpha)/ell + elldot(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)
disp_a = Constant(mesh,0.01)
def on_boundary(x):
return np.isclose(x[0],length)
dof_right_u = locate_dofs_geometrical(V_u, on_boundary)
bcr = dirichletbc(disp_a, dof_right_u, V_u)
bc_disp = [bcl, bcr]
E_du = replace(E_u,{u:du})
problem_u = LinearProblem(a = lhs(E_du), L = rhs(E_du), bcs = bc_disp,u = u, petsc_options={“ksp_type”: “preonly”, “pc_type”: “lu”, “pc_factor_mat_solver_type”: “umfpack”})
solver_u = problem_u.solve()
bc_alpha_L = dirichletbc(Constant(mesh,0.0), dof_left_u, V_alpha)
bc_alpha_R = dirichletbc(Constant(mesh,0.0), dof_right_u, V_alpha)
bcs_alpha = [bc_alpha_L,bc_alpha_R]
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)
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