ArityMismatch error

General
1

I am a beginner user of fenics. I am writing a code for phase field ductile fracture. The following is my code. I have removed the part for solution of damage variable in the staggered algorithm, and the function for return mapping algorithm.

# PHASE FIELD DUCTILE FRACTURE
# Von mises J2 plasticity
#------------------------------------------------
from dolfin import *
from matplotlib import pyplot
import numpy as np

#------------------------------------------------

#                  MESH

#------------------------------------------------
# Load mesh and plot
mesh = Mesh('mesh.xml')
plot(mesh, title="Finite element mesh")
pyplot.show()
#------------------------------------------------

#   FUNCTION SPACE, FUNCTIONS & SOLUTION VARIABLES 

#------------------------------------------------

Vp = FunctionSpace(mesh, 'CG', 1)


# DISPLACEMENT - u
# For displacement
W = VectorFunctionSpace(mesh, 'CG', 1)
v = TrialFunction(W)
u_ = TestFunction(W)
du = Function(W, name="iteration correction")

# For H+: Elastic strain energy
# Energy governing damage
WW = FunctionSpace(mesh, 'DG', 0)

# Solution space variables for u_{t} and u_{t-1} = unew, u_old
unew, uold = Function(W), Function(W)

# PLASTICITY: FOR HISTORY, FUNTIONS
WT = TensorFunctionSpace(mesh, 'CG', 1, shape = (3,3))
sig = Function(WT)
epsilon_p = Function(WT)
n_elas_p = Function(WT)
alpha = Function(Vp, name="Cumulative plastic strain")
delta_lmbda_p = Function(Vp, name="Plastic multiplier")
beta_1p = Function(Vp)
beta_2p = Function(Vp)
#------------------------------------------------

#  MATERIAL PROPERTIES & DEGD FUNCTION PARAMETERS

#------------------------------------------------
# Gc: Critical energy release rate: Mpa mm
Gc = 9.31

# Lames parameters: MPa
lmbda = 53.473e3
mu = 27.280e3

# Bulk modulus
K_n = lmbda + 2*mu/3

# Yield stress: MPa
sigma_y = 345

# Hardening modulus: MPa
H_mod = 500

# Phase field parameters
# l :mm
l = 0.02
eta = 1e-5
#------------------------------------------------

#   CONSTITUTIVE FUNCTIONS

#------------------------------------------------

# Macauly bracket
# Function required for def of stress
def mac_bracket_plus(a):
    return 0.5*(a + abs(a))

def mac_bracket_minus(a):
    return 0.5*(a - abs(a))

# Strain 
def epsilon(u):
    e = sym(grad(u))
    return as_tensor([[e[0, 0], e[0, 1], 0],
                      [e[0, 1], e[1, 1], 0],
                      [0, 0, 0]])

def grad3(u):
    k = grad(u)
    return as_tensor([[k[0, 0], k[0, 1], 0],
                      [k[0, 1], k[1, 1], 0],
                      [0, 0, 0]])


def g_degd(p):
    return (1-p)**2 + eta


# Stress elastic
def sigma_el(u):
    ep = epsilon(u)
    return g_degd(pold) * (K_n * mac_bracket_plus(tr(ep)) * Identity(3) + 2 * mu * dev(ep)) + K_n * mac_bracket_minus(tr(ep)) * Identity(3)
    


def psi_plus(u):
    return 0.5*(K_n)*(mac_bracket_plus(tr(epsilon(u))))**2+\
           mu*inner(dev(epsilon(u)),dev(epsilon(u)))

# Irreversibility condition: using conditional operator
# Function returns the maximum of tensile part of strain energy density
def H(uold,unew,Hold):
    return conditional(lt(psi_plus(uold),psi_plus(unew)),psi_plus(unew),Hold)


def sigma_tangent(uu):
    B = grad3(uu)
   # print(epsilon(dui).ufl_shape)
    return (K_n/2)* (g_degd(pold) *(B + B* tr(B)/abs(tr(B))) + (B -B * tr(B)/abs(tr(B)))) +\
            (g_degd(pold) * 2 * mu * ( beta_2p  * delta_lmbda_p - beta_1p)) * inner( n_elas_p, B) * n_elas_p +\
                   g_degd(pold) * 2 * mu * (1- beta_2p * delta_lmbda_p) * dev(B) 


# -----------------------------------------------

#  BOUNDARY CONDITIONS, NORMAL

# -----------------------------------------------

# top boundary
top = CompiledSubDomain("near(x[1], 1.0) && on_boundary")

# bottom boundary
bot = CompiledSubDomain("near(x[1], 0) && on_boundary")

# crack location
def Crack(x):
    return ((x[1] > 0.499) & (x[1] < 0.501)) and x[0] <= 0.5

# load, defined using user expression, t
load = Expression("t", t = 0.0, degree=1)

# bc, bottom
bcbot = bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)

# load applied on top u_y
# boundary condition on top
bctop = DirichletBC(W.sub(1), load, top)

# grouping disp bc together
bc_u = [bcbot, bctop]



# MeshFunction: identify the boundaries of the domain
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

# initialize boundaries before marking
boundaries.set_all(0)

# mark the top boundary
top.mark(boundaries, 1)

# to integrate over boundari 
ds = Measure("ds")(subdomain_data= boundaries)

# define normal to the surface/ boundaries
n = FacetNormal(mesh)

# -----------------------------------------------

#    VARIATIONAL FORMULATION 

# -----------------------------------------------

# Implementing Newton method manually
# K matrix
# recall the problem is
# K deltaU = f
# we will solve for deltaU and update U
a_Newton = inner(grad3(v), sigma_tangent(u_)) * dx
res = -inner( grad3(u_), sig) * dx

# -----------------------------------------------

#                   LOAD

# -----------------------------------------------

# specified as remote displacement
u_r = 0.01
# -----------------------------------------------

#           ITERATIVE SCHEME

# -----------------------------------------------

# Newton - Raphson procedure parameters
Nitermax, tol = 200, 1e-4

# time step
deltaT = 0.01
t = 0
# Staggerd solution strategy

while t<=1.0:
    t += deltaT
    if t>=0.5:
        deltaT = 0.005
    load.t = t*u_r
    iter = 0    # iteration counter
    err = 1
    A, Res = assemble_system(a_Newton, res, bc_u)
    nRes0 = Res.norm("l2")
    nRes = nRes0
    unew.assign(uold)
    while nRes/nRes0 > tol and niter < Nitermax:
        solve(A, du.vector(), Res, "mumps")
        unew.assign(unew + du)
        epsilon_p_, alpha_, sig_, n_elas_p_, beta_1p_, beta_2p_, delta_lmbda_ = return_map(unew, epsilon_p, alpha, pold)
        sig = project(sig_, WT)
        beta_1p = project(beta_1p_, Vp)
        beta_2p = project(beta_2p_, Vp)
        n_elas_p = project(n_elas-p_, WT)
        delta_lmbda_p = project(delta_lmbda_, Vp)
        A, Res = assemble_system(a_Newton, res, bc)
        nRes = Res.norm("l2")
        iter += 1
    while err > tol:
        solver_phi.solve()
        err = errornorm(pnew,pold,norm_type = 'l2',mesh = None)
    uold.assign(unew)           
    pold.assign(pnew)
    Hold.assign(project(psi(unew), WW))
    alpha = project(alpha_, Vp)

While running the code, I am getting the following error. Can anyone help me with this?
At this stage, I have not used quadrature elements; could that cause any issue?

Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
  File "/home/anupama/OneDrive/WORK_2023/PFM_nonlinearJUNE2023 (copy)/Ductile_new.py", line 286, in <module>
    A, Res = assemble_system(a_Newton, res, bc_u)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 408, in assemble_system
    A_dolfin_form = _create_dolfin_form(A_form, form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 60, in _create_dolfin_form
    return Form(form,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py", line 43, in __init__
    ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 50, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 100, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 130, in jit_build
    module, signature = dijitso.jit(jitable=ufl_object,
  File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 65, in jit_generate
    code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 142, in compile_form
    return compile_ufl_objects(forms, "form", object_names,
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in analyze_ufl_objects
    form_datas = tuple(_analyze_form(form, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in <genexpr>
    form_datas = tuple(_analyze_form(form, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 169, in _analyze_form
    form_data = compute_form_data(form,
  File "/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py", line 407, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is
  File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 177, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 159, in check_integrand_arity
    arg_tuples = map_expr_dag(rules, expr, compress=False)
  File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 36, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress,
  File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 97, in map_expr_dags
    r = handlers[v._ufl_typecode_](v)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 41, in nonlinear_operator
    raise ArityMismatch("Applying nonlinear operator {0} to expression depending on form argument {1}.".format(o._ufl_class_.__name__, t))
ufl.algorithms.check_arities.ArityMismatch: Applying nonlinear operator Abs to expression depending on form argument v_0.
2

You are solving a non-linear problem, and therefore your variational form should not be split into a and L as you are trying to do. You should write your variational for with u_ being a function, and write it on the form F==0.
If you then instead want to use a Newton method, you can form a_newton by calling derivative(F, u_, du), where du is a TrialFunction.

1 Like
3

Thank you very much for the reply.
Can I specify Jacobian manually and then use the Newton method?

4

Yes, you can. You just need to make sure that the Jacobian is linear with respect to u and v (it should be a bilinear form, forming a matrix. See for instance
https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html

1 Like
5

Thank you. I was able to resolve this issue. While running the same code, I am now getting the following error message. Can anyone suggest what the error could be

Traceback (most recent call last):
  File "/home/anupama/OneDrive/WORK_2023/ductile_july/Ductile_new.py", line 298, in <module>
    solve(A, du.vector(), Res, "mumps")
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 240, in solve
    return dolfin.la.solver.solve(*args)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/la/solver.py", line 72, in solve
    return cpp.la.solve(A, x, b, method, preconditioner)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------
6

Without producing the code that you are running, its hard to give any guidance.

7

I apologize for the confusion.
I am writing a code for phase field ductile fracture. For the plasticity part, I have followed the following fenics tutorial Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation
My code is given below.(I have removed unnecessary parts of the code) and this time I made sure jacobian is linear w.r.t test and trial functions. To make the code simple, I have only incorporated the solution of mechanical problem. ( I have omitted the solution procedure for variable phi)

mesh = Mesh('mesh1.xml')
# DAMAGE VARIABLE - PHI
# Function spaces for the problem
V = FunctionSpace(mesh, 'CG', 1)
Ve = FiniteElement("Quadrature", mesh.ufl_cell(), degree=1, quad_scheme='default')
Vp = FunctionSpace(mesh, 'CG', 1)
p, q = TrialFunction(V), TestFunction(V)

# DISPLACEMENT - u
# For displacement

W = VectorFunctionSpace(mesh, 'CG', 1)
v = TrialFunction(W)
u_ = TestFunction(W)
du = Function(W, name="iteration correction")

# For H+: Elastic strain energy
# Energy governing damage
WW = FunctionSpace(mesh, 'DG', 0)

# Solution space variables for phi_{t}, phi_{t-1} and H_{t-1} = pnew, pold, Hold
pnew, pold, Hold = Function(V), Function(V), Function(WW)

# Solution space variables for u_{t} and u_{t-1} = unew, u_old
unew, uold = Function(W), Function(W)
# PLASTICITY: FOR HISTORY, FUNTIONS

WT = TensorFunctionSpace(mesh, 'CG', 1, shape = (3,3))
sig = Function(WT)
epsilon_p = Function(WT)
n_elas_p = Function(WT)
alpha = Function(Vp, name="Cumulative plastic strain")
delta_lmbda_p = Function(Vp, name="Plastic multiplier")
beta_1p = Function(Vp)
beta_2p = Function(Vp)
# -----------------------------------------------

#  BOUNDARY CONDITIONS, NORMAL

# -----------------------------------------------

# top boundary
top = CompiledSubDomain("near(x[1], 0.5) && on_boundary")

# bottom boundary
bot = CompiledSubDomain("near(x[1], -0.5) && on_boundary")

# crack location
def Crack(x):
    return abs(x[1]) < 1e-03 and x[0] <= 0.0

# load, defined using user expression, t
load = Expression("t", t = 0.0, degree=1)

# bc, bottom
bcbot= DirichletBC(W, Constant((0.0,0.0,0.0)), bot)

# load applied on top u_y
# boundary condition on top
bctop = DirichletBC(W.sub(1), load, top)

# grouping disp bc together
bc_u = [bcbot, bctop]

# phi, bc
bc_phi = [DirichletBC(V, Constant(1.0), Crack)]
# MeshFunction: identify the boundaries of the domain
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

# initialize boundaries before marking
boundaries.set_all(0)

# mark the top boundary
top.mark(boundaries, 1)

# to integrate over boundari 
ds = Measure("ds")(subdomain_data= boundaries)

# define normal to the surface/ boundaries
n = FacetNormal(mesh)

metadata = {"quadrature_degree": 1, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
# -----------------------------------------------
#    VARIATIONAL FORMULATION 

# -----------------------------------------------

# Implementing Newton method manually

# we will solve for deltaU and update U
a_Newton = inner(grad3(v), sigma_tangent(unew, u_)) * dxm
res = -inner( epsilon(u_), sig) * dxm
# Newton - Raphson procedure parameters
Nitermax, tol = 20, 1e-8
pold = 0.001
# time step
deltaT = 0.01
t = 0.01
# Staggerd solution strategy
while t<=1.0:
    t += deltaT
    if t>=0.5:
        deltaT = 0.001
    load.t = t*u_r
    niter = 0
    err1 = 1
    A, Res = assemble_system(a_Newton, res, bc_u)
    unew.assign(uold)
    while niter < Nitermax and err1 > tol:
        solve(A, du.vector(), Res, "mumps")
        unew.assign(unew + du)
        epsilon_p_, alpha_, sig_, n_elas_p_, beta_1p_, beta_2p_, delta_lmbda_ = return_map(unew, epsilon_p, alpha, pold)
        sig = project(sig_, WT)
        beta_1p = project(beta_1p_, Vp)
        beta_2p = project(beta_2p_, Vp)
        n_elas_p = project(n_elas_p_, WT)
        alpha = project(alpha_, Vp)
        delta_lmbda_p = project(delta_lmbda_, Vp)
        A, Res = assemble_system(a_Newton, res, bc_u)
        nRes = Res.norm("l2")
        print("iteration number",niter)
        niter += 1
    err1 = errornorm(unew,uold,norm_type = 'l2',mesh = None)
    uold.assign(unew)
    conc_p << alpha
    fname.close()
end = time.time()
print('simulation time',end-start)
8

Hi @anupama, I am quite impressed by your attempt to model ductile fracture using phase field approach, I was also attempting the same challenge a few weeks ago, I assume that for phase field modelling you are adapting code from Phase field fracture implementation in FEniCS and as you have mentioned for plasticity you are following Elasto-plastic analysis of a 2D von Mises material. I want to pose some questions regarding this model -

In the paper, 2 linear problems are being solved, one is for displacement and another is for damage variable phi which are -

image
image730×165 16.3 KB

When introducing plasticity the problem that we will be solving now will be nonlinear for displacement and linear for damage variable phi.

  1. My question is you are solving only nonlinear problem and skipping the other linear problem that is needed to evaluate phi (damage variable), will not that be inconsistent with the concept you want to implement the model on ?

  2. How will you introduce a dependency of phi on the displacement like done in the first equation by multiplying image in the case when we are solving a linear problem for displacement that is the case of linear elasticity ?

I hope that you will consider my questions because at the conceptual level, how we formulate variational forms in the case when plasticity is introduced because it doesn’t remains trivial because the nonlinearity comes into play.

Thanks.

9

Hi,
I have omitted the part regarding the solution of the damage variable here. One should solve for both phi and u , in a staggered/monolithic manner.

https://www.researchgate.net/publication/338897618_Phase_field_fracture_implementation_in_FEniCS , you cannot directly use this code. They have used the hybrid formulation for brittle fracture. So, you should reformulate your variational problem according to your need.
One way is to formulate your problem as in , Phase-field modeling of ductile fracture | SpringerLink,

1 Like
10

Thanks for the insightful take. I appreciate it.

Kindly provide your mesh file (in this case .xml file) with a minimum working example so that I can reproduce the error and could possibly give any suggestion.

According to @dokken in this post PETSc Krylov Solver RuntimeError with fine mesh , this information means that the residual is zero at the first iteration, which usually indicates that the initial guess solves your problem. This was also reported for another case at the old Q&A forum.

Kindly go through PETSc Krylov Solver RuntimeError with fine mesh and PETSc Krylov solver failed to converge with big neumann boundary conditions posts where similar error is addressed, from my understanding, you should try few things like keeping the mesh coarser everywhere else except in the vicinity of crack tip and its expected crack path, try using umfpack or lu instead or mumps and also consider using iterative solvers with preconditioners, one of the setup which works for me is Bitbucket.
I hope this helps.

1 Like